Improved calculations of mean ionization states with an average-atom model

The mean ionization state (MIS) is a critical property in dense plasma and warm dense matter research, for example as an input to hydrodynamics simulations and Monte-Carlo simulations. Unfortunately, however, the best way to compute the MIS remains an open question. Average-atom (AA) models are widely-used in this context due to their computational efficiency, but as we show here, the canonical approach for calculating the MIS in AA models is typically insufficient. We therefore explore three alternative approaches to compute the MIS. Firstly, we modify the canonical approach to change the way electrons are partitioned into bound and free states; secondly, we develop a novel approach using the electron localization function; finally, we extend a method which uses the Kubo-Greenwood conductivity to our average-atom model. Through comparisons with higher-fidelity simulations and experimental data, we find that any of the three new methods usually out-performs the canonical approach, with the electron localization function and Kubo-Greenwood methods showing particular promise.


I. INTRODUCTION
Warm dense matter (WDM) is a phase of matter characterized by temperatures on the order of 1 − 100 eV and densities of 10 −2 − 10 4 g cm −3 [1,2]. Under these conditions, conventional divisions between solidstate and plasma physics are bridged and a variety of interesting phenomena emerge, including for example nonequilibrium effects [3], phase transitions [4,5], and partially ionized matter. WDM is observed in various astrophysical domains, such as exoplanets [6] and brown and white dwarfs [7,8]; furthermore, during inertial confinement fusion (ICF), materials are exposed to WDM conditions [9,10].
The mean ionization state (MIS), or equivalently the free electron density, is of particular importance in WDM. It is directly related to physical properties such as electrical conductivity, opacity, collision rates and acoustic velocities [11,12]. Furthermore, the MIS is an input parameter for various simulations including hydrodynamics [13] and Monte-Carlo simulations [14], finite-temperature pseudo-potentials for densityfunctional theory calculations [15,16], and in computing adiabats used in ICF modelling [10]. Additionally, accurate predictions of the MIS are crucial for validating and fitting models to experimental data [17,18].
In the WDM regime, it is often difficult to distinguish between 'bound' and 'free' electrons, meaning the MIS is hard to define. The ramifications of this ambiguity extend beyond direct computation of the MIS: they are relevant to recent debates regarding the ionization potential depression (IPD) effect [19][20][21][22], and further raise questions regarding the application of the Chihara decomposition [23][24][25]. These difficulties are further compounded by the variety of methods used in the modelling of WDM, running all the way from analytical models such as Stewart-Pyatt [26] and Ecker-Kroll [27] to ab initio density-functional theory (DFT) [28][29][30][31][32] and pathintegral Monte-Carlo [33,34] simulations. It is therefore of great interest to develop an approach for calculating the MIS that is consistent between different models and experimental results.
Average-atom (AA) models are a popular and successful tool in modelling the WDM regime, since they incorporate in a natural way quantum effects (typically using DFT) at a manageable computational cost [35][36][37]. There is a wide range of AA models [38], but they share in common the concept of an atom immersed in a plasma. Typically, the MIS is defined as the number of electronic states with energy above a certain threshold, where g( ) denotes the density-of-states, f FD ( ) the Fermi-Dirac (FD) distribution, and 0 the chosen energy threshold. In AA models, the threshold energy is typically chosen to be the value of the mean-field potential at the boundary of the Voronoi cell R VS (the atomic radius), 0 = v s (R VS ). Other choices for 0 , for example equating it to the chemical potential, could also be considered.
As seen in a previous work [38], the definition (1) is somewhat limited, showing large discrepancies for different choices of boundary condition and sharp discontinuities whenZ is plotted as a function of temperature or density. Furthermore, bound and free states in AA models are typically treated differently (although not always, for example Refs. [39,40]): the definition ofZ is thus both an output of and input to the model, which means any errors may self-multiply.
DFT-based molecular dynamics (DFT-MD) simulations can also be used to compute the MIS using definition (1), with the threshold energy typically assumed to start at the conduction band lower edge, 0 = c . In DFT-MD simulations, all (non-core) orbitals are treated on the same footing, which is an advantage (among others) relative to AA models. However, there are still (at least) two limitations using this definition, which are common to both DFT-MD and AA models. The first is the ambiguity about how to define the threshold energy. The second is the assumption that states can be categorized as completely bound or completely free based on their energy alone: DFT-MD results with this method have shown counter-intuitive behaviour [41] and divergence from experimental measurements [17].
Consequently, novel ways of computing the MIS have recently gained traction. For example, Bethkenhagen et al. proposed using the Kubo-Greenwood (KG) conductivity formula to measure the MIS [11]. This approach was applied to Carbon under high temperatures and gigabar pressures (and later to the metallization of helium [41]), and the resulting MIS values showed disagreement with various other methods. Interestingly, excellent agreement was seen between pressures computed with an AA model and DFT-MD under these conditions [42]. However, the MIS computed with the same AA model -using yet another definition forZ -had a systematic error relative to the DFT-MD KG result, which suggests that a more pertinent definition of theZ in an AA model might give better agreement.
In this paper, we explore three methods for computing the MIS in an AA model, and compare results with DFT-MD simulations [11] and experimental data [43][44][45]. Firstly, we apply the canonical definition (1), which (as expected) gives inconsistent results, particularly for high densities. Secondly, we modify the canonical approach such that the orbitals are no longer categorized as bound or free based on their energy. Instead, they are partitioned depending on their shell (1s, 2s, etc), an approach that was used for the (non-average-atom) XCRYSTAL model in Ref. 46. Thirdly, we introduce a novel approach which uses the electron localization function (ELF) to determine the MIS. The ELF is well-known in quantum chemistry and materials science [47,48], but has not until now been applied to study ionization in WDM. We shall see that this method yields more consistent and accurate results compared to the canonical approach. Finally, we adapt the KG method of Ref. 11 to our AA model. This approach shows excellent agreement both with DFT-MD simulations and the experimental results, but is so far limited to only one boundary condition in the AA model. Nevertheless, the ELF and KG results demonstrate that computationally efficient AA models can accurately and reliably predict the MIS across a wide range of conditions.

A. Average-atom model
The AA model we use is a generalization of the model derived in Ref. 38. We explain here the main features of this model and the differences from the one presented in Ref. 38; however, we direct readers to that paper for a detailed derivation and discussion of this AA model. In our AA model, we solve the Kohn-Sham DFT (KS-DFT) equations for a single atom consisting of a nucleus with charge Z and a fixed number of electrons N e (with N e = Z for all the systems we consider). Explicit interactions between this atom and its neighbours are ignored, and instead these interactions are implicitly accounted for via the boundary conditions imposed on the orbitals at the sphere's edge (Voronoi sphere radius, R VS ).
The spherically symmetric KS equations to be solved are given by [49] with r > (x) = max(r, x). The three terms in the potential are respectively the electron-nuclear attraction, the classical Hartree repulsion, and the exchange-correlation (xc) potential, which is equal to the functional derivative of the xc free energy. As ever, due to the dependence of the KS potential on the density n(r), the KS equations must be solved iteratively until self-consistency is reached. The density n(r) is constructed from the orbitals as where f nl ( nl , µ, T ) is the Fermi-Dirac (FD) distribution, given by we have not made any changes to avoid too much speculation.
The chemical potential µ is determined by fixing the electron number N e = 4π RVS 0 drr 2 n(r) to be equal to a pre-determined value (in this paper, N e = Z in all cases).
We impose boundary conditions on the KS orbitals X nl (r) which are intended to implicitly account for interatomic interactions. In our earlier paper [38], we argued that a physically intuitive condition was to impose smoothness of the density at the edge of the Voronoi sphere (VS), dn(r) dr Mathematically there is no unique way to enforce the above condition, but two simple choices are which we refer to respectively as "Dirichlet" and "Neumann" conditions. From a theoretical standpoint within the AA model, there is no way to unambiguously differentiate between these boundary conditions. We now note a key improvement we have made to our AA model compared to Ref. [38]. In that paper, we only solved the KS equations (2) for the 'bound' electrons, defined as those with energies below the threshold energy 0 . For the remaining 'unbound' electrons, we used the ideal approximation, which amounts to assuming a constant density for the bound electrons, n ub (r) =n. In this work, we make no distinction between 'bound' and 'unbound' orbitals during the SCF procedure: in other words, we solve the same equations (2) for all orbitals, regardless of their energy. As already mentioned, this removes the issue of the MIS being both an input to the model (via the ionization threshold 0 ) and an output of it. Moreover, as we shall soon see, the expressions for the KG conductivity and the ELF are explicitly orbitaldependent; it therefore does not make sense to calculate these properties when part of the density is constructed in an orbital-free manner as we had done in the past.
Furthermore, to extend our comparisons beyond the model described above, we have implemented the AA model proposed by Massacrier et al. [40]. In this model, the KS equations are solved for both the Dirichlet and Neumann boundary conditions, yielding energies ± nl which define the upper (Dirichlet) and lower (Neumann) limits of a band-structure. Within these limits, every energy value is permitted and the wave-function corresponding to that energy is determined. The KS equations thus become The Fermi-Dirac occupations are multiplied by the Hubbard density-of-states (DOS) function g nl ( ), defined as [50] g nl ( ) = 8 which means the density in this band-structure model is given by (12) In practice, the energy bands are discretized, and the above integral becomes a summation over energies within each band which we now denote by index k. Following some algebraic manipulation, the density can be written as where N k is the number of points used in the discretization of each energy band. The above expression closely resembles the expression for the density in plane-wave DFT codes, since it has a summation over k-points and some weighting w k (with k w k = 1), very much like the k-point mesh for reciprocal space. It is also clear to see that when the concept of bands in the AA model is not employed (i.e. when we use either the Dirichlet or Neumann conditions only), that N k = 1, w k = 1 and the above expression reduces to the ordinary expression for the density (4). The above simplification (13,14) was not shown in Ref. [40], and we therefore provide a derivation in Appendix A.

B. Counting method
As discussed in the introduction, and we shall later see in the results, the canonical definition of the MIS in AA models (1) is often erroneous, which is why we shall explore alternative approaches. In the following three sub-sections, we discuss the application of new methods -first, the counting method, which is a modification to the threshold approach, secondly, the electron localization function (ELF), and lastly, the Kubo-Greenwood conductivity -to calculating the MIS.
As seen in Eq. (1), the canonical approach to computing the MIS essentially defines electrons as bound or free depending on whether their energy exceeds some threshold 0 , typically defined as the value of the KS potential at the edge of the atomic sphere. Intuitively, this does make sense, if one imagines electrons being bound so long as their energies are below the maximum value of the KS potential, and otherwise free. However, as we shall see in the Results section, this leads to unphysical discontinuities in the MIS when the energy of an orbital crosses the threshold value, and is very sensitive to the choice of boundary conditions. Rather than making the bound-free partition dependent on some energy value, we instead propose to partition the electrons based on their shells. This method was used to compute the MIS in Ref. 46 for the XCRYSTAL model (specifically Eq. (22) and the surrounding discussion), but has not been applied (as far as we know) to average-atom models. In fact, the argument they use for this approach -"our flat potential V 0 does not share the same physical interpretation as the flat potential used in Ref. [29], as delocalized states can be found below V 0 in XCRYSTAL" -is applicable to average-atom models such as ours, in which there are no constraints on the KS potential.
The method is perhaps best illustrated with an example. Consider Aluminium at its ambient density, ρ m = 2.7 g cm −3 . It is well-known that, at room temperature, the 1s, 2s and 2p orbitals are core states, and the remaining orbitals represent free electron density. This can also be seen by inspection of the density-of-states, using for example the average-atom band-structure model.
As the temperature is increased, the character of these core states actually does not change much, as can be seen in Fig. 1. Therefore we can essentially consider these states to represent bound electron density, regardless of the temperature. Of course, as the temperature increases, the occupation of these core states will decrease as higher-energy states are occupied, causing the MIS to increase. We shall henceforth refer to this approach as the 'counting' method, and it can, in theory, be generalized to any material at a given density. The expression for the MIS in this counting method is where b denotes that subset of orbitals considered to be bound. Although we have included a k-dependence in the above sum, we have done so for generality; ideally, the bound states should be clearly identifiable as core states, in other words, their energies should not form a band and the k-index should be redundant. Clearly, the approach described above works best if orbitals can be clearly identified as being of bound or free character, as is typical for metals at their ambient density (for example). However, when this is not the case -in particular when a range of densities is spanned for a given material -the above method is likely to break down. As material density changes, the orbital character also changes significantly, bands emerge and disappear, and so on. In such scenarios, one would expect this counting method to fail. In the Results section, we shall see that this expectation is borne out.

C. Electron localization function
In this subsection, we describe the method we have developed to compute the MIS with the electron localization function (ELF). The ELF has a long history in quantum chemistry [47,48,51] as a tool for understanding atomic structure and chemical bonding. It was originally conceptualized by Becke and Edgecombe [51], who supposed that the conditional probability density -i.e., the probability of finding an electron at position r 1 given another electron with the same spin at position r 2could be used as a basis to measure electron localization. It was later generalized by Savin [47] such that any spinindependent electron density could be considered.
In KS-DFT, the expression for the (total density) electron localization function (ELF) is given by where D(r) and D 0 (r) are the electron pair density curvature (EPDC) functions for the system and for the uniform electron gas (UEG) respectively. These are given by where τ (r) is the kinetic energy density. There are in fact multiple ways to define τ (r) [52][53][54], which of course all yield the total kinetic energy when integrated over all space. The definition most commonly adopted in the ELF is [48] τ The motivation for the definition of the ELF (16) is to define electron localization in a quantitative manner, by using the EPDC of the UEG, a perfectly delocalized electron density, as a reference. The ELF is bounded in the range 0 ≤ ELF ≤ 1: a value of 1 indicates strongly localized electron density and a value of 1/2 indicates equivalence with the (delocalized) UEG.
One of the principal uses of the ELF is to calculate the number of electrons in particular shells. In the atomic picture, the spatial boundary of the shells is equated to minima in the ELF. Then, the density is integrated between minima to give the number of electrons in that shell. A visual example of this procedure is shown in Fig. 2.
We propose to use the ELF as a measure of the MIS by computing the number of electrons per shell, and assuming that any electron density beyond a particular shell is free. This presents a similar issue to the counting method described in the prior sub-section; however, as we shall see, the ELF method is advantageous when a scan over densities is performed. Nevertheless, this does introduce some ambiguity and means this approach cannot be considered a "black-box" method. We see that the 1s state is unaffected by the temperature, and the 2s and 2p states are moderately affected, but not so much to change their bound-state character. The calculation was done with the Dirichlet boundary condition, but due to the core-nature of the orbitals, the boundary condition is of minimal impact. In the application of the ELF to our AA model at moderate-to-high temperatures, we have observed that, using the normal definition of the kinetic energy density (19), the ELF's minima are often not identifiable. However, we have found that an approximate expression for the kinetic energy density τ (r), based on a second-order gradient expansion [55], yields more clearly identifiable minima in the ELF than the normal orbital-dependent expression. This approximation for τ (r) is given by which leads to the following expression for D(r), In spherical co-ordinates, this becomes In Fig. 3, we compare the ELF computed using the usual definition of the kinetic energy density (19) with the approximate form (20). We compare three temperatures: 0.01 eV, 10 eV and 100 eV, and consider both Dirichlet and Neumann boundary conditions. For 0.01 and 10 eV, we see that the shape of the ELF is in general different for the different forms of the kinetic energy density; however, the positions of the first two minima, which correspond to the boundaries of the n = 1 and n = 2 electron shells, are almost identical. On the other hand, at 100 eV, the n = 2 minimum is no longer identifiable when the orbital-based definition (19) for the kinetic energy density is used; this is in contrast to when the approximate density-based definition (20) is used, in which case the n = 2 minimum is clearly visible.
In general, we have observed the tendency for the orbital-based expression (19) to break down as temperature increases for a range of materials and densities. Consequently, we prefer to use the approximate definition (20), which does not display the same tendency, for all calculations of the MIS. As is observed in Fig. 3, particularly in the right-hand panel of this figure, this can also produce additional and unexpected minima in the ELF. It is unclear whether these minima are really physically connected to electron shells, or are simply artifacts from the average-atom model and boundary conditions. Regardless, since we assume all electron density beyond a certain shell (n ≥ 3 in this example) is free, a correct physical interpretation of these additional minima is not strictly required in this approach.

D. Kubo-Greenwood conductivity
In this sub-section, we describe the application of the Kubo-Greenwood conductivity to compute the MIS within our AA model. The Kubo-Greenwood (KG) conductivity formula for a finite system is given by [56,57] σ S1, where σ(ω) is the dynamical conductivity for two subsets S 1 and S 2 of the orbitals, V is the volume of the system under consideration, φ i are the KS orbitals and i and f i are their energies and FD occupations. For the total conductivity, S 1 and S 2 represent the complete set of orbitals.
As described in Ref. [11], Eq. (23) can be used as a proxy for the mean ionization state in combination with the Thomas-Reiche-Kuhn (TRK) sum rule [58][59][60]. This rule establishes a relationship between the KG conductivity and a certain number of electrons. For example, if we take S 1 and S 2 to both be the complete set of orbitals, then we should recover the total electron number, where σ t,t denotes the conductivity from the total, or complete, set of orbitals. We note here that the complete set of orbitals means, in theory, an infinite set of KS orbitals (i.e. not just those with non-zero occupation numbers). In practice, a sufficient number of orbitals is chosen such that the resulting electron number is equal (within reasonable tolerance) to the expected electron number. This provides a useful check of the implementation and convergence of the KG method.
To calculate the MIS, we usē where σ c,c means both orbital subsets are given by the conducting orbitals.
In the spherically symmetric AA model, the KG conductivity is given by which leads to the following expression for Z S1,S2 , Z S1,S2 = 4 nl∈S1 n l ∈S2 m∈{S1,S2} f nl − f n l n l − nl |∇ z nn ll m | 2 δ(l ± 1 − l )Θ( n l − nl ) . (27) In the above equations, ∇ z nn ll m is the z-component of the momentum integral matrix product, and Θ( n l − nl ) is the Heaviside step function. The derivation of the above expressions, and the expression for ∇ z nn ll m in terms of the radial KS orbitals and spherical harmonic functions, can be found in Appendix B.
In a conventional AA model, unlike in plane-wave DFT calculations, there is no concept of a band-structure, which is problematic for determining which subset of orbitals belongs to the conducting and valence bands. We could, for example, use a threshold energy as the dividing line between conduction and valence electrons. However, since we have the band-structure AA model at our disposal, we can use that to guide which orbitals belong in the conduction and valence bands. This is just done manually (e.g. by inspecting the DOS, see Fig. 4). Even when the conductivity is evaluated with the Dirichlet or Neumann boundary conditions, we use the band-structure model to determine the valence and conduction bands.
In Fig. 4, we plot the DOS given by the AA bandstructure model for Carbon at 100 eV and various densities. In this case, there is a clear valence band (to the left of the dotted lines) and conduction band (to the right). Through inspection of the energies, the valence band can be associated with the orbitals in the 1s band. Therefore, when evaluating the KG conductivity with the Dirichlet or Neumann condition, the 1s orbital is assigned to the valence band and all others to the conduction band. The same strategy is used in applications of the KG method in this paper.

III. RESULTS
All calculations have been performed using the opensource average atom code atoMEC [61,62]. In Ref. 61, we describe the structure of the code, together with general algorithmic and numerical details. Numerical details specific to this paper are given in the Supplemental Material [63]. We note that the following libraries are used extensively by atoMEC: NumPy [64], SciPy [65], LIBXC [66], mendeleev [67], and joblib [68].
In the following, we shall compare the four methods described for computing the MIS -the canonical or "threshold" approach (1), the counting method (15), the ELF method, and the KG approach (25) -against a higher fidelity DFT-MD benchmark and experimental data. For the threshold, counting, and ELF results, we compare the Dirichlet and Neumann boundary conditions and the band-structure model [40]. For the KG results, we use the Dirichlet boundary condition only. This is because the sum rule check for the total conductivity is observed very accurately (within 1%) across all conditions for the Dirichlet boundary condition, but not for the others. We use throughout the (spin-unpolarized) local density approximation (LDA) for the xc-functional [69]. First, in Fig. 5, we compare our results with the DFT-MD simulations for Carbon from Ref. 11. In Fig. 5. (a), we plot the MISZ using the canonical threshold method for each of the boundary conditions. We see that this method has severe limitations, especially at the highest densities, when the three boundary conditions yield completely different results. Furthermore, in the densityrange in which the DFT-MD simulations were performed, none of the AA results are remotely close to the reference result.
In Fig. 5. (b), we plot the MIS using the counting method. In order to obtain these results, we took the electron density in the 1s orbital to be bound, and everything outside it to be free. In this case, we observe that the results are quite consistent between the different boundary conditions. However, they all tend to the wrong limit as the density increases. This is a consequence of the assumption that some orbitals -in this case, the 1s orbital -are bound states for the whole range of densities. Clearly, from both an intuitive and results-based perspective, this is not the case. Consider, for example, Fig. 4: as the material density increases, the part of the DOS that comes from the 1s orbital (to the left of the left dotted line) transforms in nature from a delta-like function (signifying bound electron density) to a wide energy band (signifying free electron density).
In Fig. 5. (c), we plot the MIS obtained via the ELF method against the DFT-MD benchmark. In order to obtain these results, we took the electron density in the n = 1 sub-shell to be bound, and everything outside it to be free. We see that this approach yields a more realistic picture for the MIS, as the results from the three boundary conditions are at least consistent and capture the correct qualitative behaviour; however, they all systematically under-estimate the MIS relative to the DFT-MD result. Nevertheless, it is interesting that the ELF method, unlike the counting method, does go towards the correct limit as density increases. This suggests the ELF has some ability to inherently distinguish between components of electron density with different character. Finally, in Fig. 5. (d), we compare results from the four methods (including KG) using the Dirichlet boundary condition with the DFT-MD simulation. Here we observe very strong agreement between our AA model and the DFT-MD benchmark for the KG result, until the highest densities at which point the KG result seems to have the wrong asymptotic behaviour. In this region, the ELF method actually appears to show better agreement with the reference result. Of course, the other limitation of the KG method is that it currently only works for the Dirichlet boundary condition, and it is possible that we would see inconsistencies between the boundary conditions, were a comparison possible.
In Ref. [11], it was postulated that the AA result deviates from the DFT-MD result because the AA model does not account for the many-body interactions. Based on Fig. 5, there is encouraging evidence that if the same theory is used to calculate the MIS for the AA and DFT-MD simulations, then the agreement is much better. Physically speaking, it is perhaps not unexpected that the KG result differs from the ELF and threshold approaches. After all, the KG conductivity is a frequency or timedependent property, derived by considering the linear response of a system to a perturbation; on the other hand, the ELF and energy threshold are static properties. On that basis, we should not necessarily presume consistency between the different methods.
Next, we perform a similar set of comparisons for Beryllium in Fig. 6, this time with fixed density equal to its ambient density (ρ m = 1.85 g cm −3 ) from temperatures between 1 − 80 eV. This time, the benchmark results (shown as the 3 scattered points with error bars) are taken from an experiment, in which the free electron density n f was determined using X-ray scattering [43]. The free electron density is directly related to the MIS, where V = (4/3)πR 3 VS is the volume of the atom. Like in the prior Carbon example, we have assumed under these conditions that the electron density in the n = 1 (i.e. the 1s orbital) shell is bound, and everything outside it is free. Again, the threshold results are shown in the topleft panel (a), the counting results in the top-right (b), and the ELF results in the bottom-left (c). This time, we see better agreement between the threshold results for the different boundary conditions, although the Neumann result is significantly different from the others at low temperatures. The counting and ELF results are somewhat similar, but resolve this inconsistency at low temperatures. Whilst all three techniques seem to capture roughly the right shape of the curve and agree quite well with the highest-temperature experimental measurement, they under-estimate the MIS for the lowertemperature results.
In Fig. 6. (d), we compare all three approaches for computing the MIS (threshold, ELF and KG) with just the Dirichlet boundary condition against the experimental data. Intriguingly, the KG results are in very close agreement with the lower temperature experimental results, although slightly over-predict the free electron density at the highest temperature. The KG result for the lowest temperature (τ ≈ 2 eV) is particularly interesting, because it is the only method which correctly predicts the experimentally measured value of ≈ 2.8×10 23 cm −3 : this is higher than the value which we might naively expect if we take ambient density Beryllium to have two free electrons per atom, which corresponds to n f = 2.45 × 10 23 cm −3 .
The final comparisons we make are with a pair of experiments, both involving Aluminium at its ambient density (2.7 g cm −3 ). In the first experiment [44], the free electron density n f and electron temperature were measured. We compare our AA results using the different methods and boundary conditions with the experimental data in Fig. 7. In fact, under these conditions, the majority of the AA results actually lie within the experimental error bars, regardless of the method or boundary condition used. However, a notable exception is the prediction for n f given by the Neumann boundary condition with the threshold method (which has a sharp discontinuity at around 30 eV), further demonstrating the limitations of the threshold approach. Nevertheless, Fig. 7 indicates the AA model seems to be generally accurate under these conditions, independent of the method used to compute the MIS. It appears that the ELF method with Neumann boundary condition is in particularly good agreement with the experimental benchmark in Fig. 7. (c). Since this is the only example to show such strong agreement, we prefer not to interpret this observation, but rather attribute it to chance. In the second experiment [45], the free electron density was not itself measured, but rather the K-shell ionization energy for different charge states. We use this data indirectly in the following way to compare our methods for calculating the MIS. For a range of temperatures between 1 − 100 eV, we compute the MIS and equate it to the charge state. We then take the K-edge ionization energy as the energy required to excite the 1s orbital to the continuum (with the continuum assumed to start at c = v s (R VS ) in our model). We also follow Ref. [39] and shift the orbital energy by a constant equal to the difference between 1s − c and the experimentally measured K-shell ionization energy E 0 K at zero temperature (1559.6 eV, [70]). Therefore the ionization energy is given by where 0 c and 0 1s are the threshold energy and 1s orbital energy computed with the AA model at zero temperature. This shifting is necessary because it is wellknown that KS-DFT systematically under-predicts ionization potentials using standard xc-functionals [71,72].
This comparison is shown in Fig. 8. In Fig. 8. (a), we again see several discontinuities in the threshold re-sults (for the Neumann and Dirichlet conditions), and a systematic deviation from the experimental results for higher charge states. The reason for these discontinuities, as discussed in detail in Ref. [38], is because the KS orbital energies are temperature-dependent; if an orbital crosses the energy threshold at a certain temperature then it will change its classification from bound to free (or vice-versa) and the MIS will change instantaneously. An advantage of the band-structure model is that it is not prone to discontinuities in the MIS as a function of temperature, as can be seen in Fig. 8. This is because occupations of the non-core states in the bandstructure model are spread across a band: as the limits of the energy band change, the MIS smoothly changes. This demonstrates a significant advantage of the bandstructure model when the threshold method is used.
In Fig. 8. (b), the counting approach is seen to yield consistent agreement, both internally between the different boundary conditions, and with the experimental benchmarks. For most charge states, the results lie just inside the experimental error, with some deviation seen as the charge state goes above 6 (corresponding to temperatures 60 eV). In Fig. 8. (c), the ELF results are self-consistent between the boundary conditions but also display the same systematic deviation from the experimental data. This is likely a result of the MIS being under-estimated by these methods, as we have seen for the previous examples. However, in Fig. 8. (d), we observe that the KG results lie consistently within the experimental range. This shows further promise that the KG approach for calculating the MIS agrees very well with experimental measurements.

IV. SUMMARY AND DISCUSSIONS
In this paper, we have explored different ways of computing the mean ionization state (MIS) -an essential property in warm dense matter and dense plasmasusing a KS-DFT average-atom model. Following comparisons of the different methods with DFT-MD results and experimental data, we summarize the main findings of our paper below.
The canonical method for computing the MIS, which partitions the orbitals into bound and free states based on their energies, is generally insufficient. It often causes unphysical discontinuities, and inconsistencies between different boundary conditions. If it is to be used, it is much safer to do so with the band-structure model [40], since this avoids (at least as a function of temperature) the discontinuities.
We have explored an approach which we call the 'counting' method (which was also used for the nonaverage-atom XCRYSTAL model in Ref. [46]), where the orbitals are partitioned into bound and free states based on some pre-defined intuition. This does not suffer from the discontinuities present in the threshold method, and also yields consistent results between the boundary conditions. However, it breaks down when orbitals cannot be a priori identified as being strictly bound or free in character.
We have developed an approach which uses the electron localization function (ELF) to partition the orbitals. Like the counting method, this requires a choice by the user as to which shells should be considered bound or free; however, the shells in this case do not necessarily correspond directly to particular orbitals, and so it yields better results than the counting method when the material density is varied.
We have applied a method which uses the Kubo-Greenwood conductivity [11] to our average-atom model. This also requires a choice by the user regarding a separation of orbitals into valence and conducting bands, but the resulting MIS has a sophisticated non-linear dependence on this separation. This seems to yield the strongest agreement with DFT-MD and experimental benchmarks. However, so far we have applied it only to the Dirichlet boundary condition, since sum rules are not satisfied for the other boundary conditions. Roughly speaking, we observe two different physical situations in this paper. In one instance, Figs. 6, 7 and 8, the temperature is varied for a metallic material whose mass density is fixed to its ambient value. This case is relatively straightforward: with the exception of the canonical approach with the Dirichlet and Neumann boundary conditions, all the methods yield good agreement with the benchmark data. This is because, for metals under a wide range of temperatures, the core orbitals do not undergo much change in character so can always be treated as bound states.
The other instance, Fig. 5, in which the material density (in this paper, Carbon) is varied at fixed temperature, is far more challenging. Neither the threshold or counting method is sufficient in this case; however, both the ELF and KG methods yield promising results.
It is worth noting that the KG approach has a fundamental difference compared to the other methods, since it is based on a dynamic rather than static theory. Empirically, it seems to yield systematically higher predictions for the MIS than the other methods, and also seems closer to the experimental benchmarks. This perhaps follows from the technique used to determine the free electron density in such experiments.
Based on the previous point, it may be that the "best" method to compute the MIS depends on what is desired. If the aim is to compare or provide data for an experimental fitting, the KG approach would appear to be the best approach. However, it may be that for other purposes, such as when the MIS is used as input for hydrodynamics codes, alternative methods could be favourable. This point will benefit from further investigation in future.
As a final comment, we note that more experimental data would help identify which method is most accurate across the widest range of conditions. However, high-quality experimental measurements of the free electron density (or MIS) are not trivial to come by. The assumptions used to calculate the MIS -for example, from the ratio of the inelastic to elastic scattering in Xray scattering experiments [43] -may be more likely to break down under the "harder" case of a material whose density is varied, as described earlier. This presents a major challenge for bench-marking different approaches for calculating the MIS.
In summary, the methods and data we have presented in this paper should indicate when certain methods for computing the MIS in average-atom models work, and when they might be expected to break down. With two of the methods -the ELF and KG approaches -the results are promising for all the examples we have tested. This is of particular interest because our AA code can typically run on a laptop in the time-scale of minutes -far less computationally demanding than DFT-MD simulations. the MIS; and particularly the anonymous referee for suggesting the counting method. We are also grateful to the organizers of the "Average atom models for warm dense matter workshop" at UC Berkeley in June 2021, which motivated the idea for this paper. This work was partially supported by the Center for Advanced Systems Understanding (CASUS) which is financed by Germany's Federal Ministry of Education and Research (BMBF) and by the Saxon state government out of the State budget approved by the Saxon State Parliament. EK greatly appreciates the support of the Alexander von Humboldt Foundation.

A. DERIVATION OF w k TERMS IN BAND-STRUCTURE MODEL
The energy integral to compute the density in the band-structure model, Eq. (13), must be discretized in practise. It therefore becomes a summation over energies within each band which we now denote by index k, n(r) = 2 knl (2l + 1)δ knl g knl ( knl , ± nl ) × f knl ( knl , µ, T )|X knl (r)| 2 , g knl ( knl , ± nl ) = We now simplify the above expressions, because this simplification was not discussed in the original paper. Firstly, we note that the energy spacing in the discretization of the energy band δ knl is therefore given by where N k is the number of k points (the denominator is equal to N k − 1 because there are N k − 1 spacings for N k total points). The product δ knl ×g knl ( knl , ± nl ) therefore can be written as Next, we note that the energies in a band knl can be re-written as Substituting the above expressions into the product δ knl × g knl ( knl , ± nl ) leads to the following expression: δ knl × g knl ( knl , ± nl ) = 8 π(N k − 1) 2 k(N k − 1 − k) .
(38) It is clear the above equation is in fact independent of the quantum numbers n and l. The density n(r) thus becomes n(r) = 2 k w k nl (2l + 1)f knl ( knl , µ, T )|X knl (r)| 2 ,

B. KUBO-GREENWOOD CONDUCTIVITY IN THE AVERAGE-ATOM MODEL
In the spherically symmetric case, the KS orbitals are expanded in the form φ i (r) = φ nlm (r, θ, φ) = X nl (r)Y m l (θ, φ), and the KG conductivity (23) becomes σ S1,S2 (ω) = 2π 3V ω nlm∈S1 n l m ∈S2 (f nlm − f n l m ) | φ nlm |∇|φ n l m | 2 δ( n l m − nlm − ω) . (41) Note that, in the band-structure model, this becomes σ S1,S2 (ω) = 2π 3V ω k w k nlm∈S1 n l m ∈S2 (f knlm −f kn l m ) | φ knlm |∇|φ kn l m | 2 δ( kn l m − knlm − ω) , (42) similar to the KG conductivity in plane-wave DFT codes. For simplicity, and because we only use the KG conductivity with Dirichlet boundary condition in this paper, we shall present the equations without the k-index. Since the summation only involves orbitals with the same kvalue, it is straightforward to re-introduce this at the end of the derivation. We focus first on the integral component of the equation for σ(ω), which is given by φ nlm |∇ i |φ n l m φ n l m |∇ i |φ nlm (43) = 3 φ nlm |∇ z |φ n l m φ n l m |∇ z |φ nlm , (44) where the second equation (44) follows from (43) because the contribution from each cartesian component of the gradient is identical in spherically symmetric systems. We choose the z component because, in the traditional transformation between cartesian and spherical co-ordinates, this leads to a simpler set of equations. Let us now focus on the following term, φ n l m |∇ z |φ nlm = ∇ z nn ll mm (45) = R (d) nn ll P (2) lml m δ mm + R nn ll P (4) lml m δ mm , (46) which has been taken from Ref. 57. We do not derive the above expression, but instead direct readers to the aforementioned paper where it is derived in full.
The components of the matrix element (45) are given by where P m l (x) are the Legendre polynomials. Note there are some additional factors of 4π in the above expressions compared to Ref. 57, due to different conventions in normalization of the orbitals.
Returning to the expression for σ(ω), we now have σ S1,S2 (ω) = 2π V ω nl∈S1 n l ∈S2 m∈{S1,S2} (f nl − f n l ) |∇ z nn ll m | 2 δ( n l − nl − ω)δ(l ± 1 − l ) . (52) In the above, the double summation over m has been reduced to a single summation because of the presence of the the δ mm in ∇ z nn ll mm . Additionally, the δ(l ± 1 − l ) comes from sum rules in the evaluation of the P (2,4) integrals.