A Simple Explanation for the Observed Power Law Distribution of Line Intensity in Complex Many-Electron Atoms

It has long been observed that the number of weak lines from many-electron atoms follows a power law distribution of intensity. While computer simulations have reproduced this dependence, its origin has not yet been clarified. Here we report that the combination of two statistical models -- an exponential increase in the level density of many-electron atoms and local thermal equilibrium of the excited state population -- produces a surprisingly simple analytical explanation for this power law dependence. We find that the exponent of the power law is proportional to the electron temperature. This dependence may provide a useful diagnostic tool to extract the temperature of plasmas of complex atoms without the need to assign lines.

It has long been known that the number of weak lines emitted by many-electron atoms in plasmas follows an intensity power law. In 1982 Learner pointed out this law for the first time when measuring emission lines from a hollow cathode lamp containing iron atoms [1]. He observed that the number density of lines with a given intensity I, ρ I (I), exhibits a power law dependence on I [2], He also reported that ρ I (I) in different wavelength regions all follow this power law with the same exponent, indicating an ergodic property of the emission line distribution [1]. This work has stimulated much discussion. A theoretical study by Scheeline showed that this power law does not hold for hydrogen atom spectra [3]. In contrast, the emission spectrum from arsenic, which has a much more complex electronic structure than hydrogen, shows an intensity distribution closer to the power law, but with a different value of the exponent [4]. Bauche-Arnoult and Bauche reported a simulation with a collisional-radiative model for neutral iron atom and demonstrated that the power law dependence is again reproduced [5]. Their exponent was 17-25 % smaller than the Learner's value but the reason was not clarified.
Pain recently reviewed this power law dependence problem and presented a discussion regarding fractal dimension and quantum chaos [6]. According to his discussion, the line strength distribution evaluated under the fully quantum-chaos assumption does not explain Learner's law. As presented in his review [6], as well as in the book [7], the origin of this power law is still not understood, despite almost 40 years passing since the first report.
In this Letter, we present a surprisingly simple explanation of Learner's law. We assume local thermal equilibrium of the excited state population, and an exponential increase in the level density of complex atoms, which has been reported in several many-electron atoms and ions (e.g. [8,9]). Combining these, we show below that the number of levels with a given population follows a power law distribution. An assumption of independently and identically distributed radiative transition rates then directly gives Learner's law in the form where k is Boltzmann's constant and T e is electron temperature in the plasma. 0 is a scale parameter representing the level density growth rate against the excited energy (see Eq. (3) for its definition), which can be estimated either from the experimentally derived energy levels or from ab initio atomic structure calculations [9]. Plasma spectroscopy has been developed from simpler systems, e.g. hydrogen and rare gas atoms. It is known that comparison between intensity ratios of certain emission lines and collisional-radiative models provide us with information about plasma parameters, such as electron temperature and density [10][11][12][13]. This requires correct line identifications and accurate atomic data such as energy levels, oscillator strengths, and collision cross sections. However, accurate atomic data for open-shell atoms is difficult to obtain despite numerous demands for plasma diagnostics with complex atoms, ranging from laser produced plasmas for extreme-ultraviolet light sources [14][15][16], heavy-metal-contaminated fusion plasmas [17,18], to the emissions found after the r-process supernova (kilonova) [19,20]. Our result Eq. (2) suggests an advantage of using intensity statistics for diagnosing plasmas with many-electron atoms, where accurate ab initio simulations of such complex spectra are still difficult with currently available theory and computers.
Next we derive Eq. (2), illustrating our assumptions using Learner's example of neutral iron. Figure 1(a) shows the level density of neutral iron, ρ E (E), the number of levels with given excited energy E. This state density is evaluated from the measured energy levels taken from Atomic Spectral Database by National Institute of Standards and Technology [21]. The observed energy levels are shown by the vertical bars in the figure.
It is well known that the excited level density in the quantum many-body system increases nearly exponentially. One simple but common approximation is [9,23,24] where 0 is an atom-specific energy scale, which depends on the number of active electrons and the number, degeneracy, and distribution of single particle states in the atom [9]. It can be calculated numerically, derived from experimental energy levels, or estimated using combinatorics. Dzuba et al. presented that for open-d-or f -shell atoms the state density follows Eq. (3), at least below the ionization energy [9]. For neutral iron, we find Eq. (3) well represents the level density with 0 ≈ 1.97 ± 0.04 eV as indicated by the solid line in Fig. 1(a). Note that this value is obtained by the maximum-likelihood estimation of the simulated energy levels. Let us assume local thermal equilibrium for the excited state population. The population in state i with energy E i is given as where g i is the statistical weight of the state i (g i = 2J i +1, where J i is the total angular momentum quantum number of state i). This equilibrium is valid in plasmas with high electron density and low electron temperature [11]. By substituting Eq. (3) into Eq. (4), the number of states having the population n ∼ n + dn can be written as, where dE is the energy interval corresponding to dn, the relation of which can be obtained from Eq. (4). We assume in Eq. (6) that the statistical weight is distributed uniformly over the energy and therefore we omit it from the equation. This power law originates from the combination of one exponentially increasing variable and another exponentially decreasing variable. This is a typical mathematical structure responsible for the emergence of power laws [25]. The emission intensity I ij corresponding to the transition i → j, where j is the lower state, is proportional to the upper state population n i , the transition energy cubed ω 3 ij = (E i −E j ) 3 , and the line strength S ij between i and j states. In many-electron atoms with sufficient basis-state mixing, i.e., in quantum-chaotic systems, the probability distribution of S ij can be well approximated as uniform and independent, and modeled using the Porter-Thomas distribution p(S) ∝ 1 √ 2πS0S exp(− S 2S0 ), with a constant S 0 [8,[26][27][28][29]. This approximation is obtained by modeling the Hamiltonian with a Gaussian orthogonal ensemble. As this distribution decays considerably faster than the power law in the large S limit, we can safely approximate that S ij is a constant for all pairs of levels. Therefore, the intensity I ij is approximated as A more detailed and precise discussion can be found in the Supplementary Material [30]. The number of emission lines from state i observed in photon energy range ω ∼ ω + dω is proportional to the number of levels in this energy range, ρ E (E i − ω)dω. By considering the number of emission lines with a given intensity range I ∼ I + dI, we arrive at Eq. (2), where the integration along ω is taken over the observed photon energy range, Ω. Here the variable E is changed to I based on Eqs. (4) and (8). The factor 2 newly appears in the exponent of I compared with Eq. (7).
The exponent in Eq. (2) does not depend on Ω. This is consistent with Learner's observation that the emission line density in different wavelength regions all show the power law dependence with the same exponent [1]. Learner suggested a relation between the exponent and a constant, log 10 √ 2 [1]. In contrast, our work clearly indicates a relation with T e and the atom-specific constant 0 . By comparing the exponents in Eqs. (1) and (2), T e in Learner's experiment is estimated as (1.50 − 1) 0 /2 ≈ 0.49 eV. Although in Ref. [5] it is claimed that T e higher than 0.34 eV is not realistic, T e in hollow cathode discharges reported in literature varies from 0.2 eV to 3 eV depending on the cathode element, filler gas pressure, and discharge current [31,32]. Therefore, 0.49 eV may not be a surprising value for T e in a hollow cathode discharge.
Bauche-Arnoult and Bauche have used T e = 0.34 eV for their simulation [5], which is smaller than 0.49 eV. They obtained −1.392 ± 0.017 for the exponent [33], which is consistently smaller in magnitude than Learner's value. Our above discussion further provides an explanation for one argument in their paper, i.e., the higher the electron temperature, the larger the magnitude of the exponent [5].
We carry out an ab initio simulation of the emission spectrum of neutral iron with the flexible atomic code (FAC) [22]. FAC uses the relativistic Hartree-Fock method to calculate the electronic orbitals and configuration interaction to approximate the electron-electron interaction. The excited state population and the emission line intensity are evaluated by the collisional-radiative model implemented in FAC, where the steady state of population in the plasma is assumed. For the collisionalradiative calculation we consider spontaneous emission, electron-impact excitation, deexcitation, and ionization, as well as auto-ionization of levels above the ionization threshold, as elementary processes in plasmas. These rates are also calculated by FAC.
We assume T e = 0.34 eV and the electron density n e = 10 20 m −3 , similar to Bauche-Arnoult and Bauche [5]. We also perform the simulation with T e = 0.70 eV to observe the T e -dependence of the exponent. Note that in the FAC computations we do not explicitly adopt either of our two assumptions, namely the exponential increase of the state density and the local thermal equilibrium of the population.
The state density of a neutral iron atom computed by FAC is shown in Fig. 1(a) by a blue histogram. It shows a similar exponential dependence to the measured data, NIST ASD. Figure 1(b) shows the excited state population computed by FAC. Although we do not assume local thermal equilibrium, the population follows the exponen-  7) and (2), respectively) with the same electron temperatures used in Fig. 1(b).
tial function. The exponents for T e = 0.34 eV and 0.70 eV cases are estimated by the least-squares method to be 0.32 eV and 0.61 eV, respectively, which are similar to the electron temperature. Note that the slight difference between T e and the exponent in the population is caused by a small violation of the local thermal equilibrium in plasma [11]. Based on approximate electron impact excitation and deexcitation rates, we find that this effective temperature has a weak n e dependence and approaches T e in large n e limit. Even with a smaller electron density, n e = 10 17 m −3 , this temperature is expected to be ≈ 0.7 T e for the T e = 0.34 eV case. Details can be found in the Supplementary Material [30]. The histograms in Fig. 2(a) show the state density ρ n (n) with given population n (but scaled by n to aid visualization). The solid blue and red lines are computed according to Eq. (7) with T e = 0.32 and 0.61 eV, respectively (the same temperatures used in Fig. 1(b)). Their agreement is clear. Figure 2(b) shows the line intensity distribution ρ I (I) in the visible and infrared wavelength range (scaled by I for visualization). The solid lines show Eq. (2) with T e = 0.32 and 0.61 eV for the two cases. This also agrees with the above discussion, particularly in the first three orders studied by Learner. In Eq.
(2), we show that the exponent exclusively depends on T e and 0 , but not on other atomic data such as level energies, transition rates, and collision cross sections. The only value we require, 0 , is known to be accurately calculated with several atomic structure packages [9]. Therefore, Eq. (2) may be useful as a quick diagnostic method for many-electron atom plasmas.
As Eq. (2) is scale-free for I, the power law dependence is not affected by the system's sensitivity (see the Supplementary Material for details [30]). Thus, no system calibration is required for estimation of T e . We only need to know the dominant (in terms of the number of emission lines) atom in the plasma.
We have applied this approach to the emission spectra measured for thorium plasmas. Figure 3 (a) and (b) show spectra measured from a thorium-argon hollow cathode plasma with 75 mA discharge current [34] and 20 mA discharge current [35], respectively, with a 1-m Fourier transform spectrometer. The original data can be downloaded from Kitt Peak National Observatory website [36]. Only the spectra in 510-560 nm wavelength range are shown and analyzed in this work. Dots in each panel show the line centers and peaks detected in the two spectra.
Although not all emission lines in these spectra have been identified, we assume that most of the lines are from neutral thorium. We compute ρ I (I) from all line intensities in the wavelength range (Fig. 3(c)). The two histograms generated from the spectra show a power-law distribution. ρ I (I) for the higher current discharge shows a steeper slope.
We estimate the exponent of these distributions using the maximum-likelihood method. The optimized distributions are shown by solid lines (and their 2-σ uncertainty by colored bands) in Fig. 3 (c). They fit both histograms. Estimated values of the exponent are 1.71±0.03 and 1.64 ± 0.03 for the 75 mA and 20 mA disharges, respectively. From Eq. (2) and the value 0 ≈ 0.68 eV for neutral thorium by Dzuba et al. [9], electron temperatures for these plasmas are estimated as 0.24 ± 0.01 eV and 0.21 ± 0.01 eV, respectively. A higher T e value is estimated for the higher current discharge. Although the positive current dependence of the temperature is not trivial [31], this dependence qualitatively supports our model. Because there are no radiative rates reported for neutral thorium, it is difficult to estimate T e for this plasma by conventional methods. To our knowledge, the above procedure is the only one available to estimate T e for thorium plasmas.
Although there are significant demands to diagnose plasmas with many-electron atoms, quantitative comparison with an ab initio computer simulation model is not yet accurate enough, because of the unavailability of accurate atomic data. Our result suggests a possibility of plasma diagnostics that requires only the energy level statistics and the emission intensity statistics. Although the validity of the local thermal equilibrium assumption should be investigated further, this may open the door to a new statistical plasma spectroscopy.
In summary, we have presented a simple explanation of Learner's law, where the histogram of the emission line intensities from many-electron atoms follows a power law. We observed that the exponent is analytically represented with T e and 0 . A similar discussion should also be applicable to other fermionic many-body systems as long as the two assumptions are satisfied. Although as yet there are no reports about the emission statistics except for many-electron atoms, it is interesting to investigate other systems, such as heavy nuclei. This is left for future study.  In the main text, we used "∝" in most of the equations and ignored constant factors (e.g., Eqs. 2, 3, and 6) to simplify the discussion. Furthermore, we made a rather drastic approximation in Eq. (8): the constant radiative transition rate. In this Supplemental Material, we present explicit formulae and explain the approximation in more detail. As we will see, even if we consider the probability distribution of the transition rate, we arrive at the same result.
Let us redefine explicitly the level density ρ E (E) and the population n i at state i having the excited energy E i as follows, and where ρ 0 and n 0 are constants. g is the averaged statistical weight for all the state, which we assumed constant over all the levels in the main text. The distribution of g i does not affect the result as long as the distribution is uniform and independent of energy. The number of states having the population n ∼ n + dn can be written as Let us consider the number of emission lines found in the photon energy range of ω ∼ ω + dω from level i with 0 ≤ ω ≤ E i . Because of the finite wavefunction mixing, the number of emission lines equals to the number of levels within The number of emission lines from the excited states existing within the excited energy E ∼ E + dE is The intensity of the emission line corresponding to the transition from the state i to state j, I ij , is written as I ij = A ij n i , where A ij is the radiative transition rate from state i to state j. A ij relates to the line strength S ij , A ij = γω 3 ij S ij , with ω ij = E i − E j and γ = 16π 3 e 2 3ν0h 4 c 3 , with e the elementary charge, ν 0 vacuum permittivity, h is Planck's constant, and c light speed.
In order to simplify the discussion, for the moment let us assume A ij has a solid relation A ij = γω 3 ij S 0 for all i and j pairs with a constant S 0 , and relax this assumption later. In this case the intensity only depends on the population of state i. The number of lines with given intensity I c under this constant S 0 assumption is where dI c is the intensity range that corresponds to the energy range dE, which can be found from Eq. (S2) and I c = ω 3 γS 0 n. Note here already that the dependence of I c in the distribution follows the final distribution (Eq. (2)). Now we relax the constant line strength assumption and consider the stochastic nature of S ij . Lets be the fluctuation in S, S =sS 0 . The probability distribution ofs may be written as ps = 1 √ 2πs exp(−s/2), which is a Porter-Thomas distribution with mean 1 [8,[26][27][28]. As we will see later, the actual shape of ps is not important as long as it decays exponentially in the larges limit. More importantly, we here assume that ps is independent and identical for all the transitions. With this assumption, the intensity I becomes a product of two independent random variables, I =sI c .
Let us consider only the probability distribution of I c , rather than the number of lines. Since the distribution of I c follows the power law, the probability distribution of I c is written as a Pareto distribution, where α = 2kT e / 0 and I min is the minimum value of the intensity. We will take a limit of I min → 0 later but in order to make the distribution integrable, we now assume this value is sufficiently small. The probability distribution of the intensity is written as If ps(s) decays faster thans −α−1 in larges limit (such as the exponential decay) and I I min , the integral in the right hand side converges and does not depend on I. With sufficiently small I min (which corresponds to sufficiently large E for the upper state), p(I|ω) becomes the power law distribution except for the very small I region. This distribution actually has the same shape as Eq. (S8) where we assume that S is constant.
Let us consider the sensitivity of the observation system, ξ(ω). The density distribution of the observed intensity I = ξ(ω)I is The system's sensitivity only affects the scale of the distribution, but not the exponent of the power law. A histogram may be computed from an emission spectrum observed in finite wavelength range, Ω. It corresponds to the integration of Eq. (S14) over the observation wavelength range, which results in the same power law. Therefore, no system calibration is necessary to compute the intensity histogram. This ergodic property essentially comes from the fact that the power law is scale invariant.

BIAS IN THE Te ESTIMATION
The T e estimation method newly proposed in this work is based on the local thermal equilibrium assumption. As a slight temperature difference can be found between the simulated and estimated T e values in the main text, this assumption is not always satisfied perfectly. Here, we present a simple analytical model to predict the estimation bias.
We consider plasmas with weak radiation field, where the photo excitation of atoms is negligible. We further neglect ionization and recombination for the sake of the simplicity, and only consider electron collision excitation and deexcitation, as well as radiative deexcitation. Let us assume that the excited level population of an atom is approximated by Boltzmann's distribution with an effective temperature T eff , where T eff may be different from T e . We seek the value of T eff where the population balance is established under electron-collisional excitations, deexcitations and radiative deexcitations. Let r ij (T e ) be the collisional excitation rate from i to j states. Based on the detailed balance principle, the inverse process, i.e., the collisional deexcitation rate from j to i state can be written as follows, where ω ij is the energy difference between i and j states. Here, we neglect the statistical weight difference in states i and j as in the previous section. The net population influx from j to i states is written as γ e ji = −n i r ij (T e ) + n j r ji (T e ) (S18) The total population flux into state i is Γ e i = j γ e ji . We approximate it by continous integration taking the level density ρ n into account, where We use an approximate form of r ij , with β = where m e the electron mass, a 0 Bohr radius, E H Rydberg unit of energy, and ξ ≈ 0.6 the integrated gaunt factor [37]. Here, we substituted the relation of the oscillator strength f ij = 16π 2 3e 2 meh 2 c 2 g ω ij S ij to the original formula found in Ref. [37]. Substituting Eqs. (S21) and (S19) into Eq. (S20), Γ e i can be written as where x = 1 + 0 kTe − 0 kT eff . Note that we neglected the variation of S ij and assume it as constant S 0 . E i 0 is also assumed in the above derivation.
The population influx by radiative transition into state i can be written as Γ rad indicates the summation over the state j having with E j > E i (E j < E i ). We approximate it by a continuous integration as follows: Although its solution is not analytically tractable, it can be easily solved numerically. We note that the dependence on E i is not significant if E i is a few times larger than 0 . We assume E i ≈ 4 0 in the following numerical evaluation, which corresponds to the first ionization energy of neutral iron χ Fe = 7.90 eV with 0 = 1.97 eV. It should be also noted that T eff 0 /2 is required to keep the approximations valid that are used for deriving Eq. (S25). The n e dependence of T eff with 0 = 1.97 eV, which is the same value used in the main text for neutral iron, is shown in Fig. S1 (a). T eff approaches to T e in the large n e limit, i.e., when the relative importance of the radiative decay becomes negligible. T eff decreases in a lower n e plasma. The vertical bars in the figure show threshold densities, at which T eff = 0.9 T e . Lower T e values give lower threshold densities, however even below the threshold, the n e dependence of T eff is weak.
We also carry out ab initio simulations with FAC under several n e and T e conditions. T eff is estimated from the simulated excited state population. The estimated T eff values are shown by markers in Fig. S1 (a). Note that FAC utilizes the distorted wave approximation for simulating the electron impact excitation rate, which takes the mixing of wavefunctions into account, i.e., this approximation is more robust than Eq. (S21).
A qualitative agreement between the analytical model and FAC is clear, especially in the threshold density. This indicates the validity of the above approximations. Although there is still a slight discrepancy in low n e region, the agreement is surprising because only 0 is used as an atomic data in the analytical model. We may even be able to correct the estimation bias on T e based on this model if we have a rough value of n e in the plasma.
In Fig. S1 (b), we show T eff values for neutral thorium with 0 = 0.68 eV. For a realistic electron density in the plasma (n e = 10 18 -10 20 m −3 ) and for T eff ≈ 0.2 eV as estimated in the main text, T eff is within 30% of T e . This small n e dependence indicates the sensitivity for the temperature measurement. Since the emissivity is roughly proportional to n e , if n e varies over the measurement region and time but T e does not, the exponent will correctly represent the temperature. This is a unique feature of this plasma diagnostic technique with many-electron atom emission, while many of the existing methods must assume a single density and temperature over the measurement range and duration.

CRITERION OF LOCAL THERMAL EQUILIBRIUM FOR MANY-ELECTRON ATOMS
For hydrogen-like atomic ions, Griem has proposed a validity criterion of local thermal equilibrium [38], n e 7 × 10 24 p −17/2 z 7 kT e where z the nuclear charge, p the principal quantum number, and χ H = 13.6 eV the ionization energy of neutral hydrogen. In this section, we derive a similar criterion for many-electron atoms from Eq. (S25).
We define the boundary density where T eff 0.9T e . With the assumption 0.1 kT eff / 0 < 1/2, Eq. (S25) may be reduced to For the T e = 0.7 eV case, this boundary density is 2 × 10 21 m −3 , which is similar to the numerically derived one (5 × 10 20 m −3 , shown by vertical bars in Fig. S1). This criterion has a different T e dependence to the Griem criterion Eq. (S26). However we can make a connection with the Griem criterion by noting that the emitter charge state distribution will itself change with temperature. In this case 0 changes with T e , but 0 kTe and z 2 χH kTe may be approximately constant. Here z is promoted to an "effective charge" such that z 2 χ H is the ionization energy of the emitter. From the case with neutral iron and T e = 0.7 eV, 0 kTe ≈ 3 and z 2 χH kTe ≈ 10. Since the chaotic behavior in many-electron atoms starts around E 0 , the corresponding effective principal quantum number for the chaotic states in neutral iron may be p ≈ χ H /(χ Fe − 0 ) = 1.5. With this effective principal quantum number, Griem's boundary density for many-electron atoms becomes with T e in eV. Both the boundary densities scale as T 7/2 e consistently, but they should be understood as providing only a rough order-of-magnitude estimate.
By contrast, the consistency with the criterion proposed by McWhirter [39], which is based on the maximum level spacing, is not obvious, because our result does not depend on the absolute level spacing ≈ 1/ρ 0 , i.e., ρ 0 is canceled out in Eq. (S25). More precise and quantitative discussion of the validity condition for local thermal equilibrium in many-electron atoms is left to future study.
Let us consider the validity of local thermal equilibrium in some typical many-electron atom plasmas. In laser produced plasmas for ultraviolet light sources, the temperature is typically ≈ 10 2 eV while the density varies in the range 10 25 − 10 27 m −3 , depending on the incident laser wavelength [40,41]. This parameter range is around the boundary Eq. (S29). Although many parameters are unknown for the laser produced plasmas, including the actual value of 0 for the dominant emitter ion, the density dependence on the spectral shape (showing a broader spectrum from more dense plasmas) has been observed [40].
The divertor plasmas of tokamak nuclear fusion reactors is also close to the boundary (T e ≈ 10 0 eV, n e ≈ 10 20 m −3 ) [42,43]. On the other hand, parameters of tokamak core plasmas stay typically very far from the threshold, n e ≈ 10 20 m −3 and T e ≈ 10 4 eV; local thermal equilibrium is not expected for the core plasma.