text:
# An iterative technique for solving equations of statistical equilibrium

## 1 Introduction

In order to extract all the information contained in an astronomical spectrum, the structure of the emitting layers must be modelled and the resulting spectrum computed. In doing this, various technical challenges must be overcome, including in particular that of solving the equations of statistical equilibrium for appropriately detailed atomic models of the relevant elements. Now, in some cases, quite simple atomic models suffice, as when predicting the strengths of collisionally-excited forbidden lines in the spectra of gaseous nebulae (e.g., de Robertis et al. 1987). But when treating line formation in stellar atmospheres and winds or in novae and supernovae ejecta, the important contibutions of the iron group elements Sc–Ni require atomic models with numbers of energy levels $`N10^3`$ per atom (e.g., Hubeny & Lanz 1995). The direct solution of the equations of statistical equilibrium with a standard linear equation-solving package is then decidedly problematic. As Press et al. (1992, p.23) note, solving systems with $`N`$ as large as several hundred is feasible with double precision arithmetic but beyond this the limiting factor is generally machine time. This limitation arises because the required computing time scales as $`N^3`$ (e.g., Press et al. 1992, p.34).

Classically, this computational problem was completely circumvented by simply assuming that the levels are in local thermodynamic equilibrium (LTE); and essentially the same strategy can in be used to approximate reality more closely if the Boltzmann–Saha LTE formulae are replaced by formulae that take account of the dominant NLTE effect in the problem of interest. This strategy has been followed in Monte Carlo codes for stellar winds (Abbott & Lucy 1985; Lucy & Abbott 1993) and supernovae (Lucy 1987; Mazzali & Lucy 1993; Lucy 1999) where the dominant NLTE effect is dilution of the radiation field due to spherical extension. But this approach, despite its theoretical basis, has been criticized by Hillier & Miller (1998) as giving an ionization structure and line source functions that are ad hoc or, at best, inconsistent with the radiation field. On the other hand, Springmann & Puls (1998) find an encouraging degree of agreement between the ionization formula and accurate NLTE calculations.

The solution strategy preferred by Hillier & Miller (1998) and indeed by most other investigators is that pioneered by Anderson (1989) in which a reduced atomic model is constructed by consolidating neighbouring energy levels into single superlevels. In this way, $`N`$ is reduced to $`50`$, so that solution with a standard package is possible even with single precision arithmetic (Press et al. 1992, p.23). Then, from the the reduced atom’s level populations, those for the complete atom – needed for the transfer calculation – are obtained by assuming that all the sublevels of a given superlevel have the same departure coefficient.

In the earliest implementations of Anderson’s seminal idea, the partitioning into superlevels was based solely on excitation energy (Anderson 1989, Dreizler & Werner 1993). But subsequently Hubeny & Lanz (1995) introduced the constraint that all sublevels of a given superlevel must be of the same parity. This eliminates the awkward problem of radiative transitions between a superlevel’s sublevels as well as enhancing the validity of the assumption of identical departure coefficients. Nevertheless, for some diagnostically-important emission lines, Hillier & Miller (1998) still find it necessary to solve for the populations of the individual emitting levels.

Evidently, current treatments of the statistical equilibria of complex atoms represent interim procedures pending more powerful computers or better solution strategies. An example of the latter is the subject of this paper.

The technique described and tested herein arose in an ongoing effort to extend the range and sophistication of problems treatable with Monte Carlo methods (Lucy 1999a,b). But since the technique itself is more widely applicable, it is presented here separated from its Monte Carlo origins.

## 2 An iterative scheme

In this section, an iterative scheme is described for determining the statistical equilibrium of a model atom with large N in a specified environment – i.e., the ambient radiation field is assumed known as are the electron density and temperature. Because of these assumptions, this scheme will often be used within an outer iteration loop determining the thermal structure of the emitting layers as well as the internal radiation field.

### 2.1 Basic idea

In the problems of primary interest, the ambient radiation field, although not black body, nevertheless derives from thermal processes. For such a non-pathological radiation field, the dominant ions of a given species have heavily populated low levels, especially the ground state and low-lying metastable levels, and sparsely populated high levels. In this circumstance, the populations of high levels are largely determined by their coupling to the low levels of the same ion via bound-bound (b-b) transitions and to the low levels of the upper ion via free-bound (f-b) transitions. Accordingly, to a first approximation, the population of a high level is independent of the populations of other high levels of the same or of other ions. This then suggests that the populations of high levels be estimated one by one with a simple relaxation scheme that brings such levels into statistical equilibrium with the readily estimated populations of low levels. Then, having thus derived approximate populations for high levels, we can obtain improved populations for the low levels with the same relaxation scheme.

The simple, fast iterative technique just described often yields accurate level populations for a single ion with far less computer time than required by the direct solution of the basic equations (Li, McCray & Sunyaev 1993). But in some circumstances, intial convergence is followed by divergence, and so the technique lacks the robustness needed if it is to be included in a code of wide applicability. Accordingly, this relaxation procedure is here incorporated into the superlevel technique discussed in Section 1. By combining these procedures, the desired robustness is achieved without sacrificing speed. Moreover, in contrast to the superlevel technique, there is no limit to the achievable accuracy.

### 2.2 Equations of statistical equilibrium

Level $`i`$ of the complete model atom is in statistical equilibrium if

$$R_\mathrm{}in_{\mathrm{}}(R_i\mathrm{}+R_{iu})n_i+R_{ui}n_u=0$$
(1)
where a summation convention is adopted for upper ($`u>i`$) and lower ($`\mathrm{}<i`$) levels. The quantity $`R_{ij}`$ denotes the total rate coefficient for the transition $`ij`$, and $`n_i`$ denotes the relative population of level $`i`$ – i.e., we impose the normalization constraint

$$\underset{i}{}n_i=1$$
(2)
This normalization constraint together with statistical equilibrium equations for $`N1`$ levels constitute a system of $`N`$ equations in the $`N`$ unknowns $`n_i`$. If the ambient conditions are given, the coefficients $`R_{ij}`$ are independent of $`n_i`$ and the system is linear.

Note that equation (1) applies to all levels of all ions included in the complete atomic model. Accordingly, $`R_{ij}`$ refers to b-b transitions when $`i`$ and $`j`$ denote levels of the same ion and to b-f or f-b transitions when they denote levels of different ions.

### 2.3 Superlevels

Levels consolidated into a superlevel will be referred to as sublevels of the parent superlevel, and the complete set of superlevels defines the reduced atomic model. Rate coefficients coupling superlevels determine the statistical equilibrium of the reduced atom exactly as do the $`R_{ij}`$ for the complete atom – see equation (1). Moreover, these rate coefficients must be defined so that, if the populations of their sublevels are known exactly, the statistical equilibrium of the reduced atom yields the exact populations for the superlevels. Clearly, this is a necessary condition if the iterative scheme is to converge to the exact solution for the complete atom.

#### 2.3.1 Coupling of superlevels

Let $`j=j_1,j_2,\mathrm{}`$ be an index denoting the sublevels of superlevel $`J`$ and let $`n_j`$ be the exact populations of these sublevels. The population of the superlevel is then naturally defined to be

$$n_J=\underset{j}{}n_j$$
(3)
from which it follows that

$$\underset{J}{}n_J=1$$
(4)
Now consider transitions between superlevels $`J`$ and $`K`$. The rate $`R_{JK}n_J`$ at which the hypothetical transition $`JK`$ transfers population from superlevel $`J`$ to superlevel $`K`$ must equal (for the exact solution) the sum of the rates $`R_{jk}n_j`$ at which the actual transitions $`jk`$ transfer population from the sublevels of $`J`$ to the sublevels of $`K`$. Thus

$$R_{JK}n_J=\underset{j}{}\underset{k}{}R_{jk}n_j$$
(5)
and correspondingly for the inverse coefficient $`R_{KJ}`$.

With rate coefficients for superlevels thus defined, the equations of statistical equilibrium for the complete model atom, when summed over the sublevels of individual superlevels, yield the corresponding equilibrium equations for the reduced atom. Accordingly, if $`n_i`$ is the exact solution vector for equations (1) and (2), then the corresponding superlevel vector $`n_J`$ from equation (4) is the exact solution for the reduced atom.

#### 2.3.2 Positivity

A complete set of approximate rate coefficients $`R_{JK}`$ for the reduced atom can be computed from equation (5) if the vector of population ratios $`r_j=n_j/n_J`$ is estimated for each superlevel $`J`$. Provided these estimates $`r_j`$ are $`>0`$ for all $`j`$, the coefficients $`R_{JK}`$ will be $`>0`$, and so the equations of statistical equilibrium for the reduced atom, having exactly the form of equations (1) and (2), will give $`n_J>0`$ for all J.

### 2.4 Iteration steps

The iterative scheme for solving the equations of statistical equilibrium for the complete atomic model now proceeds step by step as follows:

1) From the current estimate for the level populations $`n_i`$ of the complete atom, the population ratios $`r_j=n_j/n_J`$ for each superlevel’s sublevels are computed and used in equation (5) to derive rate coefficients for transitions between superlevels.

2) The equations of statistical equilibrium for the reduced atomic model are solved directly with a standard linear equation-solving package, thus obtaining improved populations $`n_J`$ for the superlevels.

3) The populations of each superlevel’s sublevels are rescaled by setting $`n_j=n_Jr_j`$, where $`n_J`$ are the updated values from step 2) and $`r_j`$ the ratios from step 1). With this rescaling, the current estimate $`n_i`$ for the complete atom now satisfies the normalization constraint.

4) The populations of each ion’s energy levels are adjusteded sequentially from highest to lowest according to the following relaxation scheme

$$n_i=\frac{R_\mathrm{}in_{\mathrm{}}+R_{ui}n_u}{R_i\mathrm{}+R_{iu}}$$
(6)
where, as in equation (1), summation over the indices $`\mathrm{}`$ and $`u`$ is implied.

5) The vector $`n_i`$ is renormalized.

6) The current estimate $`n_i`$ is accepted if the adopted convergence criterion is satisfied.

7) If the criterion is not satisfied, a further iteration is initiated by inputting $`n_i`$ at step 1).

### 2.5 Comments

The following comments are intended to clarify the working of the iterative scheme:

a) To start the iterations, an initial estimate for the level populations is required. An obvious choice is the Saha–Boltzmann LTE population, and this is how the numerical examples of Section 3 are initiated. But when the scheme is implemented within an outer iteration loop for the thermal structure, better initial estimates are available after the first outer iteration.

b) From the initial estimate, only the ratios $`r_j=n_j/n_J`$ are actually used. Thus the scheme may be regarded as iterating on these ratios.

c) Equation (6) gives the value of $`n_i`$ that satisfies equation (1) for level $`i`$, given the current estimates of the vectors $`n_{\mathrm{}}`$ and $`n_u`$.

d) The value of $`n_i`$ from equation (6) is used in $`n_u`$ when $`n_{i1}`$ is calculated. Thus, improved values are used immediately rather than waiting for the next iteration cycle.

e) If we estimate the level populations of an ion, we are typically fairly accurate for the well-populated low levels and progressively less so for higher levels. This consideration is behind the downward sweep of the relaxation step. The least accurately determined levels are then adjusted first, and so their improved populations are available for the summations $`R_{ui}n_u`$ when later adjusting the lower levels of the same ion.

f) Step 4) clearly preserves the positivity of $`n_i`$. Together with the remarks in Section 2.3.2, this proves that the iterative scheme as a whole always returns $`n_i>0`$ for all $`i`$ provided only that their initial values are $`>0`$.

g) Transitions between sublevels of the same parent superlevel do not appear in equation (5) and so make no direct contribution to the solution vector $`n_J`$ obtained at step 2). But all such transitions appear in equation (6). Accordingly, after the first iteration, these transitions contribute indirectly to $`n_J`$ as the relaxation step adjusts the ratios $`r_j`$ and hence also the rate coefficients $`R_{JK}`$.

h) In view of the physical meaning of both the numerator and denominator in equation (6), no problems with loss of accuracy due to the accumulation of roundoff errors will arise from the relaxation step even if the scheme is applied to atomic models with extremely large $`N`$. Accordingly, such problems should not arise for the iterative scheme provided the reduced atom has an appropriately limited $`N`$.

i) If the scheme converges, the particular choice of partitioning into superlevels becomes irrelevant. In contrast, with the standard superlevel technique, the partioning determines and limits the accuracy attained.

j) Steps 4) and 5) may be repeated several times before proceeding to step 6).

### 2.6 Convergence

If the iterative scheme is initiated with (or reaches) the exact solution $`n_i^{(x)}`$ of equations (1) and (2), then, as discussed in Section 2.3.1, steps 1) and 2) return the exact populations of the superlevels. Step 3) then recreates $`n_i^{(x)}`$, and this is left unchanged by step 4), since every level is found to be already in statistical equilibrium. On the other hand, any vector $`n_in_i^{(x)}`$, and which therefore does not satisfy equations (1) and (2), will be changed at step 4). Accordingly, the iterative scheme converges if and only if it reaches the exact statistical equilibrium solution for the complete atomic model.

Note that if the rate coefficient $`R_{ij}`$ is changed for the transition between any pair of sublevels belonging to the same superlevel, the solution for the complete atom changes but that for the reduced atom may not – cf. comment g) of Section 2.5. Accordingly, there may exist vectors $`n_in_i^{(x)}`$ that are unchanged by steps 1) to 3). Equation (6) is therefore fundamental to the ‘if and only if’ assertion in the above convergence statement.

For practical applications, the rate of convergence to $`n_i^{(x)}`$ is of prime concern. Unfortunately, no general quantitative results have been established on this question. But in the limiting case where each superlevel comprises just one sublevel, the iterative scheme reduces to the direct solution of equations (1) and (2) and so converges at the first iteration. As we increasingly depart from this limit by consolidating more and more levels into fewer and fewer superlevels, the rate of convergence will surely slow and, eventually, the scheme might even diverge. Fortunately, numerical experience indicates that, with the degree of consolidation typical of current implementations of the superlevel technique, the scheme does converge. Moreover, the rate of convergence is such that solutions of the desired accuracy are obtained with a huge reduction of computer time relative to that required by direct solution.

## 3 Numerical experiments

In this section, the iterative scheme described in Section 2 is applied to model atoms having the completeness of level structure needed for modern spectral synthesis calculations. No particular attempt has been made to compile the best modern data for these atoms since the intent is merely to illustrate the typical performance of the scheme.

### 3.1 Convergence criteria

In practical applications of iterative schemes, the exact solution is not known and so convergence is often judged from the corrections to the previous iterate. For atoms with large $`N`$, it is convenient to summarize the corrections to all the levels in a single quantity. Here we weight all levels equally and compute the mean of the absolute fractional changes to the individual levels’ populations. Thus the adopted measure of the correction at the $`r`$th iteration is

$$\delta ^{(r)}=\frac{1}{N}\underset{i}{}\frac{|n_i^{(r)}n_i^{(r1)}|}{n_i^{(r1)}}$$
(7)
By now, a well-recognized danger in judging convergence on the basis of small iterative corrections is that an acceptance criterion of the form $`\delta ^{(r)}<\delta `$ does not necessarily imply that the typical fractional error is $`\delta `$. The remaining errors will in fact be much larger if the scheme converges slowly. Accordingly, in the following numerical experiments, for which the exact solutions $`n_i^{(x)}`$ are known, the quantity $`\delta ^{(r)}`$ is compared to

$$ϵ^{(r)}=\frac{1}{N}\underset{i}{}\frac{|n_i^{(r)}n_i^{(x)}|}{n_i^{(x)}}$$
(8)
In this way, we may hope to establish values of $`\delta `$ required to achieve a specified absolute accuracy.

### 3.2 Fe II

As a first test, the scheme is applied to the 394-level model Fe II ion used previously (Lucy 1999b) to investigate treatments of line emissivity in a Monte Carlo spectral synthesis code for supernovae. Indeed, it was these emissivity experiments that led to the relaxation scheme that forms step 4) of the present technique.

The energy levels of the Fe II ion and the f-values for permitted transitions were extracted from the Kurucz–Bell (1995) compilation by M.Lennon (Munich). Einstein A-values for forbidden transitions are from Quinet et al. (1996) and Nussbaumer & Storey (1988). Collision strengths are from Zhang & Pradhan (1995), with the van Regemorter (1962) formula as default for permitted transtions and the value 0.05 as default for forbidden transitions.

As previously (Lucy 1999b), the ambient radiation field is taken to be $`WB_\nu (T_b)`$ with $`T_b=12500K`$ and dilution factor $`W=0.067`$ corresponding to the point $`r=2R`$ above the photosphere. But now instead of the low density limit, electron excitation is included with $`T_e=10000K`$ and $`N_e=10^7cm^3`$.

In order to solve for the Fe II ion’s level populations in these ambient conditions with the iteration scheme of Section 2, the 394 levels must be partitioned into superlevels. Because of the importance of the lowest levels, the 16 levels of the terms $`a^6D`$, $`a^4F`$, $`a^4D`$ and $`a^4P`$ are treated individually. Thereafter, taking account of the physical arguments of Hubeny & Lanz (1995), levels are consolidated into a single superlevel if they arise from terms having the same multiplicity and parity and are not widely separated in energy. The resulting reduced ion has 43 levels.

#### 3.2.1 Iterations

Figure 1 shows the performance of the iterative scheme for the specified ambient conditions and with the chosen partitioning. The initial estimate – the zeroth iterate – is the Boltzmann population at $`T_e`$ whose accuracy relative to the exact solution $`n_i^{(x)}`$ is plotted at $`r=0`$. Thereafter, success in converging to $`n_i^{(x)}`$ is measured by $`ϵ^{(r)}`$ and the iterative corrections to $`n_i`$ by $`\delta ^{(r)}`$. From the plotted values of $`ϵ^{(r)}`$, we see that there is a rapid initial improvement followed by slow but steady subsequent improvements, with $`ϵ^{(r)}`$ then decreasing by about 0.1 dex per iteration.

Because of this relatively slow convergence, a large number of iterations are needed if very high accuracy is demanded. Thus, while only 4 iterations give $`ϵ<10^2`$, 12 are required to achieve $`ϵ<10^3`$, 21 for $`ϵ<10^4`$ and 32 for $`ϵ<10^5`$. But given the parametric and structural uncertainties of astrophysical models and the limited accuracy of atomic data, there is no practical need for high accuracy at present. Solutions with $`ϵ10^2`$ are entirely adequate for current diagnostic investigations of observed spectra.

The quantity $`ϵ^{(r)}`$ measures the typical accuracy achieved for the level populations. But the largest fractional error among all 394 levels is also recorded at each iteration. For the test plotted in Figure 1, $`ϵ_{max}<10^1`$ at $`r=3`$, $`<10^2`$ at $`r=11`$ and $`<10^3`$ at $`r=19`$.

In practical applications, the exact solution is not known and so convergence must be decided using a quantity such as $`\delta ^{(r)}`$. Figure 1 shows that, after a few iterations, $`\delta ^{(r)}`$ is significantly smaller than $`ϵ^{(r)}`$, a situation familiar from iteration schemes for radiative transfer problems. In fact, for large $`r`$, we find that $`\delta ^{(r)}/ϵ^{(r)}0.25`$. Accordingly, this test suggests that a convergence criterion of the form $`\delta ^{(r)}<\gamma ϵ`$ with $`\gamma 0.25`$ will give a solution for which the typical fractional error of $`n_i`$ levels is $`ϵ`$. Unfortunately, the speed of convergence and the relationship between $`\delta ^{(r)}`$ and $`ϵ^{(r)}`$ are problem-specific and undoubtedly also partitioning-specific. This particular test illustrates the relatively difficult case near a photosphere where significant radiative excitation of high levels occurs. On the other hand, in very dilute, UV-strong radiation fields, where the normal levels are severely depopulated, the convergence is rapid.

To illustrate this point, the above test is repeated with parameters $`T_b=40000K`$, $`W=7.9\times 10^{13}`$, $`T_e=10000K`$ and $`N_e=4000cm^3`$, representing conditions at a typical point in the Orion Nebula (Lucy 1995). In this case, the scheme achieves $`ϵ^{(r)}<10^3`$ already at $`r=5`$, at which point $`ϵ^{(r)}`$ is decreasing by about 0.64 dex per iteration. Notice that this extremely rapid convergence occurs despite the very poor initial estimate, namely the Boltzmann population at $`T_e`$.

#### 3.2.2 Asymptotic behaviour

The nearly constant logarithmic decrements in $`\delta ^{(r)}`$ and $`ϵ^{(r)}`$ for large $`r`$ seen in Figure 1 indicate that the iterative scheme has a simple asymptotic behaviour, and this becomes strikingly evident when the plot is extended to $`r40`$, for then the two sets of points become linear and parallel. Although no rigorous results have been established, inspection of residuals suggests that

$$\frac{n_i^{(r)}n_i^{(x)}}{n_i^{(r1)}n_i^{(x)}}\alpha asr\mathrm{}$$
(9)
with the constant $`\alpha 0.81`$ for this Fe II test problem.

A useful application of this numerically-suggested asymptotic behaviour is to estimate $`ϵ^{(r)}`$ even when the exact solution $`n_i^{(x)}`$ is not known. It readily follows from equation (9) that as $`r\mathrm{}`$

$$\frac{\delta ^{(r)}}{\delta ^{(r1)}}\alpha andϵ^{(r)}\frac{\alpha }{1\alpha }\delta ^{(r)}$$
(10)
Together these formulae allow $`ϵ^{(r)}`$ to be estimated for $`r2`$. Comparison of these estimates with the exact values of $`ϵ^{(r)}`$ reveals acceptable accuracy for $`r4`$ and the expected high accuracy for large $`r`$. At $`r=4`$, the estimate is $`0.72\times `$ the exact value.

Another possible application is to predict $`n_i^{(x)}`$ from consecutive iterates $`n_i^{(r1)}`$ and $`n_i^{(r)}`$ and thereby achieve accelerated convergence. Thus, after the $`r`$th iteration, we would replace $`n_i^{(r)}`$ by the estimate

$$n_i^{(x)}=n_i^{(r)}+\frac{\alpha }{1\alpha }[n_i^{(r)}n_i^{(r1)}]$$
(11)
and then continue the iterations at step 1). But this modification is not in fact recommended since it jeopardizes the scheme’s robustness and, in any case, by the time the asymptotic behaviour is established, satisfactory accuracy has already been achieved. Nevertheless, if and when highly accurate solutions are required for atoms with large $`N`$, this device for accelerating convergence could be exploited, perhaps stabilized with an undercorrection factor.

#### 3.2.3 Computer time

This iterative scheme has been designed so that, when treating model atoms with $`N10^3`$, the full $`N^3`$ scaling of direct solutions is avoided. Clearly, $`N^3`$ scaling still applies for step 2), but with $`N`$ here reduced from 394 to 43, the time required by the direct solution for the reduced atom is reduced dramatically – by a factor $`1.3\times 10^3`$ – compared to that required for the complete atom. This leaves the initial calculation of the coefficients $`R_{ij}`$ and, during each iteration, the calculation of the superlevel coefficients $`R_{JK}`$ and the relaxation step, all of which scale as $`N^2`$, as the major users of CPU.

The actual times required on a Sun Ultra 1 with clock speed 170 MHz and with 256 Mb of memory are as follows. The direct solution for the 394-level ion takes 21.15s, of which 0.70s is the time required to calculate the $`R_{ij}`$.

In contrast, the time required for the iterative scheme to obtain a solution with 1 percent accuracy is 1.30s. This comprises the 0.70s for the $`R_{ij}`$ and 3 iterations, each of which takes 0.20s. Within an iteration, the breakdown is 0.13s for the $`R_{JK}`$, 0.02s for the direct statistical equilibrium solution for the reduced ion, and 0.05s for the relaxation step.

Although not providing acceptable accuracy, the time required by the superlevel technique is 0.85s. This comprises the times required for calculating the $`R_{ij}`$ and $`R_{JK}`$ coefficients in addition to the trivial time required for the direct solution of the reduced ion.

Timing experiments have also been carried out to test the efficacy of performing $`m`$ relaxation and normalization steps per iteration cycle – see comment j) of Section 2.5. The Fe II test was repeated with $`m=14`$, and the results compared by plotting the accuracy parameter $`ϵ^{(r)}`$ not against $`r`$ but against $`t^{(r)}`$, the elapsed computer time. This comparison shows that the scheme with $`m=2or\mathrm{\hspace{0.33em}3}`$ is initially marginally superior but that for $`r16`$ the standard scheme with $`m=1`$ is best.

### 3.3 Ni I–IV

In view of the scheme’s success in treating a detailed model Fe II ion, we now test its ability to treat ionization as well as excitation by applying it to a multi-ion model of Ni, an atom closely comparable to Fe with regard to the complexity of its level structure.

The model Ni atom has 1266 levels, comprising a 186-level model for Ni I, 717 levels for Ni II, 344 for Ni III, and a severely truncated 19-level model for Ni IV. These energy levels and the f-values for permitted transitions between them were extracted from the Kurucz–Bell (1995) compilation. Data for forbidden transtions are from Garstang (1964) for Ni I, from Nussbaumer & Storey (1982) for Ni II, and from Garstang (1958) and Osterbrock (1992) for Ni III. The van Regemorter (1962) formula is used to compute collision strengths for permitted transitions and 1.0 is the default value for forbidden transitions. Also included are photoionizations to and radiative recombinations from the well-populated lower levels of the ions Ni II–IV. The photoionization cross sections are assumed to follow Seaton’s formula with $`\beta =2`$ and $`s=2`$ (Pauldrach et al. 1994) and to have hydrogenic threshold cross sections. The corresponding total radiative recombination coefficients are derived from Milne’s relation. Collisional ionization is neglected.

The partitioning of the Ni ions into superlevels follows the procedure adopted for Fe II, with several of the lowest levels left single and higher levels grouped into superlevels when they have the same multiplicity and parity and are not widely separated in energy. Specifically, the Ni I levels are consolidated into 24 superlevels with the first 7 being single, the corresponding numbers are 33 and 8 for Ni II, 27 and 9 for Ni III, and 2 and 0 for the truncated Ni IV ion. The net result of this partitioning is to reduce a model atom with 1266 levels to one with 86 superlevels.

As for the Fe II test problem, the ambient radiation field is taken to be $`WB_\nu (T_b)`$ but now with $`T_b=7500K`$ and dilution factor $`W=0.5`$. This value of $`W`$ corresponds to a point at the photosphere ($`r=R`$), and thus the importance of radiative excitation of high levels is enhanced relative to the Fe II test. Electron collisional excitation is included with $`T_e=6500K`$ and $`N_e=10^8cm^3`$.

#### 3.3.1 Iterations

Because of the huge value of $`N`$, no attempt has been made to derive the exact solution $`n_i^{(x)}`$ with a linear equation- solving package. Instead, $`n_i^{(x)}`$ is taken to be the solution obtained after 50 iterations of the iterative scheme, and this is used when calculating $`ϵ^{(r)}`$ from equation (8)

Figure 2, drawn to the same scale as Figure 1, shows the performance of the iterative scheme for the specified ambient conditions and with the chosen partitioning of the Ni atom. The initial estimate is the Saha–Boltzmann population at $`T_e`$. From the plotted values of $`ϵ^{(r)}`$, we see that there is again a rapid initial improvement followed by slow but steady subsequent improvements, with $`ϵ^{(r)}`$ then decreasing by about 0.07 dex per iteration, a somewhat slower rate of convergence for large $`r`$ than that seen in Figure 1.

As in the Fe II test, the largest fractional error among all 1266 levels of the complete atom is also recorded at each iteration. For the test plotted in Figure 2, $`ϵ_{max}<10^1`$ at $`r=6`$, $`<10^2`$ at $`r=18`$ and $`<10^3`$ at $`r=32`$.

The nearly constant logarithmic decrements in $`\delta ^{(r)}`$ and $`ϵ^{(r)}`$ for large $`r`$ evident in Figure 2 indicate that the iterative scheme’s simple asymptotic behaviour when applied to a single ion also applies for a multi-ion atomic model. Inspection of residuals confirms that equation (10) again holds asymptotically, but with slightly different coefficients $`\alpha `$ for the different ions.

#### 3.3.2 Error distributions

In Figure 3, the distribution of the absolute values of the logarithmic errors in the level populations obtained with the superlevel technique – i.e., after step 3) of the first iteration – is compared with the corresponding error distributions after the first three iterations of the iterative scheme.

Now, given the accuracy expected of modern diagnostic investigations of astronomical spectra, errors in level populations of 0.1 dex are surely the maximum acceptable. But for the superlevel solution plotted in Figure 3, no less than 37.4 percent of the levels have errors exceeding 0.1 dex, and there is an extended high-error tail to beyond 0.3 dex. Thus, despite a physically-sensible partitioning and a relatively large $`N`$ for the reduced atom, the superlevel solution, if used in a diagnostic analysis, would significantly limit the accuracy of that analysis. This conclusion is in accord with the finding of Hillier and Miller (1998) that the superlevel technique is not sufficiently accurate for diagnostically-important emission lines.

The comparison in Figure 3 of the errors of the superlevel solution with those after 1 iteration of the iterative scheme shows the dramatic effect of the relaxation step in immediately decreasing the number of levels with substantial errors. Levels with errors exceeding 0.1 dex drop to 12.2 percent in consequence of the relaxation step, and then drop to 2.9 and 0.2 percent with two further iterations. Clearly, in this case, even with as few as 3 iterations, the iterative scheme provides a solution that could be used in a diagnostic analysis with the confident expectation that errors in solving for statistical equilibrium are unlikely to be compromising the accuracy of the results.

#### 3.3.3 Computer time

By scaling from the Fe II case, we can estimate that the direct solution for the 1266-level would take 678s, assuming that the build up of round off errors did not in fact prevent a solution. This is to be compared with an iteration time of 4.66s. The breakdown of this is 2.08s to calculate the $`R_{JK}`$, 0.18s for the direct solution of the 86-level reduced atom, and 1.40s for the relaxation step. With just a few iterations needed to achieve the desired accuracy, there is evidently a vast saving in computer time. Moreover, the direct solution is in any case not a realistic option.

### 3.4 Line formation

In the above test problems, the ambient radiation field is a featureless continuum. But typically $`J_\nu `$ will be strongly affected by line formation, and so it is important to investigate how this impacts on the performance of the scheme.

The simplest test incorporating line formation assumes that the point considered is in an expanding flow and treats line formation in the Sobolev approximation. The b-b rate coefficients are then modified as follows (Klein & Castor 1978): line trapping is included by replacing the Einstein coefficient $`A_{ul}`$ by $`A_{ul}\beta _{ul}`$, where $`\beta _{ul}`$ is the Sobolev escape probability; and the profile-averaged mean intensity, needed for the rates of radiative excitation and de-excitation, is $`\beta _{ul}J_{lu}`$, where $`J_{lu}`$ is the mean intensity of the unattenuated continuum

With these modifications, the rate coefficients $`R_{ij}`$ depend on the Sobolev optical depth and therefore on $`n_i`$. In consequence, the statistical equilibrium problem is no longer linear. Fortunately, simple repeated back substitutions (Lucy 1999b) allow $`n_i^{(x)}`$ to be determined with high accuracy in $`5`$ iterations.

A second consequence of the dependence of the $`R_{ij}`$ on $`n_i`$ is that, in applying the iterative scheme, the $`R_{ij}`$ must be recalculated at the beginning of each iteration cycle before calculating the coefficients $`R_{JK}`$ for the superlevels.

#### 3.4.1 Test problem

The effect of line formation on the iterative scheme is now investigated by repeating the test of Section 3.2, but we now suppose that the point considered is in a supernova’s homologously expanding envelope, that the elapsed time since explosion is 13 days, and that the number density of Fe II ions is $`5.1\times 10^5cm^3`$. These two additional parameters determine the Sobolev optical depths of the b-b transitions and are such that the converged solution has numerous optically-thick UV lines, with optical depths up to $`2\times 10^4`$.

As with the previous tests, convergence is monitored by computing $`ϵ^{(r)}`$, and this reveals faster convergence than when line formation is ignored, with $`ϵ^{(r)}`$ now dropping by 0.15 dex per iteration at large $`r`$ as against 0.1 dex previously. Because of this faster convergence, only 2 iterations are required to achieve $`ϵ<10^2`$, 4 for $`ϵ<10^3`$, 9 for $`ϵ<10^4`$ and 16 for $`ϵ<10^5`$.

As in the more extreme Orion Nebula case of Section 3.2.1, the faster convergence in this test can be attributed to the reduced populations of high normal levels. When high levels are depleted, either by dilution or by line formation, the validity of the assumption motivating the relaxation step – that high levels are negligibly coupled to other high levels – is strengthened – see Section 2.1.

## 4 Conclusion

The aim of this paper has been to describe and test an iterative technique for obtaining approximate solutions to the equations of statistical equilibrium with the accuracy required by modern diagnostic analyses of astronomical spectra. In particular, the technique is designed to treat model atoms whose numbers of levels $`N`$ are so large that exact solutions with standard software packages are precluded, either because of excessive computer time or because of an anticipated loss of accuracy due to the accumulation of round-off errors. By combining Anderson’s (1989) superlevel technique with a simple relaxation method, solutions of the desired accuracy are derived in $``$ 3–5 iterations with a huge saving of computer time compared to direct solutions, and for values of $`N`$ for which direct solutions are not a feasible option.

As mentioned in Section 1, this technique was developed for Monte Carlo codes. But the most promising immediate application would seem to be to the already-existing Accelerated Lambda Iteration (ALI) codes that use superlevel partitioning to avoid treating atomic models with excessively large $`N`$. In view of Figure 3 and the earlier work of Hillier & Miller (1998), the accuracy of these codes almost certainly falls short of minimum requirements for the quantitative analyses of observed spectra. Fortunately, the iterative technique described herein would seem to be readily incorporated into ALI codes and would bring them to the desired accuracy with little impact on their complexity, robustness or use of computer time.


format: text

text:
# Short-time dynamics of random-bond Potts ferromagnet with continuous self-dual quenched disorders1footnote 11footnote 1Supported by the NNSF of China, 19975041 and 10074055, and by the Deutsche Forschungsgemeinschaft, DFG Schu 95/9-3.

## I Introduction

Ferromagnetic systems with quenched randomness are studied intensively for their critical properties affected by an addition of disorders Cardy96 ; Jaco98 ; Folk00 . The Harris criterion was claimed 27 years ago that if the specific heat critical exponent $`\alpha `$ of the pure system is positive, then quenched disorder is a relevant perturbation at the second-order critical point and it causes changes in critical exponents Harris . Then, following the work of Imry and Wortis Imry79  who argued that quenched disorder could smooth of a first-order phase transition and thus produce a new second-order phase transition, the introduction of randomness to pure systems originally undergoing a first-order transition has been comprehensively considered Hui89 ; Aizen89 ; Cardy99 . The theory was firstly numerically checked with the Monte Carlo (MC) method by Chen, Ferrenberg and Landau (CFL) Chen92  who studied the 8-state Potts model with a self-dual random-bond disorder (see Eqs.(2) and (3)). Recently most of the work for such disordered systems is carried out by intensive MC simulations Wise95 ; Kardar95 ; Kim96 ; Yasar98 ; Chate98 ; Olson99  for examining how a phase transition is modified by the quenched disorder coupling to the local energy density for the $`q`$-state random-bond Potts ferromagnet (RBPF).

The $`q`$-state RBPF, an interesting framework to study the influence of impurity on pure systems, is described by the Hamiltonian,

$`\beta H={\displaystyle \underset{<i,j>}{}}K_{ij}\delta _{\sigma _i\sigma _j},K_{ij}>0,`$ (1)
where $`\beta =1/k_BT`$ is the inverse temperature, the spin $`\sigma `$ can take the values 1,$`\mathrm{},q`$, $`\delta `$ stands for the Kronecker delta function, and the nearest-neighbor interactions $`<i,j>`$ are considered. Usually the dimensionless couplings $`K_{ij}`$ are selected from two positive (ferromagnetic) values $`K_1`$ and $`K_2=rK_1`$, with a strong to weak coupling ratio $`r=K_2/K_1`$, called a disorder amplitude, according to a bimodal distribution,

$`P(K)=p\delta (KK_1)+(1p)\delta (KK_2).`$ (2)
With $`p=0.5`$, the system is self-dual and the exact critical point can be determined by Kin81 ,

$`(e^{K_c}1)(e^{K_c^{}}1)=q,`$ (3)
where $`K_c`$ and $`K_c^{}`$ are the corresponding critical values of $`K_1`$ and $`K_2`$, respectively. The value $`r=1`$ corresponds to the pure case where the critical point is located at $`K_c=\text{log}(1+\sqrt{q})`$ and the phase transitions are first-order for $`q>4`$. As $`\alpha >0`$ for all $`q>2`$ the randomness acts as a relevant perturbation and, for $`q>4`$, it even changes the nature of the transition from first to second order.

Up to present, the MC studies are usually carried out in equilibrium and concentrated on disorder amplitudes in the range of $`r=2`$-20 adapted to numerical analysis. They give good estimates of the disordered fixed point exponents Cardy99 ; Picco98 . In Table 1 we present the magnetic scaling index $`\eta =2\beta /\nu `$ of the 2D 8-state RBPF obtained by different groups Cardy97 ; Chate98a ; Chate99 ; Picco98 ; Ying00 . In a recent work, Olson and Young (OY) Olson99  used their special self-dual disorder distribution and performed the MC study of multiscaling properties of the correlation functions for several values of $`q`$. Their results for the magnetic exponent are very interesting to examine the universality of the RBPF in the crossover regime from the pure fixed point to a percolation-like limit Picco98 . Cardy and Jacobsen Cardy97  studied the RBPF based on the connectivity transfer matrix (TM) formalism, and their estimates of the critical exponents lead to a continuous variation of $`\eta `$ with $`q`$, which is in sharp disagreement with CFL’s results Chen92 .

In this paper, we present a MC study to verify the dynamic scaling of the RBPF and estimate the critical exponents in the short-time dynamics (STD) Janss89 ; Zheng98  for the OY’s continuous distribution of the disorder interactions. It is a self-dual scheme to introduce quenched randomness and defined by the following probability

$`P_X(x)={\displaystyle \frac{2\sqrt{q}}{\pi (1x)^2+qx^2}},`$ (4)
for the Boltzmann factor $`xe^K`$. It is obvious that the model has strong disorder as the distribution function is nonzero everywhere at $`x(0,1)`$ ($`K\mathrm{}`$, $`K0`$) and has a finite weight at its two limit points $`x=0`$ and $`x=1`$. Another advantage of this distribution is that it is easy to generate random numbers with the probability $`P_X(x)`$ through

$`x={\displaystyle \frac{1}{1+\sqrt{q}\mathrm{tan}(\pi R_x/2)}},`$ (5)
where $`R_x`$ is a random number with uniform distribution between 0 and 1.

In particular, we will investigate the critical behavior affected by introduction of such continuous quenched randomness to clarify the universality class of the RBPF in the crossover regime. In simulations we mostly choose the model with $`q=8`$ which is known to have a strong first-order phase transition without disorders. We hope to confirm the new second-order phase transition induced by the random disorder in order to show the influence of quenched impurities on the first-order systems and to check whether an Ising-like universality would be satisfied from the STD approach.

## II Short-time Dynamics

For a long time, it was believed that universality and scaling relations can be found only in equilibrium or in the long-time regime. In Refs. Janss89 ; Zheng98 , however it was discovered that for a $`O(N)`$ vector magnetic system in states with a very high temperature $`TT_c`$, when it is suddenly quenched to the critical temperature $`T_c`$, and evolves according to a dynamics of model A, a universal dynamic scaling behavior emerges already within the short-time regime,

$$M^{(k)}(t,\tau ,L,m_0)=b^{k\beta /\nu }M^{(k)}(b^zt,b^{1/\nu }\tau ,b^1L,b^{x_0}m_0),$$
(6)
where $`M^{(k)}`$ is the $`k`$th moment of magnetization, $`t,\tau =(TT_c)/T_c,L`$ and $`b`$ are time, reduced temperature, lattice size and scaling factor, respectively. $`\beta `$ and $`\nu `$ are the well known static critical exponents. The $`x_0`$, a new independent exponent, is the scaling dimension of initial magnetization $`m_0`$. This dynamic scaling form is generalized from finite size scaling in the equilibrium stages. The dynamic MC simulations have been successfully performed in such non-equilibrium processes to estimate the critical point and the critical exponents ($`K_c`$,$`z`$,$`\beta `$,$`\nu `$,$`\theta `$) Li95 ; Zheng99 ; Ying98 ; Ying01 .

The MC simulation begins with a study of the time evolution in the initial stage of the dynamic relaxation starting at very high temperature and small magnetization ($`m_00`$). For a sufficiently large lattice ($`L\mathrm{}`$), from Eq.(6) by setting $`\tau =0`$ and $`b=t^{1/z}`$, it is easy to derive that

$`M^{(k)}(t,m_0)=t^{k\beta /\nu z}M^{(k)}(1,t^{x_0/z}m_0).`$ (7)
When $`k=1`$ we get the most important scaling relation on which our measurements of the critical exponent $`\theta `$ are based,

$`M(t)m_0t^\theta ,\theta =(x_0\beta /\nu )/z.`$ (8)
In most cases the magnetization undergoes an initial increase at the critical point $`K_c`$ after a microscopic time $`t_{mic}`$.

In our MC sweeps, the time evolution of $`M(t)`$ is calculated through the definition

$`M(t)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\left[{\displaystyle \frac{qM_ON}{q1}}\right].`$ (9)
Here $`M_O=\text{max}(M_1,M_2,\mathrm{},M_q)`$ with $`M_i`$ being the number of spins in the $`i`$th state among $`q`$ states. $`N=L^2`$ is the number of spins on the suqare lattice. We use $`L`$ up to 128. $`<\mathrm{}>`$ denotes thermal averages over independent initial states and random number sequences, and $`[\mathrm{}]`$ the disorder averages over quenched randomness distributions. The unit of time $`t`$ is defined as a MC sweep over all spins on the lattice.

The susceptibility plays an important role in the equilibrium. Its finite size behavior has been often used to determine the critical temperature and the critical exponents $`\gamma /\nu `$ and $`\beta /\nu `$ Chen92 . For the STD approach, the time-dependent susceptibility (the second moment of the magnetization) is also interesting and important. For the random-bond Potts model, it is defined as

$`M^{(2)}(t)`$ $`=`$ $`{\displaystyle \frac{1}{N}}\left[\left(<M^2(t)><M(t)>^2\right)\right].`$ (10)
One studies its scaling behavior from $`m_0=0`$ in short-time relaxation. Because the spatial correlation length in the beginning of the evolution is small compared with the lattice size $`L^d`$ the second moment behaves as $`M^{(2)}(t,L)L^d`$. Then the finite size scaling Eq.(6) induces a power-law behavior at the critical temperature,

$`M^{(2)}(t)t^y,y=(d\eta )/z.`$ (11)
In the above considerations the dynamic relaxation was assumed to start from disordered states with $`m_0`$ small enough. Another interesting and important process is the dynamic relaxation from a completely ordered state. The initial magnetization is taken exactly at its other fixed point $`m_0=1`$, where a critical scaling form

$`M^{(k)}(t,\tau ,L)=b^{k\beta /\nu }M^{(k)}(b^zt,b^{1/\nu }\tau ,b^1L)`$ (12)
is expected Ito93 ; Zheng98 . This scaling form looks to be the same as the dynamic scaling form in the long-time regime, however, it is now assumed already valid in the macroscopic short-time regime. For the magnetization itself, $`b=t^{1/z}`$ yields, for a sufficiently large lattice and at the critical point ($`\tau =0`$), a power-law decay behavior of

$`M(t,\tau )=t^{c_1},c_1=\beta /\nu z.`$ (13)
The formula can be used to calculate the critical exponents $`\eta =2\beta /\nu `$ or $`z`$.

We further observe an evolution of the Binder cumulant $`U(t,L)=M^{(2)}(t,L)/M^2(t,L)1`$. A similar power-law behavior at criticality induced from the Eq.(12) shows that

$`U(t,L)t^{c_u},c_u=d/z`$ (14)
on large enough lattices. Here, unlike the relaxation from the disordered state, the fluctuations caused by the initial configurations are much smaller. In practical simulations, these measurements of the critical exponents and critical temperature are better in quality than those from the dynamical relaxation starting from disordered states.

## III MC Results

In simulations, the heat-bath algorithm is used for MC updating, and up to 500 samples are chosen as disorder averages and 200-500 initial configurations and/or random number sequences are taken as thermal averages for each disorder realization (so total MC averages are over 100,000). Statistical errors are simply estimated by performing three groups of averages with different random number sequences as well as independent initial states. We first focus our attention on evolutions of the magnetization both from the initial states with small magnetization ($`m_00`$) and complete ordered states ($`m_0=1`$).

At the critical point ($`\tau =0`$), achieved by generating the Boltzmann factors according to the distribution of Eq.(5), the curves $`M(t)`$ should show power-law behavior as described by Eq.(8) or Eq.(13) dependent on the initial states. They, however, should deviate from the power-law scaling behavior for a small but nonzero $`\tau `$ characterized by $`M(t,\tau )=m_0t^\theta M(t^{1/\nu z}\tau )`$ or $`M(t,\tau )=t^{c_1}F(t^{1/\nu z}\tau )`$ with respect to the $`m_00`$ or $`m_0=1`$ initial states respectively. Now a question arises how to introduce a distribution which deviates from the critical point with distance $`\tau `$? We can first generate the Boltzmann factors, $`xe^K`$, from the self-dual distribution $`P_X(x)`$ in Eq.(4), then modify them by the replacement $`xx^{1/(1+\tau )}`$ away from criticality Olson99 . We observe that the best power-law behavior appeared at the same critical point for both the random and ordered initial states. This fact indicates a second-order phase transition induced at $`K_c`$, as argued by Schülke and Zheng Schu00  since $`K_{order}<K_c<K_{random}`$ if a first-order phase transition happened. We exhibit in Fig.1 and Fig.2 the double-logarithmic plots of magnetization versus the MC time on a $`L=64`$ lattice to show this behavior obviously.

Next, we carry out simulations at the exact critical point $`K_c`$ to investigate the magnetization and Binder cumulant, starting from $`m_0=1`$ ordered state. Such state is a ground state where all the spins orient in the same direction (e.g., $`\sigma _i`$ = 1, $`i=1,N`$). In Figs.3 and 4, $`M(t)`$ and $`U(t)`$ are presented and values of the scaling indices $`c_1=\beta /\nu z`$ and $`c_u=d/z`$ can be estimated from slopes of the curves in $`t=`$(100,1000). The results are listed in Table 2. Here we find that, unlike for relaxations from the disordered state, the fluctuations caused by ordered initial configurations are much smaller. As a result, these measurements of the critical exponents based on Eqs.(13) and (14) are better in quality than those from disordered states on Eqs.(8) and (11). So we take only 200 MC configurations for the thermal averages in the MC simulations. Now, from the results of $`c_1`$ and $`c_u`$, the critical exponents $`\eta `$ and $`z`$ can be calculated based on Eqs.(13) and (14), and their values are summarized in Table 2.

To estimate the critical exponent $`\theta `$, the random initial states are prepared with small magnetizations $`m_0=0.04`$ – 0.01. The curves of $`M(t)`$ for different $`m_0`$ on $`128^2`$ lattice with a double-logarithmic scale are displayed in Fig.5 which exhibits the finite $`m_0`$ effect on $`M(t)`$ and a very nice power-law behavior at $`m_00`$. Thus the exponent $`\theta `$ can be estimated from slopes of the curves in $`t`$=(50,1000) by using Eq.(8) after an extrapolation to $`m_00`$. Its values are shown also in Table 2.

As mentioned above, for a determination of the dynamic exponent $`z`$ and the magnetic exponent $`\eta `$ a dynamic process starting from the ordered state ($`m_0=1`$) is more favorable since the thermal fluctuation is much smaller than for random states. Therefore we have mainly studied the short-time dynamic process from such an ordered state to calculate $`z`$ and $`\eta `$. On the other hand it is also desirable to measure these exponents from the $`m_0=0`$ initial states Li95 . In our work we now take the evolution for the second moment of magnetization $`M^2(t)`$ from the $`m_0=0`$ to verify the results of $`\eta `$ and $`z`$. By input of values of $`\eta `$ and $`z`$ already obtained from the indices $`c_1`$ and $`c_u`$, we ckeck if they also satify the relation $`y=(d\eta )/z`$ in Eq.(11). As expected we find they are consistent with less than 2% error for $`L96`$.

## IV Conclusions

We have made a first attempt to apply the STD approach to the 2D RBPF with a continuous distribution of quenched disorders. Our simulations verified that a second order phase transition has been induced. Then the universal dynamic scaling behavior in the short-time dynamics is used to estimate the exponents $`\theta `$ and $`z`$, as well the magnetic critical exponent $`\eta `$. It is found that they really violate the Ising-like universality class where $`\theta =0.191(1)`$ and $`\eta =0.240(15)`$ Li95 . Our result of the dynamical exponent $`z=3.41(6)`$ is much bigger than that for systems without disorder, and it causes the microscopic time $`t_{mic}`$ to be extended up to one hundred for the Binder cumulant from the ordered initial state. The value for the exponent $`\eta `$ nearly completely overlaps with those listed in Table 1, except that it is somewhat smaller than that given in Ref.Olson99 . More accurate calculation remains to be done to solve this controversy.

## V acknowledgements

We acknowledge stimulating discussions with Prof. L. Schülke. H.P.Y. would like to thank the Heinrich-Hertz-Stiftung for a fellowship and be grateful to the hospitality of the University of Siegen where the MC simulations were performed.