Solution of the statistical equations

The system Eq. 27 for $[n_{\rm new}]$ is non-linear with respect to the $n_{i,\rm new}$ and n_e because the coefficients of the $\ifmmode{[R_{ij}^*]}\else\hbox{$[R_{ij}^*]$}\fi$ and $\ifmmode{[R_{ji}^*]}\else\hbox{$[R_{ji}^*]$}\fi$-operators are quadratic in $n_{i,\rm new}$and the dependence of the Saha-Boltzmann factors and the collisional rates on the electron density, respectively. The system is closed by the abundance and charge conservation equations. To simplify the iteration scheme, and to take advantage of the fact that not all levels strongly influence all radiative transitions, we use a linearized and splitted iteration scheme for the solution of Eq. 27. This scheme has the further advantage that many different elements in different ionization stages and even molecules can be treated consistently. A problem where this is important is the modeling of nova and supernova atmospheres, where there are typically very large temperature gradients within the line forming region of the atmosphere.

To linearize Eq. 27, we follow [33] and replace terms of the form $n_{j,\rm new} [R_{ji}^{*}][n_{\rm new}]$ in Eq. 27 by $n_{j,\rm old} [R_{ji}^{*}][n_{\rm new}]$:

This removes the major part of the non-linearity of Eq. 27 but the modified system is still non-linear with respect to n_e and still has the high dimensionality of the original system. However, as has been noted before, not all levels are strongly coupled to all other levels. Equation 29 can be solved for each element (or groups of elements if they are coupled tightly) separately if the electron density is given. Therefore, we split the electron density calculation from the rate equation solution so that the n_e can be considered as given during the rate equation solution process and changes in the electron density are then accounted for in an outer iteration to find a consistent solution of the rate equations and the electron densities.

The most important advantage of this method is that it requires the solution of large linear systems and low-dimensional non-linear system (for the electron density). Thus, its solution is more stable and uses much less computer resources (time and memory) than the direct solution of the original non-linear equations. This allows us to treat many more levels with this method then with more conventional algorithms. Using a nested iteration scheme like the one described here will slow down the convergence of the iterations, but this is more than offset for by the possibility of calculating much larger models with less memory. Since we are able to solve a separate equation for each group of elements, we can trivially parallelize the solution by distributing the groups among the available processors. Convergence acceleration methods can in principle be used, but they frequently lead to convergence instabilities in the nested iterations for the solution of the statistical equilibrium equations.

We have so far assumed that the electron density n_e is given. Although this is a good assumption if only trace elements are considered, the electron density may be sensitive to non-LTE effects. This can be taken into account by using either a fixed point iteration scheme for the electron density or, if many species or molecules are included in the non-LTE equation of state, by a modification of the LTE partition functions to include the effects of non-LTE in the ionization equilibrium. The latter method replaces the partition function, $Q=\sum g_i\exp(-E_i/kT)$, with its non-LTE generalization, $Q_{\rm NLTE} = \sum b_i g_i\exp(-E_i/kT)$, and uses $Q_{\rm NLTE}$ in the solution of the ionization/dissociation equilibrium equation. We use this method because of the large number of elements with various ionization stages as well as molecules and condensation of dust grains included in statistical equilibrium calculations (and not all of them in non-LTE).

Our iteration scheme for the solution of the multi-level non-LTE problem can be summarized as follows: (1) for given n_i and n_e, solve the radiative transfer equation at each wavelength point and update the radiative rates and the approximate rate operator, (2) solve the linear system Eq. 29 for each group for a given electron density, (3) compute new electron densities (by either fixed point iteration or the generalized partition function method), (4) if the electron density has not converged to the prescribed accuracy, go back to step 2, otherwise go to step 1. The iterations are repeated until a prescribed accuracy for the n_e and the n_i is reached. It is important to account for coherent scattering processes during the solution of the wavelength dependent radiative transfer equation, it explicitly removes a global coupling from the iterations.