Charge State Distributions in Dense Plasmas

Charge state distributions in hot, dense plasmas are a key ingredient in the calculation of spectral quantities like the opacity. However, they are challenging to calculate, as models like Saha-Boltzmann become unreliable for dense, quantum plasmas. Here we present a new variational model for the charge state distribution, along with a simple model for the energy of the configurations that includes the orbital relaxation effect. Comparison with other methods reveals generally good agreement with average atom based calculations, the breakdown of the Saha-Boltzmann method, and mixed agreement with a chemical model. We conclude that the new model gives a relatively inexpensive, but reasonably high fidelity method of calculating the charge state distribution in hot dense plasmas, in local thermodynamic equilibrium.

Finite-temperature density functional theory (DFT) based models of hot dense plasmas have proven to be very useful for calculating material properties such as conductivity and equation of state [16][17][18].They can also be used to calculate opacity [18][19][20][21].However, due to the Fermi-Dirac occupation of states and the use of Kohn-Sham orbitals, such calculations miss the important physics that arises from considering the plasma to be composed of individual ions characterized by configurations with integer occupation numbers.This physical picture naturally includes the concept of a charge state distribution, which describes the total number of ions in each ionization stage, i.e. those ions that possess the same integer number of bound electrons.Furthermore, this picture allows for the calculation of many more lines (radiative transitions) that appear in the measured spectra.
Several models have been presented that combine DFT-like treatments of free electrons with an integeroccupation treatment of bound orbitals [22][23][24][25][26] to give the energies of the configurations.Using these energies, a population model is then used to determine which of the configurations exist in the plasma and with what probability.Alternative to these models are approaches that use isolated-ion electronic structure calculations, and then introduce the plasma effects, such as bound state pressure ionization, via a secondary model [27].
Here we present a new variational method, which uses ions screened by free electrons, to predict the charge state distribution.Because of the free electron screening, the model is relatively simple to implement and is computationally cheap, but still includes the effect of orbital relaxation.The variational model uses a constrained free energy minimization to determine the charge state distribution.The key constraint is the average ionization, here provided by the Tartarus average atom model [28].It is argued that the use of free electron screening is justified (as opposed to solving the Schrödinger equation for the continuum states) since in either case there are more important physical approximations being made.We compare predictions of the charge state distribution to other models, finding good agreement with other average atom based models that do not use this free electron approximation.Disagreements with the ChemEOS model [27,29,30] are discussed.

A. Variational Model for Charge State Distribution
The physical model is of ions and an electron gas in spheres of equal volume V, given by the average volume per ion, and an overall average ionization Z.We require each sphere to be charge neutral, which should be a reasonable approximation for LTE plasmas.If an ion corresponding to configuration x has Q x bound electrons, the remaining electrons are considered to be free, with a density n 0 x , such that Z x = Q x + Zx , where Z x is the nuclear charge of the ion, and Zx = V n 0 x is the number of ionized electrons for the ion.The free energy F per ion of a plasma composed of ions with a neutralizing electron gas is where the sum is over all electronic configurations of the ions in the plasma, and we are defining a configuration to mean a set of integer occupation numbers for orbitals defined by their principal, orbital angular, magnetic, and spin quantum numbers.W x is the probability that configuration x exists, T is the temperature of the plasma, F eg is the free energy of the electron gas 1 , and U x is the configuration average internal energy of configuration x.We have used Hartree atomic units in which h = m e = k B = 1.The entropy term in equation ( 1), W x ln W x , is the entropy of mixing due to configuration x, and due to the bound electrons.The occupations of electrons in bound orbitals, for a given configuration x, are an input to the model.The complete list of the configurations to be included in the sum over x is not known a priori, but is determined by an exhaustive search which we will discuss in the next section.
While the bound orbital occupancy is fixed by the choice of configuration, in this model the free states are occupied according to Fermi-Dirac statistics.This breaks the consistency between free and bound electrons so important for thermodynamic consistency of equation of state (EOS) models [31], but it is of less importance for photon spectra due to the (generally) large number of configurations [25].In average atom based EOS models a large amount of effort has been devoted to accurate treatment of the free orbitals [28,[31][32][33].However, for spectral resolution, an important physical effect completely missing in these single center models for dense plasmas is multiple scattering, which destroys the atomic eigenstate character of loosely bound orbitals [34].In light of this, it seems like unnecessary effort to apply the same level of rigor and numerical accuracy to calculating the free orbitals when solving for the ion's electronic structure.Even if one were to solve for the continuum orbitals in the same potential as the bound orbitals, consistency with bound orbitals would still be broken and important physical effects would still be missing.
The main physical effect that we wish to capture is the influence of plasma screening on the ion's eigenvectors and eigenvalues.This can be included by approximating the free electrons as a truly free electron gas, i.e. the free electron density is homogeneous, given by where F 1 2 is the Fermi integral (definition in reference [28]), µ is the electron chemical potential, and Vx is a local variation in potential taking a value such that the correct Zx of ion x is recovered.
To ensure that the correct average ionization Z for the plasma is recovered, a constrained minimization of the 1 For which we take the average value, ignoring the variation from ion to ion.
free energy is used.Let Ω be the constrained free energy where B and C are Lagrange multipliers.The first constraint ensures that the probabilities sum to 1, the second ensures that the average ionization of the plasma is Z.
Minimizing with respect to W x and enforcing the first constraint gives where is the partition function, and the sum is now over configurations with different energies U x , with a configuration now being defined by a set of integer occupation numbers together with their principal and orbital angular momentum quantum numbers, and g x is the degeneracy of such a configuration.The second constraint is enforced by varying C until is satisfied, which is done numerically.
To close this model, Z must be provided.As there is no unique definition of this quantity, one has some freedom.We have found that using an average atom model with a 'chemical' definition works well.In this definition the number of bound electrons is the number of electrons in orbitals that asymptotically decay, and so are localized around the ion, consistent with the ion concept above.We use the Tartarus average atom model of reference [28].In the notation of that reference, the chemical definition of ionization is " Z".Note that using the alternative definition of average ionization given in reference [28] ("Z * "), is inconsistent with the present model's definition of an "ion", and in practice using Z * in place of Z sometimes causes equation (6) to have no solution.

B. Energy of an ion
The configuration average energy of an ion, U x , corresponding to an integer occupation number configuration x, can be approximated in the local density approximation (LDA) by the expression where U el x is the electrostatic energy, ∆U k x is the kinetic energy of the electrons without a free electron contribution and ∆U xc x is the exchange and correlation internal energy, again with the free electron contribution removed.U el x is given by where the electrostatic potential is and the ion electron density is where P i,x (r) is an eigenfunction of the radial Schrödinger equation, and the sum is over all bound eigenstates of the configuration, with integer occupation q i x , with 0 ≤ q i x ≤ 2(2l i + 1), and l i is the orbital angular momentum quantum number of the eigenstate.
The kinetic energy of the electrons in the configuration is given by where i x is the eigenenergy of eigenstate i in configuration x, and This normalization-like integral will be equal to 1 for deeply bound orbitals, but will be less than 1 for weakly bound orbitals whose eigenfunctions have not decayed to zero by the ion sphere radius R.
The exchange and correlation energy is given by and we have used the Perdew-Zunger functional [35] for our calculations.This results in an effective interaction potential with V xc the corresponding exchange and correlation potential.Note that the eigenfunctions are normalized over all space and the potential V ef f x (r) is formally assumed to be constant outside the spheres for this purpose.
The configuration energies are improved by replacing the LDA energy for the bound orbitals with the Hartree-Fock interaction energy [36].This provides a correction that improves configuration energies [36].In this approximation the energy correction ∆U HF x is given by where with n b x (r) = n x (r) − n 0 x , and ( Z = 6), at 100 eV and 2.7 g/cm 3 .Also shown is the 2s eigenfunction for a neutral, isolated aluminum atom.
V i,j is the shell-shell interaction energy [36].We also include the (generally small) correlation correction given in reference [36].It is worth pointing out that this Hartree-Fock correction is not self-consistent as it uses the LDA eigenfunctions.

C. Discussion of some model approximations
In comparison to other models of charge state distribution there are various differences.A central difference from references [22,24,25], which present models that also calculate density and temperature dependent eigenfunctions and eigenvalues, is that in those models the ion sphere sizes are varied from configuration to configuration so that the free electron density at the edge of the spheres is the same.Here we allow those free electron densities to vary from ion to ion, while keeping the ion sphere volume and the chemical potential the same.It does not seem probable that the amount of 'space' an ion takes up will be strongly correlated with its ionization stage (at a given plasma density), as typical ionization/recombination times are much faster than ion motion time scales.The concept of an ion sphere is itself a rather crude approximation in dense plasmas, with loosely bound valence states being strongly affected by neighboring ions [34].Hence, such models' descriptions of weakly bound states can only be approximately correct in any case.Therefore, our assumption of one ion sphere size is probably equally justified, while being simpler computationally, as there is no search for the unknown sphere sizes.
The ChemEOS model [27,29,30] starts from a different perspective.It uses a calculation of eigenstates and eigenvalues of isolated ions.The infinite bound spectrum of isolated ions leads to a divergent internal partition function and thus free energy.In ChemEOS, this is resolved with an ionization model that combines a hardsphere free-volume model with a plasma microfield model for the destruction of bound orbitals.In the present model the partition function converges because there is a finite list of configurations that exist.This is due to the electronic structure being calculated self-consistently with the plasma effects, in contrast to the ChemEOS model.
In some models average atom eigenfunctions are used as approximations to the ion eigenfunctions [25,37].In figure 1 the 2s eigenfunction from the average atom model Tartarus, for aluminum at 100 eV and 2.7 g/cm 3 , is compared to that from two configurations, 1s 2 2s 2 and 1s 2 2s 2 2p 3 , at the same temperature and density.The average atom eigenfunction is quite close to the configuration wave functions, supporting the use of these wave functions as approximations for the configuration eigenfunctions [25,37].A disadvantage of such an approximation is that, because the average atom has a finite number of bound orbitals, if an excited state configuration is needed which requires an eigenfunction that does not exist in the average atom model, one does not have an approximate eigenfunction to use.By performing selfconsistent field calculations for each configuration, the present model is able to recover excited-state orbitals which are pressure ionized out of the average atom, but which can appear in an integer-occupation ion due to orbital relaxation.

FIG. 3: (Color online)
Energies of all singly or doubly excited ions that exist in the present model, for ion stages III through XIII, for a solid density aluminum plasma at 100 eV.The three distinct bands correspond to ions with the no excitations of the 1s orbital (1s 2 . ..), a single excitation (1s 1 . ..), or empty 1s orbital (1s 0 . ..).
Also shown in figure 1 is the 2s eigenfunction from a calculation of a neutral, ground state, isolated aluminum atom.This is different from the other eigenfunctions due to density and temperature effects.While the eigenstate is well confined to within the ion sphere at 2.7 g/cm 3  (where the radius of the sphere is ∼ 3 a B ), the screening due to other electrons changes from that due only to other bound electrons, to screening from free electrons as well as bound electrons.
The present configuration model includes the effect known as orbital relaxation, where the eigenvalue of a particular nl orbital varies from configuration to configuration.This is demonstrated in figure 2, where the eigenvalues for the 2s and 2p orbitals are shown for the ground state configuration, for a range of ionization stages of aluminum.Going to higher ionization, the orbitals become more tightly bound.This is due to free electrons being less effective at screening the nucleus than bound electrons, as bound electrons are, on average, closer to the nucleus.Note that the eigenvalue is shown for the 2p orbital of the 1s 2 2s 2 configuration despite the fact that it is not occupied, and it does not directly affect the calculation of the energy of the ion.
In figure 3 we show the energies of all configurations that exist for ion stages III through XIII, according to the present model, for a solid density aluminum plasma at 100 eV.To determine whether or not an ion configuration exists, the model is run with the desired set of integer occupation numbers, and if a self-consistent field solution can be found for that configuration, then the configuration exists.As one does not know the final list of configurations that exist a priori, many more configu-  [27,29,30], and red triangles are from reference [24], where we have used their model labeled 'DFT'.Also shown are open black circles which come from the present variational model, but using the ChemEOS value for the average ionization, and a Saha-Boltzmann (SB) calculation using the present ion configurations and energies.
rations than those that are found to exist must be tested.In the figure, three distinct bands are predicted, corresponding to a fully occupied, singly occupied, or empty 1s orbital.In the Tartarus model the eigenenergy of the 1s orbital is −61.6 E H , roughly equal to the gap between the observed bands.Note that a fully ionized atom would have zero energy on this scale.In the figure, ion stages with more ionized electrons have energies closer to this limit.

III. RESULTS
In figure 4 we shown the charge state distribution for an aluminum plasma at solid density and a temperature of 100 eV.Shown is the result from the present model compared to the recent model of Faussurier and Blancard [24] (FB) and to the ChemEOS model [27,29,30].The FB model uses varying ion-sphere sizes to maintain the same free electron density at the edge of the spheres, has built in plasma screening, and uses a Saha-Boltzmann formalism to determine the charge state distribution and average ionization.The agreement between the present model and the FB model is quite good.The FB model corresponds to a lower average ionization ( Z = 7.60) than the present model ( Z = 7.90).This may be due to the Saha-Boltzmann model, which is known to underestimate Z at high material densities because it is based on classical statistical mechanics.We have also implemented the Saha-Boltzmann model using the ion configurations and energies of the present model.This result is also shown in figure 4, as open blue circles ( Z = 7.45).While it does not agree exactly with the FB model (due to the different ion models used), it is shifted to lower ionizations than the full variational model, indicating that the Saha-Boltzmann model is, at least partially, the cause of the lower average ionization.Also shown in figure 4 is the result from the ChemEOS model.It gives a significantly lower average ionization ( Z = 6.95) than the present variational model ( Z = 7.90), but the shape of the distribution is similar.The present average ionization comes from the Tartarus average atom model, whereas Z for ChemEOS comes from the above mentioned hard-sphere and microfield model.We have also used the ChemEOS Z to constrain the present variational model, and the result is also shown in figure 4, as open black circles.There is very agreement with ChemEOS, indicating that the different methods of calculating the average ionization are indeed the cause of the differences, while the shape of the charge distributions is similar in both models for this case.
In figure 5 charge state distributions for aluminum plasmas at temperatures of 10 and 100 eV, and three orders of magnitude in density, ranging from solid density to 1/100 th of solid, are shown.Shown are the present variational model using ion energies with and without the Hartree-Fock correction.This correction makes little difference to any of the cases shown because it amounts to a roughly constant shift in configuration energy.It should be important, however, when calculating spectra [37].Also shown in the figure are the results from the ChemEOS model and our own Saha-Boltzmann calculation.Agreement between the variational model and ChemEOS is generally very good, for all cases.The exception is the 100 eV, solid density case discussed previously with figure 4. The Saha-Boltzmann model also compares well to the variational model at 100 eV, but is less reliable at the lower temperatures, and higher densities, as expected.
In figure 6 we show charge state distributions for iron plasmas at a temperature of 40 eV.We compare the present model to the model of reference [25], which is based on an average atom model with a variational formula for the charge state distributions.That model uses the average atom eigenfunctions and the free state density for all ion configurations, with a varying ion-sphere size.The level of agreement is quite good, and only small differences are observed.Such differences would translate to small but visible effects in spectra.
Also shown in figure 6 are the results of the ChemEOS model, which consistently predicts a lower average ionization than the present model, especially as density increases.This is due to the very different way that Z is calculated.In the present model plasma effects are included self-consistently in the calculation of Z, in contrast to the ChemEOS model.Thus, we may expect the present model to be more realistic for this quantity where plasma effects are more significant.However, as Z is not uniquely definable, we must be cautious before drawing strong conclusions.
Lastly, shown in figure 6 is the result from our own calculation using the Saha-Boltzmann formalism.At low density agreement between it and the full variational model is good, but gets worse as density increases.This is to be expected, as the Saha-Boltzmann model becomes steadily less reliable as density increases, for fixed temperature.

IV. CONCLUSIONS
We have presented a new variational model for calculating the charge state distribution in LTE plasmas.The model is based on a free energy minimization, constrained to give the desired free electron density.The inputs to the model are the energies of the configurations that exist, and the average ionization Z.For Z we use the average atom model Tartarus.For the configuration energies we developed a free-electron screened ion model, which is used to determine whether a configuration exists in the plasma, and in which each configuration has a finite list of bound orbitals (occupied and unoccupied), due to free electron screening.
Comparison with another variational model, reference [25], reveals good agreement.An important difference with that model is that the present model includes orbital relaxation, which will be important for spectra calculations.Comparison with the recent model of reference [24] also results in good agreement.That model uses the Saha-Boltzmann formalism, and using our own implementation of this, we demonstrated its breakdown at high densities and lower temperatures.Finally, comparison with the ChemEOS model [27,29,30] reveals mixed agreement.Typically, agreement was better for lower densities, as expected.

2 FIG. 2 :
FIG.2:(Color online) Eigenvalues for the 2s and 2p orbitals for the ground states of ion stages VII through X, for aluminum at solid density and 100 eV.The upper green lines are for the 2p orbital, while the lower blue lines correspond to the more tightly bound 2s orbitals.The fact that the 2s or 2p orbitals have different energies in different configurations is called orbital relaxation.

FIG. 4 :
FIG.4: (Color online) Charge state distributions for aluminum at 100 eV and solid density.Filled black circles are current variational model with HF energy for bound states, filled green squares are the ChemEOS model[27,29,30], and red triangles are from reference[24], where we have used their model labeled 'DFT'.Also shown are open black circles which come from the present variational model, but using the ChemEOS value for the average ionization, and a Saha-Boltzmann (SB) calculation using the present ion configurations and energies.

FIG. 5 :
FIG. 5: (Color online) Charge state distributions for aluminum plasmas at various temperatures and densities.Filled black circles are current variational model with HF energy for bound states, red crosses are current variational model with LDA energy for bound states, open blue circles are from Saha-Boltzmann model with HF ionization energies from the present ion model, and filled green squares are the ChemEOS model.

3 FIG. 6 :
FIG. 6: (Color online) Charge state distributions for iron plasmas at various densities a temperature of 40 eV.Filled black circles are current variational model with HF energy for bound states, open blue circles are from Saha-Boltzmann model with HF ionization energies, filled green squares are the ChemEOS model, and open red triangles are from reference[25], where we have used the model labeled 'eq.(43)'.