15.3: Rate Equations for Statistical Equilibrium

The condition that the sum of all transitions into and out of any specific level must be zero implies that there is no net change of any level populations. This means that we can write an expression that describes the flow into and out of each level, incorporating the detailed physics that governs the flow from one level to another. These expressions are known as the rate equations for statistical equilibrium. The unknowns are the level populations for each energy level which will appear in every expression for which a transition between the respective states is allowed. Thus we have a system of n simultaneous equations for the level populations of n states. Unfortunately, as we saw in estimating the rates of collisional ionization and photoionization, it is necessary to know the radiation field to determine the coefficients in the rate equations. Thus any solution will require self-consistency between the radiative transfer solution and the statistical equilibrium solution. Fortunately, a method for the solution of the radiative transfer and statistical equilibrium equations can be integrated easily in the iterative algorithm used to model the atmosphere (see Chapter 12). All that is required is to determine the source function in the line appropriate for the non-LTE state.

Since an atom has an infinite number of allowed states as well as an infinite number of continuum states that must be considered, some practical limit will have to be found. For the purpose of showing how the rate equations can be developed, we consider two simple cases.

a. Two-Level Atom

It is possible to describe the transitions between two bound states we did for photo- and collisional ionization. Indeed, for the radiative processes, basically we have already done so in (Section 11.3) through the use of the Einstein coefficients. However, since we are dealing with only two levels, we must be careful to describe exactly what happens to a photon that is absorbed by the transition from level 1 to level 2. Since the level is not arbitrarily sharp, there may be some redistribution of energy within the level. Now since the effects of non-LTE will affect the level populations at various depths within the atmosphere, we expect these effects will affect the line profile as well as the line strength. Thus, we must be clear as to what other effects might change the line profile. For that reason, we assume complete redistribution of the line radiation. This is not an essential assumption, but rather a convenient one.

If we define the probability of the absorption of a photon at frequency \(\nu’\) by \[\phi\left(\nu^{\prime}\right) d \nu^{\prime}=\left[\int_{-\infty}^{+\infty} R\left(\nu^{\prime}, \nu\right) d \nu\right] d \nu^{\prime}\label{15.2.1}\]

and the probability of reemission of a photon at frequency \(\nu\) as \[\psi(\nu) d \nu=\left[\int_{-\infty}^{+\infty} R\left(\nu^{\prime}, \nu\right) d \nu^{\prime}\right] d \nu\label{15.2.2}\]

then the concept of the redistribution function describes to what extent these photons are correlated in frequency. In Chapter 9, we introduced a fairly general notion of complete redistribution by stating that \(\nu'\) and \(\nu\) would not be correlated. Thus, \[\phi\left(\nu^{\prime}\right)=\psi(\nu)\label{15.2.3}\]

Under the assumption of complete redistribution, we need only count radiative transitions by assuming that specific emissions are unrelated to particular absorptions. However, since the upward radiative transitions in the atom will depend on the availability of photons, we will have to develop an equation of radiative transfer for the two-level atom.

Equation of Radiative Transfer for the Two-Level Atom In Chapter 11 [equations \ref{11.3.6} and \ref{11.3.7}] we described the emission and absorption coefficients, \(\mathrm{j_\nu}\), and \(\kappa_\nu\), respectively, in terms of the Einstein coefficients. Using these expressions, or alternatively just balancing the radiative absorptions and emissions, we can write an equation of radiative transfer as \[\mu \frac{d I_\nu}{d x}=-N_1 B_{12}\left(\frac{h \nu}{4 \pi}\right) \phi_\nu I_\nu+N_2\left(A_{21}+B_{21} I_\nu\right)\left(\frac{h \nu}{4 \pi}\right) \phi_\nu\label{15.2.4}\]

This process of balancing the transitions into and out of levels is common to any order of approximation in dealing with statistical equilibrium. As long as all the processes are taken into account, we will obtain an expression like equation \ref{15.2.4} for the transfer equation for multilevel atoms [see equation \ref{15.2.25}]. Equation \ref{15.2.4} can take on a somewhat more familiar form if we define \[d \tau_\nu=\frac{-\left(N_1 B_{12}-N_2 B_{21}\right) h \nu d x}{4 \pi}\label{15.2.5}\]

Then the equation of transfer becomes \[\mu \frac{d I_\nu}{d \tau_\nu}=\phi_\nu\left(I_\nu-S_{\ell}\right)\label{15.2.6}\]

where \[S_{\ell} \equiv \frac{N_2 A_{21}}{N_1 B_{12}-N_2 B_{21}}\label{15.2.7}\]

Making use of the relationships between the Einstein coefficients determined in Chapter 11 [equation \ref{11.3.5}], we can further write \[S_{\ell}=\frac{2 h \nu^3}{c^2}\left(\frac{N_1 g_2}{N_2 g_1}-1\right)^{-1}\label{15.2.8}\]

Under conditions of LTE \[\frac{N_1 g_2}{N_2 g_1}=e^{\epsilon_i /(k T)}=e^{h \nu_{12} /(k T)}\label{15.2.9}\]

so that we recover the expected result for the source function, namely \[\label{15.2.10}\]

Two-Level-Atom Statistical Equilibrium Equations The solution to equation \ref{15.2.6} will provide us with a value of the radiation field required to determine the number of radiative transitions. Thus the total number of upward transitions in the two-level atom is \[N_{1 \rightarrow 2}=N_1 B_{12} \int \phi_\nu J_\nu d \nu+N_1 N_e \Omega_{12}\label{15.2.11}\]

Similarly, the number of downward transitions is \[N_{2 \rightarrow 1}=N_2 A_{21}+N_2 B_{21} \int \phi_\nu J_\nu d \nu+N_2 N_e \Omega_{21}\label{15.2.12}\]

The requirement that the level populations be stationary means that \[N_{1 \rightarrow 2}=N_{2 \rightarrow 1}\label{15.2.13}\]

so that the ratio of level populations is \[\frac{N_1}{N_2}=\frac{A_{21}+B_{21} \int \phi_\nu J_\nu d \nu+N_e \Omega_{21}}{B_{12} \int \phi_\nu J_\nu d \nu+N_e \Omega_{12}}\label{15.2.14}\]

Now consider a situation where there is no radiation field and the collisions are driven by particles characterized by a maxwellian energy distribution. Under these conditions, the principle of detailed balancing requires that \[N_1 \boldsymbol{\Omega}_{12}=N_2 \boldsymbol{\Omega}_{21}\label{15.2.15}\]

or \[\Omega_{12}=\frac{g_2}{g_1} \Omega_{21} e^{-h \nu /(k T)}\label{15.2.16}\]

This argument is similar to that used to obtain the relationships between the Einstein coefficients and since the collision coefficients depend basically on atomic constants, equation \ref{15.2.16} must hold under fairly arbitrary conditions. Specifically, the result will be unaffected by the presence of a radiation field. Thus we may use it and the relations between the Einstein coefficients [equations \ref{11.3.5}] to write the line source function as \[S_{\ell}=\frac{\int \phi_\nu J_\nu d \nu+\frac{N_e \Omega_{21}}{A_{21}} \frac{2 h \nu^3}{c^2} e^{-h \nu /(k T)}}{1+\frac{N_e \Omega_{21}}{A_{21}}\left(1-e^{-h \nu /(k T)}\right)}\label{15.2.17}\]

If we let \[\epsilon=\frac{N_e \Omega_{21}\left(1-e^{-h \nu /(k T)}\right) / A_{21}}{1+N_e \Omega_{21}\left(1-e^{-h \nu /(k T)}\right) / A_{21}}\label{15.2.18}\]

then the source function takes on the more familiar form \[S_\ell=\epsilon B_\nu+(1-\epsilon) \int \phi_{\nu^{\prime}} J_{\nu^{\prime}}, d \nu^{\prime}\label{15.2.19}\]

The quantity ε is, in some sense, a measure of the departure from LTE and is sometimes called the departure coefficient. A similar method for describing the departures from LTE suffered by an atom is to define \[b_j=\frac{N_j}{\tilde{N}_j}\label{15.2.20}\]

where \(\mathrm{N_j}\) is the level population expected in LTE so that \(b_\mathrm{j}\) is just the ratio of the actual population to that given by the Saha-Boltzmann formula. From that definition, equation \ref{15.2.14}, and the relations among the Einstein coefficients we get \[1-\frac{b_1}{b_2}=\frac{\int \phi_\nu J_\nu d \nu-B_\nu}{\left(1-e^{-h \nu /(k T)}\right)^{-1} \int \phi_\nu J_\nu d \nu+\left(N_e \Omega_{21} / A_{21}\right) B_\nu(T)}\label{15.2.21}\]

b. Two-Level Atom plus Continuum

The addition of a continuum increases the algebraic difficulties of the above analysis. However, the concepts of generating the statistical equilibrium equations are virtually the same. Now three levels must be considered. We must keep track of transitions to the continuum as well as the two discrete energy levels. Again, we assume complete redistribution within the line so that the line source function is given by equation \ref{15.2.8}, and the problem is to find the ratio of the populations of the two levels.

We begin by writing the rate equations for each level which balance all transitions into the level with those to the other level and the continuum. For level 1, \[\begin{aligned}
& N_1\left(B_{12} \int \phi_\nu J_\nu d \nu+N_e \Omega_{12}+R_{1 k}+N_e \Omega_{1 k}\right) \\
& \quad=N_2\left(A_{21}+B_{21} \int \phi_\nu J_\nu d \nu+N_e \Omega_{21}\right)+N_1^*\left(R_{k 1}+N \Omega_{1 k}\right)
\end{aligned}\label{15.2.22}\]

The parameter \(\mathrm{R_ik}\) is the photoionization rate defined in equation \ref{15.1.1}, while \(\mathrm{R_ki}\) is the analogous rate of photorecombination. When the parameter Ω contains the subscript k, it refers to collisional transitions to or from the continuum. The term on the left-hand side describes all the types of transitions from level 1 which are photo-and collisional excitations followed by the two terms representing photo- and collisional, ionizations respectively. The two large terms on the right-hand side contain all the transitions into level 1. The first involves spontaneous and stimulated radiative emissions followed by collisionly stimulated emissions. The second term describes the recombinations from the continuum. The parameter \(\mathrm{N_i}^*\) will in general represent those ions that have been ionized from the ith state.

We may write a similar equation \[\begin{aligned}
& N_2\left(A_{21}+B_{21} \int \phi_\nu J_\nu d \nu+N_e \Omega_{21}+R_{2 k}+N_e \Omega_{2 k}\right) \\
& \quad=N_1\left(B_{12} \int \phi_\nu J_\nu d \nu+N_e \Omega_{12}\right)+N_2^*\left(R_{k 2}+N_e \Omega_{2 k}\right)
\end{aligned}\label{15.2.23}\]

for level 2 by following the same prescription for the meaning of the various terms. Again letting the terms on the left-hand side represent transitions out of the two levels while terms on the right-hand side denote inbound transitions, we find the rate equation for the continuum is \[\begin{aligned}
& N_1^*\left(R_{k 1}+N_e \Omega_{1 k}\right)+N_2^*\left(R_{k 2}+N_e \Omega_{2 k}\right) \\
& \quad=N_1\left(R_{1 k}+N_e \Omega_{1 k}\right)+N_2\left(R_{2 k}+N_e \Omega_{2 k}\right)
\end{aligned}\label{15.2.24}\]

However, this equation is not linearly independent from the other two and can be generated simply by adding equations \ref{15.2.22} and \ref{15.2.23}. This is an expression of continuity and will always be the case regardless of how many levels are considered. There will always be one less independent rate equation than there are levels. An electron that leaves one state must enter another, so its departure is not independent from its arrival. If all allowed levels are counted, as they must be if the equations are to be complete, this interdependence of arrivals and departures of specific transitions will make the rate equation for one level redundant. Noting that the same kind of symmetry described by equation \ref{15.2.15} also holds for the collisional ionization and recombination coefficients, we may solve equations \ref{15.2.22} and \ref{15.2.23} for the population ratio required for the source function given by equation \ref{15.2.8}. The algebra is considerably more involved than for the two levels alone and yields³ a source function of the form \[\begin{aligned}
S_{\ell} & =\frac{\int \phi_\nu J_\nu d \nu+\tilde{\epsilon} B_\nu(T)+\eta B^*}{1+\tilde{\epsilon}+\eta} \\
\tilde{\epsilon} & =\frac{N_e \Omega_{21}}{A_{21}}\left(1-e^{-h \nu /(k T)}\right) \\
\eta & =\frac{1}{A_{21}} \frac{\left(R_{2 k}+N_e \Omega_{2 k}\right) N_1^*\left(R_{k 1}+N_e \Omega_{1 k}\right)-\left(g_1 / g_2\right)\left(R_{1 k}+N_e \Omega_{1 k}\right) N_2^*\left(R_{2 k}+N_e \Omega_{2 k}\right)}{N_1^*\left(R_{k 1}+N_e \Omega_{1 k}\right)+N_2^*\left(R_{k 2}+N_e \Omega_{2 k}\right)} \\
B^* & =\frac{2 h \nu^3}{c^2}\left[\frac{N_1^* g_2\left(R_{2 k}+N_e \Omega_{2 k}\right)\left(R_{k 1}+N_e \Omega_{1 k}\right)}{N_2^* g_1\left(R_{1 k}+N_e \Omega_{1 k}\right)\left(R_{k 2}+N_e \Omega_{2 k}\right)}-1\right]^{-1}
\end{aligned}\label{15.2.25}\]

If the terms involving ε dominate the source function, the lin is said to be collisionly dominated, while if the terms involving η are the largest, the line is said to be dominated by photoionization. If \(\bar{\varepsilon} B_\nu(T)>\eta B^*\) but \(\eta>\bar{\varepsilon}\) (or vice versa), the line is said to be mixed. Some examples of lines in the solar spectrum that fall into these categories are given in Table 15.1.

Table 15.1 Types of Solar Spectral Lines

Collisionly Dominated	Dominated by Photoionization
Resonance lines of singly ionized metals Resonance lines of hydrogen Resonance lines of nonmetals	Resonance lines of neutral metals Balmer lines of hydrogen

c. Multilevel Atom

A great deal of effort has gone into approximating the actual case of many levels of excitation by setting up and solving the rate equations for three and four levels or approximating any particular transition of interest by an "equivalent two level atom" (see D.Mihalas³, pp. 391-394). However, the advent of modern, swift computers has made most of these approximations obsolete. Instead, one considers an n-level atom (with continuum) and solves the rate equations directly. We have already indicated that this procedure can be integrated into the standard algorithm for generating a model atmosphere quite easily. Consider the generalization of equations \ref{15.2.22} through \ref{15.2.24}. Simply writing equations for each level, by balancing the transitions into the level with those out of the level, will yield a set of equations which are linear in the level populations. However, as we have already indicated, these equations are redundant by one. So far we have only needed population ratios for the source function, but if we are to find the population levels themselves, we will need an additional constraint. The most obvious constraint is that the total number of atoms and ions must add up to the abundance specified for the atmosphere. Mihalas⁴ suggests using charge conservation, which is a logically equivalent constraint. Whatever additional constraint is chosen, it should be linear in the level populations so that the linear nature of the equations is not lost.

It is clear that the equations are irrevocably coupled to the radiation field through the photoexcitation and ionization terms. It is this coupling that led to the rather messy expressions for the source functions of the two-level atom. However, if one takes the radiation field and electron density as known, then the rate equations have the form \[\mathrm{A} \vec{N}=\vec{B}\label{15.2.26}\]

where A is a matrix whose elements are the coefficients multiplying the population levels and \(\vec{N}\) is a vector whose elements are the populations of the energy levels for all species considered in the calculation. The only nonzero element of the constant vector \(\vec{B}\) arises from the additional continuity constraint that replaced the redundant level equation. These equations are fairly sparse and can be solved quickly and accurately by well-known techniques.

Since the standard procedure for the construction of a model atmosphere is an iterative one wherein an initial guess for the temperature distribution gives rise to the atmospheric structure, which in turn allows for the solution of the equation of radiative transfer, the solution of the rate equations can readily be included in this process. The usual procedure is to construct a model atmosphere in LTE that nearly satisfies radiative equilibrium. At some predetermined level of accuracy, the rate equations are substituted for the Saha-Boltzmann excitation and ionization equations by using the existing structure (electron density and temperature distribution) and radiation field. The resulting population levels are then used to calculate opacities and the atmospheric structure for the next iteration. One may even chose to use an iterative algorithm for the solution of the linear equations, for an initial guess of the LTE populations will probably be quite close to the correct populations for many of the levels that are included. The number of levels of excitation that should be included is somewhat dictated by the problem of interest. Depending on the state of ionization, four levels are usually enough to provide sufficient accuracy. However, some codes routinely employ as many as eight. One criterion of use is to include as many levels as is necessary to reach those whose level populations are adequately given by the Saha-Boltzmann ionization-excitation formula.

Many authors consider the model to be a non-LTE model if hydrogen alone has been treated by means of rate equations while everything else is obtained from the Saha-Boltzmann formula. For the structure of normal stellar atmospheres, this is usually sufficient. However, should specific spectral lines be of interest, one should consider whether the level populations of the element in question should also be determined from a non-LTE calculation. This decision will largely be determined by the conditions under which the line is formed. As a rule of thumb, if the line occurs in the red or infra-red spectral region, consideration should be given to a non-LTE calculation. The hotter the star, the more this consideration becomes imperative.

c. Thermalization Length

Before we turn to the solution of the equation of radiative transfer for lines affected by non-LTE effects, we should an additional concept which helps characterize the physical processes that lead to departures from LTE. It is known as the thermalization length. In LTE all the properties of the gas are determined by the local values of the state variables. However, as soon as radiative processes become important in establishing the populations of the energy levels of the gas, the problem becomes global. Let \(l\) be the mean free path of a photon between absorptions or scatterings and \(\mathscr{L}\) be the mean free path between collisional destructions. If scatterings dominate over collisions, then \(\mathscr{L}\) will not be a "straight line" distance through the atmosphere. Indeed, \(\mathscr{L} \gg l\) if \(\mathrm{A_{ij}} \gg \mathrm{C_{ij}}\). That is, if the probability of radiative de-excitation is very much greater than the probability of collisional de-excitation, then an average photon will have to travel much farther to be destroyed by a collision than by a radiative interaction. However, \(\mathscr{L} \gg l\) if \(\mathrm{C_{ij}} \gg \mathrm{A_{ij}}\). In this instance, all photons that interact radiatively will be destroyed by collisions.

If the flow of photons is dominated by scatterings, then the character of the radiation field will be determined by photons that originate within a sphere of radius \(\mathscr{L}\) rather than \(l\). However, (\mathscr{L}\) should be regarded as an upper limit because many radiative interactions are pure absorptions that result in the thermalization of the photon as surely as any collisional interaction. In the case when \(\mathscr{L} \gg l\), some photons will travel a straight-line distance equal to (\mathscr{L}\), but not many. A better estimate for an average length traveled before the photon is thermalized would include other interactions through the notion of a "random walk". If n is the ratio of radiative to collisional interactions, then a better estimate of the thermalization length would be \[l_{\mathrm{th}}=\ell \sqrt{n}=\ell \sqrt{\mathscr{L} \mid \ell}=\sqrt{\mathscr{L} \ell}\label{15.2.27}\]

If the range of temperature is large over a distance corresponding to the thermalization length \(l_\mathrm{th}\), then the local radiation field will be characterized by a temperature quite different from the local kinetic gas temperature. These departures of the radiation field from the local equilibrium temperature will ultimately force the gas out of thermodynamic equilibrium. Clearly, the greatest variation in temperature within the thermalization sphere will occur as one approaches the boundary of the atmosphere. Thus it is no surprise that these departures increase near the boundary. Let us now turn to the effects of non-LTE on the transfer of radiation.