Tutorial: Exact Solutions for the Populations of the n-level Ion

Essentially three parameters fully determine the physical nature of an ionized nebula: the electron temperature (T_e), the electron density (n_e), and the level populations of the ions (N_i). From the latter, ratios of emission rates and emission line intensities can be calculated (see, e.g., Aller 1984; Osterbrock & Ferland 2006). The 3-parameter family forms a closed set whereby knowledge of any two of them allows for a determination of the third, since

In principle, once this equation is solved, all other secondary physical quantities such as ionization parameters, ionic abundances, and effective temperatures can then be calculated. Prior to the development of matrix solutions to the problem using the master equation method of Martin et al. (1996) or numerical approximations (Mendoza 1983), the precise functional form of f had remained elusive and unknown to astronomers and atomic physicists.

In parallel, great progress has been made in observational astrophysics, with the development of empirical tools to estimate the nebular values of n_e (Menzel et al. 1941; Copetti & Writzl, 2002; Shaw & Dufour, 1995), T_e (Shaw & Dufour, 1995; Aller et al. 1949; Seaton 1954; Spitzer 1948) and N_i (Seaton 1975). In particular, certain ratios of forbidden lines emitted by the p², p³, and p⁴ configuration ions present in the nebular gas allow the integrated electron density and electron temperature to be estimated point-to-point across projected images, while sums of forbidden lines emitted by different stages of ionization of the gas have been found to weakly correlate with ionic abundances (Shaw & Dufour, 1995; Pagel et al. 1979; Taylor & Díaz 2007).

The early work in this field was performed by Menzel, Hebb, Aller, Spitzer, Seaton, and Osterbrock (Menzel et al. 1941; Aller et al. 1949; Spitzer, 1948; Hebb & Menzel, 1940; Seaton & Osterbrock, 1957; Seaton 1960) who obtained initial estimates of the level populations of 3-level ions. An exact analytic solution for the 3-level ion was finally worked out just over 30 years ago by Seaton (1975). The computer age then allowed for the development of various higher order numerical approximations. For example, a first-order numerical approximation for the 5-level ion is currently used in nebular analysis software, such as the code FIVEL by De Robertis et al. (1987), the 3D Monte Carlo photoionization code MOCASSIN by Ercolano et al. (2003), or state of the art codes that use Breit-Pauli R-matrices to approximate ions containaing even higher stages of ionization. Rodriguez (2002) has produced an approximate 34-level model of [Fe III] and Pelan & Berrington (2001) have approximated the [Fe IV] ion to 180 levels. Despite these computational achievements, only the master equation method (Martin et al. 1996) has been advanced as a way to calculate exact analytical solutions. In this tutorial, we will guide the reader through the method using a simple matrix prescription.

The theoretical context for the physics of collisionally-excited ions is the thermal equilibrium present in ionized gases. The temperature in a static nebula is fixed by the equilibrium between heating by photoionization and cooling by recombination, free-free radiation (bremsstrahlung), and line radiation. Menzel and coworkers (Menzel et al. 1941; Hebb & Menzel 1940) first set up the equation of energy conservation in ionized nebulae, whereby heating due to the energy absorbed by photoionization of the gas (G) is balanced by cooling due to the energy liberated in capture and subsequent recombination events (L_R), free-free Bremsstrahlung emission (L_ff), and energy emitted in collisionally-excited radiative decays (L_C),

Closed forms for G, L_R, and L_ff can be found in standard texts on the physics of ionized gaseous nebula (see, e.g., Aller 1984; or Osterbrock & Ferland, 2006). The aim of this article is to show how to calculate a closed form for L_C

the cooling due to the emission of collisionally-excited lines (CELs) with level populations N_i, transition probabilities A_ij, de-excitation rate coefficients q_ij, and energy level potential differences hν_ij. In this tutorial, we will focus on the low-density regime below the critical density , where collisional de-excitation of level i is negligible—the thermal regime of the vast majority of planetary nebulae.

In many gaseous nebulae, mechanical or magnetic energy from hydromagnetic waves is dissipated in the gas and there can exist regions where the mean temperature is raised such that atomic levels are excited through collisions with electrons. However, in the low density plasmas of many nebulae and ionized hydrogen (H II) regions, the most frequently observed optical lines come from the forbidden lines that violate the Laporte Parity Rule (Aller 1984). For these lines, we need to know the transition probabilities (A_ij) and collision strengths (Ω_ij) for the excitation of metastable levels. Forbidden lines cannot be observed in the laboratory and so we have to rely entirely on theory to determine A_ij and Ω_ij (see, e.g., Mendoza 1983). In what follows, we use the matrix master equation to rigorously solve for the populations of the n-level ion that give rise to these emitted forbidden lines.

One should start with the formal general equations of statisitcal equilibrium for an ion with excited levels (i) above ground (Aller 1984, eqs. 5–32), which can be written (Appendix A)

and which are subject to the total ion density condition that the sum overall level populations (i.e., over all stages of ionization) is equal to the total number density (N) of the ions [cm^-3]

Here we have adopted the excitation potential convention E_i > E_j such that the transition i → j corresponds to de-excitation and A_ij is the transition probability per second, the Seaton variable (Seaton 1954) for electron temperature T_e [K] and electron density n_e [cm^-3], Ω_ji is the velocity-averaged collision strength (Seaton 1968), the statistical weights are ω_k = (2k + 1), and the atomic constant K = 8.629 × 10^-6.

The full statistical equilibrium is therefore described by a set of (n - 1) simultaneous equations in n unknowns for the excited levels above the ground level (i = 1), with closure ensured through the addition of the total ion density condition. In matrix form,

such that the elements α_ij are equal to the coefficients of N_i. The first row with α_1j = 1 is the total ion density condition. Introducing the matrix of coefficients (α_ij), the vector y of level populations (N_i), and the vector b for the righthand side, the matrix master equation can be written simply and succinctly as

Provided that (which is always the case since ions have nonnegative activation levels above ground), then we can use symbolic manipulation software such as MATLAB or Mathematica to obtain and hence the level population vector

Since b₁ = N and b_i>1 = 0, the sought-after level populations are identically given by

It is at this point that we can identify the function f in equation (1). We know from the above definitions that the population of the ith level N_i is given by

If we use now introduce the matrix of cofactors of , then its transpose (the adjugate matrix) allows for an algebraic solution for the matrix inverse

such that the level populations will have the general form

or equivalently,

as implied by equation (1).

Equations (4)–(9) allow the astrophysicist to model (exactly) an ion of as many levels as required.

3.1. The 5-Level Ion

Ions having p², p³, and p⁴ electron configurations all have five low-lying energy levels. For such ions, collisional and radiative transitions can occur between any of the levels, and excitation and de-excitation cross sections as well as collision strengths exist between all pairs of levels. A central assumption that is often made is that line emission from these five levels alone provides a good approximation to the expected emission lines of the full n-level ion. The justification (see, e.g., Osterbrock & Ferland 2006) is that higher levels in these ions are not significantly populated through collisions, recombinations, or other mechanisms. Although this is true for weakly ionized, low-density nebular plasmas such as planetary nebulae and HII regions, active galactic nuclei (AGNs) and supernovae (SNs) are much more strongly ionized, meaning that higher level transitions are likely to be significant to their emission line spectra—as evidenced by the tremendous numerical modeling work of the IRON and OPACITY projects.

For a 5-level ion in a steady state with E₅ > E₄ > E₃ > E₂ > E₁, the total number density condition for each ion species is given by N₁ + N₂ + N₃ + N₄ + N₅ = N and the equations of statistical equilibrium (4) give rise to the following four exact level population equations for the excited levels above ground

The elements α_ij of are then equal to the coefficients of the level populations N_i and are listed in Appendix B. We have deliberately not inserted the values of A_ji, Ω_ji, and _ij up until now as they are estimated values that will evolve with the accuracy of quantum mechanical calculations currently being performed by the IRON and OPACITY projects.

The matrix equation having the coefficients α_ij is exact for the 5-level ion. At this stage, the problem is reduced to finding the matrix inverse. This can be done numerically as in Martin et al. (1996) or using singular value decomposition. With the advent of symbolic mathematics software such as MATLAB and MAPLE, it is also possible to find the inverse matrix analytically so as to obtain exact expressions for the level populations. For example, the full solution includes terms up to (Kx)⁴, and while cumbersome and too large to reproduce here, is exact. However, in Appendix C we provide a general public license (GPL) MATLAB code for generation of the symbolically-obtained analytic solutions for comparison with exact numerical solutions. In general, an n-level ion will contain terms up to (Kx)^n-1. The inverse matrix solution for the 3-level ion includes terms only up to (Kx)² and we show in § 3.2 how the known analytic solution (Seaton 1975) can be reproduced.

3.2. The 3-Level Ion

For ions of type p² such as the [O III] ion shown in Figure 1, suppression of all terms with indices equal to 4 and 5 reduces the general solution for the 5-level ion to that of the 3-level ion. This amounts to removing the fine splitting physics of the ³ P level. As a consequence, the excitation potentials E₂ and E₃ of the 3-level ion correspond identically with E₄ and E₅ of the 5-level ion while E₁ = 0 eV is the ground level. Similarly, a suitable suppression of indices allows the 5-level representations of the p³ and p⁴ ions to be reduced to their exact 3-level counterparts. In every case, the following exact level populations for the 3-level ion are obtained

with denominator D,

In order to compare the exact matrix master equation solution for the n = 3 level ion with the known analytic solution of Seaton (1975), we consider the ratio of the level populations N₃ and N₂ (although nowadays the inverse N₂/N₃ for the ratio of nebular to auroral lines is customary)

Extracting the term (Ω₁₃₃₁) from the numerator and (Ω₁₂₂₁) from the denominator and canceling like terms we obtain

where

This is the same result as that obtained by Seaton for the exact 3-level ion (Seaton 1975, eq. 1.11).

Zoom In Zoom Out Reset image size

By considering the dimensions of the physical variables (there are 7 metric units: mass M[Kg], length L[m], time T[s], temperature θ[K], electrical current I[A], concentration N[mol], and light intensity J[cd]), a further theoretical check is provided by the requirement for dimensional homogeneity which for a set of physical variables {x_i} requires that in order to satisfy causality (Taylor et al. 2008). The dimensions of the various physical quantities and parameters presented in our derivation are

such that

All equations obtained from the matrix master equation approach were thoroughly checked for dimensional homogeneity with these dimensional relations and no inconsistencies were found.

3.3. Empirical Checks

In order to compare the algebraically-obtained for the level populations of the 5-level ion with those obtained from observational measurements and numerical computer simulations, it is necessary to know the values of the measured emission line intensities of an ion I(i,k) and the average nebular electron density n_e and electron temperature T_e, as well as the quantum-mechanically calculated values of A_ji, _ij(T_e), Ω_ji(T_e), and x(n_e,T_e). The line intensities (being the only direct observable) are proportional to the line emission rate j(i,k) of line photons resulting from downward transitions i → k,

where X^l is the ion species; e.g., O⁺² and ΔE = hν_ik corresponds to the energy difference between excitation potentials E_i and E_k. Hence, a ratio of line intensities will involve only the ratio of level populations N_i and ionic constants A_ik and ΔE_ik. So, for example, the ratio of auroral to nebular lines of [O III] which we denote

for the 3-level ion is given by

and for the 5-level ion is

Therefore, by measuring observationally, it is possible to compare its value with that obtained using the level populations (9) provided by the matrix master equation method using the values of A_ji, _ij(T_e), Ω_ji(T_e), and x(n_e,T_e) evaluated at the measured value of electron density and temperature.

3.4. Spectrophotometry of the Spherical Planetary Nebula A39

The spherical planetary nebula A39 was discovered by George Abell in his 1957 survey of the Southern Hemisphere (Abell 1966) and has recently been subject to thorough spectrophotometry (Jacoby et al. 2000). From these detailed measurements made by George Jacoby and coworkers at Kitt Peak, the reddening-corrected line ratio averaged over the whole nebula was measured to be (Jacoby et al. 2000)

with all line intensities I(λ) relative to I(H_β). The electron density and electron temperature were estimated to be n_e ≈ 30 cm^-3 T_e ≈ 15400 K (Jacoby et al. 2000). Inserting the most up-to-date values of the atomic and ionic constants A_ij, Ω_i,j, and _ji (see Osterbrock & Ferland 2006 and references therein) calculated at n_e = 30 cm^-3 and T_e = 15400 K, we obtain for the exact 5-level [O III] ion and for the Seaton 3-level [O III] ion. Furthermore, the first-order numerical approximation to the 5-level [O III] ion provided by the code TEMDEN (De Robertis et al. 1987) gives at this electron density and temperature; all well within the observational standard error.

3.5. Three-Dimensional Photoionization Modeling of A39

Using the 3D Monte Carlo photoionization code MOCASSIN (Ercolano et al. 2003) and a best-fit model to the observed spectrum with a 3-shell density profile (Figure 2), one obtains the estimate at an average nebular electron density 〈n_e〉 = 11.94 cm^-3 and temperature 〈T_e〉 = 15896 K. Once again, inserting the values of the atomic and ionic constants A_ij, Ω_i,j, and _ji calculated at n_e = 11.94 cm^-3 and T_e = 15896 K, we obtained the values for the exact 5-level [O III] ion, for the Seaton 3-level [O III] ion, and with TEMDEN. Incidentally, a rerun of MOCASSIN using the value of averaged electron density 〈n_e〉 = 30 cm^-3, suggested from the observational study of Jacoby and coworkers, yielded a slightly lower average nebular temperature 〈T_e〉 = 14997 K. In this particular case, MOCASSIN obtained the much higher value , in agreement also with high values for the exact 5-level ion , the Seaton 3-level ion , and TEMDEN . To check that such an increase is to be expected at lower electron temperatures, we calculated the value of using the tabulated values for [O III] (Lennon & Burke 1994) of A_ij, Ω_i,j, and _ji at an electron density n_e = 30 cm^-3 and temperature T_e = 10000 K. Here, we found that the exact 5-level ion gives , the Seaton 3-level ion gives , and the approximate 5-level ion calculation of TEMDEN gives ; all dramatically larger and much more dispersed. The level populations, as expected, are very sensitive to electron temperature. All of these results are collected and presented in Table 1 for ease of comparison.

Zoom In Zoom Out Reset image size

At electron densities of n_e = 30 cm^-3, all of the results are within the observational standard error apart from the low-temperature case (T_e = 10000) where, even here, there is consistency between the exact 5-level ion, the Seaton 3-level ion, and the numerical approximation to the 5-level ion provided by TEMDEN. In the case of the best-fit 3D photoionization model with MOCASSIN using a 3-shell density profile having an average electron density of 〈n_e〉 = 11.94 cm^-3, the results are more dispersed with the code producing a value of closer to the observations. This may be due to uncertainties in the estimated values of the atomic constants A_ij, Ω_ij, and _ji. What is interesting is that the value of from the exact 5-level ion is consistently higher than both the exact Seaton 3-level ion and the numerical approximation to the 5-level ion of TEMDEN, suggesting that the extra terms are not negligible, particularly at lower electron temperatures.

This tutorial presents the master equation method for the calculation of the populations of the n-level ion from the general equations of statistical equilibrium. The matrix formulation presented here, together with the algebraic power offered by symbolic manipulation software, means that it is possible to obtain an exact solution for the populations of ions of an arbitrary number of levels undergoing collisional excitation and de-excitation for comparison with other exact numerical inversions. This is particularly important in studies of strongly ionized ions in high-energy objects such as AGNs or SNs.

The three-fold comparison of theory with observations, simulations, and numerical approximations gives consistent results and we invite the reader to extend the application to other ions (besides [O iii] considered here). The setting up of the n-level problem is straightforward and we hope that these results, in conjuction with the supplementary MATLAB code, will help make the incorporation of higher order transitions into existing astrophysics software easier for the study of very high-ionization objects such as AGNs and SNs. In addition, the interested reader may also like to take up the challenge of including more complex atomic processes such as fluorescence and recombination mechanisms.

M. T. would like to thank Luis Colina Robledo for kindly making available a copy of Aller's 1984 book, Physics of Thermal Gaseous Nebulae, José Cernicharo for originally facilitating this work, the Greek National Scholarships Foundation (IKY) for their financial support, and to the members of DAMIR-CSIC (Madrid) and ISARS-NOA (Athens) for their hospitality. This work is partially funded by project AYA2007-67965-C03-02 of the Spanish MCINN. This work is dedicated to the memory of Donald E. Osterbrock, Bernard E. J. Pagel, and Michael J. Seaton who have done so much to advance this important field.

• Figure 1
• Figure 2

Assuming that level populations N_i are determined by spontaneous emission (with transition probability A_ij s^-1), then the formal general equations of statistical equilibrium for an ion with excited levels (i) above ground having excitation rate coefficients (q_ji) and de-excitation rate coefficients (q_ij) to the other levels j ≠ i are given by Menzel et al. (1941: eqs. 5–32)

subject to the total ion density condition that the sum over all level populations (i.e., over all stages of ionization) is equal to the total number density (N) of the ions [cm^-3]

Here we have adopted the excitation potential convention E_i > E_j such that the transition i → j corresponds to de-excitation. The next step is to note that q_ij times the electron density is equal to the collisional excitation rate c_ij [s^-1]

where K = 8.629 × 10^-6, the Seaton variable (Copetti & Writzl 2002) for electron temperature T_e [K] and electron density n_e [cm^-3], ω_i = (2i + 1) are the statistical weights, and where Ω_ji is the velocity-averaged collision strength (Shaw & Dufour 1995)

for colliding electrons having initial kinetic energy E = mu²/2. The collision strengths must be calculated quantum mechanically, and consist, in general, of a part that varies slowly with energy superimposed with resonance contributions that vary rapidly. The fact that integration over a broad Maxwellian distribution of electron energies tends to dramatically smooth out the variations means that Ω_ji is now known to be fairly insensitive to temperature (Aller et al. 1949). Similarly, the collisional deexcitation rate

is related to the collisional excitation rate through

with _ij = e^{-(E_i-E_j)/k_B T_e}, and where E_i and E_j are the excitation potentials of the relevant levels. The radiative transition probabilities A_ij are independent constants being inversely proportional to the occupancy lifetimes (and hence the energy) of the upper levels. Replacing all n_e q_ji and n_e q_ij terms in the equations of thermal equilibrium for the excited levels by equations (A1)–(A5) for the collisional excitation and de-excitation rates, and dividing through by Kx( ≠ 0), we obtain the equation of statistical equilibrium (4) used in the main text,

From the four level population equations for the excited levels above ground (eq. [14]), the elements α_ij of are equal to the coefficients of the level populations N_i. For the 5-level ion these are:

Here, we provide a general public license (GPL) MATLAB code for the generation of the matrix master equation solution of n-level ions using symbolic mathematics inversions. The code includes relevant atomic data for [O III] and, using the prescription given by equations (4)–(9), can be easily adapted for an ion of n-levels or for other ion species.

% PROGRAM: n_level_ion.m% % m-file containing the symbolic computations for calculation% of the populations of the n-level ion with full matrix components% provided for the exact 5-level ion and the Seaton 3-level ion.% % Written by Michael Taylor 09 October 2007% % Code distributed under the General Public License (GPL)syms O12 O13 O23 E23; syms O21 O31 O41 O51 O32 O42 O52 O43 O53 O54;syms A21 A31 A41 A51 A32 A42 A52 A43 A53 A54;syms E21 E31 E41 E51 E32 E42 E52 E43 E53 E54;syms K x;syms A y b I N;syms Level N_total;% Physical constants (erg, Angstrom) ang=1e-10;erg=1e-07;ev=1.60217733e-19;c=2.99792458e+08/ang;h=6.62607554e-34/erg;K_b=1.380658e-23;% 5-level [OIII]: level energies above ground (eV) (NIST) L1=0.0;L2=0.0140323;L3=0.0379607;L4=2.513565;L5=5.354349;L21=L2-L1;L31=L3-L1;L41=L4-L1;L51=L5-L1;L32=L3-L2;L42=L4-L2;L52=L5-L2;L43=L4-L3;L53=L5-L3;L54=L5-L4;% 3_level ion: Exact solution by symbolic manipulation % b=[N,0,0]';% A=[1,1,1;% O12*E21/3,-A21/(K*x)-(O12+O32)/5,A32/(K*x)+O32*E23/7;% O13*E31/3,O23*E32/5,-(A31+A32)/(K*x)-(O13+O23)/7];% Level=inv(A)*b;% N_total=Level(1)+Level(2)+Level(3);% simple(Level);% pretty(ans);% 5_level ion: Exact solution by symbolic manipulation b=[N,0,0,0,0]';A=[1,1,1,1,1 ; O21*E21/3,-A21/(K*x)-(O21+O32+O42+O52)/5,A32/(K*x)+O32*E32/7,A42/(K*x) +O42*E42/9,A52/(K*x)+O52*E52/11; O31*E31/3,O32*E32/5,-(A31+A32)/(K*x)-(O31+O32+O43+O53)/7,A43/(K*x) +O43*E43/9,A53/(K*x)+O53*E53/11; O41*E41/3,O42*E42/5,O43*E43/7,-(A41+A42+A43)/(K*x)-(O41+O42+O43+O54)/9, A54/(K*x)+O54*E54/11; O51*E51/3,O52*E52/5,O53*E53/7,O54*E54/9,-(A51+A52+A53+A54)/(K*x) -(O51+O52+O53+O54)/11];Level=inv(A)*b;N_total=Level(1)+Level(2)+Level(3)+Level(4)+Level(5);simple(Level);pretty(ans);% Physical conditions for Table 1 ne=30.0;% ne=11.94;Te=15400;% Te=14997;% Te=15896;% Te=10000;% Te=15848;N=1;x=ne/(Te^(.5));K=8.629e-06;E21=exp(-L21/(K_b*Te));E31=exp(-L31/(K_b*Te));E41=exp(-L41/(K_b*Te));E51=exp(-L51/(K_b*Te));E32=exp(-L32/(K_b*Te));E42=exp(-L42/(K_b*Te));E52=exp(-L52/(K_b*Te));E43=exp(-L43/(K_b*Te));E53=exp(-L53/(K_b*Te));E54=exp(-L54/(K_b*Te));% Einstein rate coefficients (NIST) A21=2.61e-05;A31=3.17e-11;A41=1.690e-06;A51=0.0;A32=9.76e-05;A42=6.995e-03;A52=2.268e-01;A43=2.041e-02;A53=6.091e-04;A54=1.561;% Collision cross-sections': numerical approximation formulae (TEMDEN) t4=Te*1.0e-4;C1=1.835+t4*(0.3981-t4*0.06);C2=0.2127+t4*(0.0767-0.013*t4);O21=0.4825+t4*(0.0806-t4*0.022);O31=0.2397+t4*(0.0381-t4*0.007);O41=C1/9; O51=C2/9;O32=1.1325+t4*(0.203-t4*0.05);O42=C1/3;O52=C2/3;O43=C1*5/9;O53=C2*5/9;O54=0.3763+t4*(0.3375-t4*0.105);% 5-level ion calculation n1=subs(Level(1));n2=subs(Level(2));n3=subs(Level(3));n4=subs(Level(4));n5=subs(Level(5));R_OIII=(n4*A42*L42+n4*A43*L43)/(n5*A54*L54);% END: n_level_ion.m