First-principles derivation and properties of density-functional average-atom models

Finite-temperature Kohn--Sham density-functional theory (KS-DFT) is a widely-used method in warm dense matter (WDM) simulations and diagnostics. Unfortunately, full KS-DFT-molecular dynamics models scale unfavourably with temperature and there remains uncertainty regarding the performance of existing approximate exchange-correlation (XC) functionals under WDM conditions. Of particular concern is the expected explicit dependence of the XC functional on temperature, which is absent from most approximations. Average-atom (AA) models, which significantly reduce the computational cost of KS-DFT calculations, have therefore become an integral part of WDM modelling. In this paper, we present a derivation of a first-principles AA model from the fully-interacting many-body Hamiltonian, carefully analysing the assumptions made and terms neglected in this reduction. We explore the impact of different choices within this model -- such as boundary conditions and XC functionals -- on common properties in WDM, for example equation-of-state data, ionization degree and the behaviour of the frontier energy levels. Furthermore, drawing upon insights from ground-state KS-DFT, we discuss the likely sources of error in KS-AA models and possible strategies for mitigating such errors.

The theoretical description of WDM is particularly challenging: on the one hand, established plasma physics methods do not sufficiently account for quantum effects and strong coupling in WDM; on the other hand, the length, time and temperature scales of WDM often render popular approaches from condensed-matter physics computationally impractical. More formally, these difficulties can be understood in terms of several dimensionless parameters; in particular, the Coulomb coupling parameter Γ i,e and electron degeneracy parameter Θ e , * t.callow@hzdr.de † eli.kraisler@mail.huji.ac.il ‡ a.cangi@hzdr.de which are defined as where E kin and E pot are respectively the average kinetic and potential energies; τ the temperature; E F the Fermi energy; and the subscripts i and e respectively refer to nuclei and electrons. These parameters are of order unity in the WDM phase, Γ, Θ ∼ 1, which corresponds to a phase of matter lying somewhere between a classically ionized plasma and a strongly-correlated condensedmatter system. At typical WDM temperatures and densities, the electron coupling parameter, Γ e , is approximately equal to the density parameter r s [33], defined as [34] Γ e ≈ r s = 3 4πn e where n e is the number density of free electrons. We note that the Fermi energy for a non-interacting system of electrons can also be expressed using the above definition of r s [35], From the above relationship, and because E kin ≈ E F at low temperatures, the approximation Γ e ≈ r s (3) is valid in WDM (for r s ∼ 1). Finally, the classical Coulomb coupling parameter of nuclei is defined as with Ze the nuclear charge, and a i the mean inter-ionic distance. In Fig. 1, we show a rough schematic of some typical phenomena in the WDM regime, and the approximate region of WDM phase space which we later explore with our model. Notwithstanding the difficulties mentioned earlier, applying Kohn-Sham density functional theory (KS-DFT) [37,38] has recently led to promising results in the simulation of WDM [39,40]. KS-DFT is a well-established, successful method [41][42][43][44][45][46] for predicting the electronic structure of materials, from single atoms and small molecules to nanoparticles, periodic solids and large biomolecules. Within KS-DFT the interacting manybody problem is tackled by mapping it onto a fictitious, non-interacting system [38] which yields the same electronic density as the interacting problem. Usually, the systems are treated at zero temperature. The formal generalization of KS-DFT to finite temperature was established by Mermin [47]. In WDM simulations, KS-DFT is used to calculate forces on nuclei, which are then time-evolved through molecular dynamics techniques [40,48,49]. The primary target is to calculate the equation of state (EOS) which relates the atomic density, energy, temperature and pressure of a material. The EOS data can be used, for example, to compute the Hugoniot curve [40], which describes the possible final state from a given initial state after a shock wave, relevant to many WDM experiments. Furthermore, the electrical and thermal conductivities of WDM are calculated via Onsager coefficients [49] from KS orbitals and eigenvalues.
In principle, the KS method is an exact approach, but in practice, the exchange-correlation (XC) energy functional must be approximated. There is no systematic approach for the development of XC functionals, and thus a plethora of zero-temperature approximations exist; see, for example, Refs. 50-54 for discussions on this subject. In the extension to finite temperatures, functional con-struction is further complicated by the fact that, in principle, the XC functional should depend explicitly on the temperature [55]. However, the nature of this explicit temperature-dependence remains unclear and, in fact, is usually neglected in standard calculations [56,57]. Hence, the temperature dependence is only crudely included through the implicit temperature-dependence in the density, as the KS orbitals are occupied according to Fermi-Dirac statistics. The theoretical development of temperature-dependent XC approximations has recently found new momentum [57][58][59][60][61][62][63][64][65]. Investigations of the electron liquid [66] and uniform electron gas [67][68][69][70] provide insights into constructing local [56,[71][72][73] and generalized gradient approximations [74,75] to the temperaturedependent XC contribution. An exact inclusion of the exchange energy at finite temperature [76,77] has also been achieved using the optimized-effective-potential method, however, with the drawback of higher computational cost than that of the standard KS method. These are constructive steps towards developing accurate, reliable and computationally affordable temperature-dependent XC functionals.
The aforementioned computational cost is a further challenge for finite-temperature KS-DFT, because a large number of KS states must be accounted for under WDM conditions [78]. Despite this limitation, and the abovecited difficulties with developing suitable XC approximations, KS-DFT is currently the predominant firstprinciples method for simulations of large systems in this thermodynamic regime; alternative approaches such as orbital-free DFT (OF-DFT) [79][80][81] or path-integral Monte-Carlo (PIMC) [36,[82][83][84][85] tend either to not be sufficiently accurate (OF-DFT [86,87]) or too expensive (PIMC, especially at lower temperatures).
Consequently, the development of methods which obtain close to KS-DFT accuracy at reduced computational cost is an active area of research. Recently, there have been some promising developments in this area, such as the development of surrogate models using machine learning [88], stochastic methods [89][90][91], and approaches to reduce the cost of core-electron calculations [92,93].
In this paper, we consider an alternative approach known as an average-atom (AA) model: the premise of such a model is that the full system of interacting electrons and nuclei is partitioned into a set of Voronoi spheres, each containing a central nucleus, and the full electronic calculation is reduced (under certain approximations which we discuss later) to a calculation for a single atom. This concept, which has clear computational advantages, has a long history in plasma physics and electronic structure theory. The earliest AA models [94][95][96][97][98] were based on the Thomas-Fermi (TF) approximation [99,100] and modifications thereof; subsequent models built on this premise by adopting a mixed KS-DFT and TF approach for the bound and continuum electrons respectively [101], then treating the full spectrum (discrete and continuous) via KS-DFT [102], and later incorporating effects from outside the central atom such as ionic correlations [103][104][105][106]. AA models continue to be extensively developed and used under a variety of different approaches and assumptions [107][108][109][110][111][112][113][114][115][116][117][118][119][120][121][122][123][124][125]; we also mention here the recent works in Refs. 126 and 127 which attempt to bridge the gap between AA models and full KS-DFT via novel approaches.
AA models are a well established and successful tool in the plasma physics community as they produce results of useful accuracy at a fraction of the cost of full KS-DFT simulations. However, their evolution from the initial TF based models has been largely driven by organic developments, and there have been relatively few attempts to derive an AA model starting from a fully quantummechanical perspective (with some notable exceptions, for example Refs. 104, 105, 108, 112, and 117). As a result, AA models do not always follow the same conventions: for example, AA models which differ in their choice of boundary condition for solving the KS equations have been observed to yield quite different results [116]. By contrast, solid-state KS-DFT codes have differences in their numerical implementation, but generally follow a common set of theoretical assumptions and thus are expected to yield the same results for the same set of inputs (atomic configuration, temperature, etc) [128]. Establishing a similar framework and hierarchy for AA models would improve understanding of the limits under which they might be expected to give reasonable results in the WDM regime, in particular for τ E F /k B .
This motivates the main results of this paper, namely a systematic, first-principles derivation of an AA model starting from the many-body Hamiltonian of coupled electrons and nuclei, and the comparison of some fundamental results when the model is solved using KS-DFT. The paper is structured as follows: in Section II, we start with the full many-body Hamiltonian of electrons and nuclei and reduce it to an effective atomic Hamiltonian for a classical nucleus and surrounding electron density. We mention explicitly all the approximations that have to be done during the derivation and discuss their possible impact. We then introduce finite-temperature KS theory and apply it to minimize the grand free potential for this AA model, which requires particular consideration of the boundary conditions and treatment of unbound (continuum) electrons. Following a discussion on numerical implementation in Section III, we compute some common properties for a range of temperatures, densities, boundary conditions and XC functionals in Section IV. Finally, we discuss the implications of our results for AA models, and finite-temperature KS-DFT more generally, in Section V.

II. THEORY
In this section, we first reduce the many-body Hamiltonian of interacting electrons and nuclei to an effective single-atom Hamiltonian (Part II A), analysing the assumptions and approximations used in this reduction.
Then, in Part II B, we briefly review the grand canonical ensemble which can be used to describe a quantum system of electrons in thermal equilibrium with a reservoir (the nuclei); and in Part II C we explain how the problem of minimizing the grand free energy by solving for all the interacting states is greatly simplified by finite-temperature KS-DFT. Next, in Part II D, we apply KS-DFT to the reduced average-atom Hamiltonian and deduce suitable boundary conditions for the electron density within this model. Finally, in Part II E, we discuss our treatment of the unbound electrons and how the boundary condition on the density can be realised through a number of different boundary conditions for the bound KS orbitals, from which we consider two simple choices.
A. Many-body Hamiltonian of coupled electrons and nuclei in the dilute gas limit We begin with the full many-body Hamiltonian of interacting electrons and nuclei, The individual components of the Hamiltonian are defined aŝ In the above, we have assumed that the system is composed of a single element with nuclear mass M , charge Z and electron number N e . N n denotes the number of nuclei in the system,thus N n · N e is the total number of electrons. r i and R I are the positions of the i th electron and the I th nucleus respectively. We have also assumed that there is no external field applied to the system. Note that here and below we adopt Hartree atomic units, = e = m e = a 0 = 1.
We work within the Born-Oppenheimer (BO) approximation [129], which assumes the electrons react instantaneously to any changes in the positions of the nuclei due to their relatively small masses, m e M . In fact, though the BO approximation is used extensively in WDM simulations with KS-DFT [1], it is likely to be prone to inaccuracies in the WDM regime, due to strong non-adiabatic effects from excited states, core electron chemistry, and so forth [130,131]. While we persist with the BO approximation in our derivation, we discuss some possibilities for incorporating non-adiabatic couplings between electrons and nuclei in Section V.
Having fixed the nuclear co-ordinates in the Hamiltonian with the BO approximation, we transform the vectors r i (which act on the electron wavefunction) as r i = R I + x iI , where x iI is determined by performing a Voronoi decomposition of space. In other words, each vector r i is now defined relative to the closest nuclear co-ordinate R I . We wish to make clear that this is sim-ply a relabelling of the terms in the Hamiltonian and the notion of a "closest nucleus" does not imply any assumptions regarding the electron density distribution. Following this transformation, and ignoring the nuclear kinetic energy to the BO approximation, the Hamiltonian can be re-written aŝ where N (I) e denotes the number of vectors x iI closest to the I th nucleus. Next, we decompose the electron-nuclear and electron-electron interactions in the above expression into those parts which contain interactions between charges in a single Voronoi cell, and those involving inter-cell interactions. The Hamiltonian (10) thus becomeŝ Next, we consider the decomposition of the Hamiltonian into two parts: an "average" termĤ av , and an "inhomogeneous" termĤ in , i.e.Ĥ =Ĥ av +Ĥ in . The average componentĤ av is constructed by considering the expectation value Ψ|Ĥ|Ψ (where Ψ is an anti-symmetric wave-function), for the particular case in which the nuclei are distributed exactly evenly in space. In this case, the electron density is identical in each of the Voronoi partitions, and therefore the expectation value Ψ|Ĥ|Ψ is equal to the expectation value of the average Hamiltonian Ψ|Ĥ av |Ψ . We then make two further assumptions regarding the average Hamiltonian. Firstly, we assume that each nucleus is associated with the same number of vectors r i following the Voronoi decomposition, which means that N (I) e = N e . Secondly, we assume that the vectors x iI are independent of any nuclear positions; in other words, x iI can be written simply as x i . By further transformation of the co-ordinate system such that the I th nucleus lies at the origin, R I = 0, the average Hamiltonian becomesĤ The inhomogeneous term,Ĥ in is the difference between the full HamiltonianĤ of Eq. (11) andĤ av of Eq. (12). In our derivation, we choose to neglect it. This is a reasonable assumption if the nuclear distribution is relatively uniform (for the fictitious system considered in the previous paragraph in which the nuclear distribution is exactly uniform, Ψ|Ĥ in |Ψ = 0). The contribution of H in for non-uniform systems can in principle be taken into account via perturbation theory. As a conceptual remark, we note that the average Hamiltonian (12) is akin to the Hamiltonian in a typical KS-DFT simulation for a periodic system, in which some electrons in a unit cell (here containing a single nucleus) interact between themselves and their periodically repeating images.
We shall now split the re-formulated Hamiltonian (12) into two parts as follows, where The first componentĤ at el (14) is equivalent to the electronic Hamiltonian of a single atom and we shall return to it later in the paper. We first treat the termsŴ J (15) which make up the second component of the Hamiltonian (13). We shall expand perturbatively in powers of |y iJ |, since by the Voronoi decomposition of space, the inequality |y iJ | ≤ 1 2 holds strictly, because the distance between an electron and its nearest nucleus cannot exceed half the distance between two nuclei. Moreover, for electrons located in cells far from the central cell, and for electrons tightly bound to the central nucleus, |y iJ | is much lower than 1, which further justifies a power expansion.
We re-writeŴ J in the form and expand both terms in the above expression using the binomial expansion for 1/ √ 1 + . A full derivation can be found in Appendix A. The zeroth-order term in this expansion is equal tô where the final equality holds for systems with neutral charge, N e = Z. We now consider the first and second order terms in the expansion ofŴ J , which turn out to be generally non-vanishing and equal tô where we introduced the notation In principle, one could continue this perturbative expansion to include even higher order terms. However, in our AA model, we neglect all the coupling termsŴ (k) J , with k 1, leaving only an atomic Hamiltonian. Considering the interaction terms we neglect in our model already provides insight regarding the limits under which we might expect it to be accurate. Of course, we would expect high accuracy in the limit of a dilute gas (|y iJ | 1). Moreover, for an approximately uniform nuclear distribution, which implies a highly symmetric electronic distribution, we would expect the expectation value Y J to be close to zero. This suggests the model is likely to even be accurate in the high-density limit, in which the Voronoi cells and their enclosed electronic distribution will be spherically symmetric to a large extent. We can also expect accurate results at high temperatures, when there is significant ionization and the electron kinetic energy dominates over interaction effects. In addition, the detailed derivation performed here shows how the AA model can be made more accurate in the future, by including higher-order terms ofŴ J and perturbatively treating the inhomogeneous term,Ĥ in .
We have therefore reduced the full many-body Hamiltonian of electrons and nuclei, defined by equations (6)(7)(8)(9), to N n ×Ĥ at el , whereĤ at el is the Hamiltonian for a system of electrons interacting with a single fixed nucleus. We note that, although we ignore inter-cell interaction terms in the Hamiltonian, these interactions will be partially accounted for by the choice of boundary conditions; this is what distinguishes this model from a single isolated atom. In the literature, such a reduction is commonly referred to as an AA model: in neglecting completely thê W J terms, our model falls within a class of AA models known as ion-sphere models [101,109,114,116]; other socalled ion-correlation models [106,119] attempt to model to some extent these coupling contributions via the introduction of a uniform background field [102,103,120,121]. Under the umbrella of ion-sphere models there exists an abundance of models, differing for example in how the atomic Schrödinger equation is solved for the bound and unbound electrons, choice of boundary conditions and more besides; we shall draw comparisons with appropriate examples of such models later in this paper.

B. Finite-temperature quantum systems
We are interested in applying the above model to systems at finite temperatures; such systems are described by a statistical ensemble of states. In this sub-section, we review some basic theory regarding quantum statistical ensembles.
We restrict our analysis to systems at some fixed temperature in thermal equilibrium with a reservoir, which (like the nuclei) is treated classically. In such an ensemble, known as the grand canonical ensemble, the grand canonical HamiltonianΩ and associated grand canonical potential (or grand free energy) Ω play the role of the Hamiltonian and energy in a zero-temperature calculation,Ω whereŜ andN are the entropy and number particle operators;Γ is a statistical operator which can be used to determine average observable values of operators; and µ is the chemical potential, defined as the change in energy when a particle is added or removed from the system [132].Ŝ andΓ are defined bŷ where |Ψ k,N are orthonormal eigenstates of the operator Γ , where k is the principal quantum number and N is the number of electrons in the system. w k,N are statistical weights, which satisfy k,N w k,N = 1. One example for the use of the operatorΓ is for finding the ensemble density: wheren(r) = N i=1 δ(r − r i ) is the density operator and n k,N (r) = Ψ k,N |n(r) |Ψ k,N is the density of the (k, N ) th state.
When the system of electrons and nuclei is in thermal equilibrium with the reservoir, the grand canonical potential is minimized with respect toΓ. In this case the equilibrium weights w k,N are given by with E k,N = Ψ k,N |Ĥ |Ψ k,N the equilibrium eigenvalues of the states |Ψ k,N , β = 1/k B τ , and Z is the partition function, We also note that the grand canonical potential is often expressed in terms of the partition function, C. Finite-temperature KS-DFT Whilst there are many first-principles techniques to determine ground-state electronic properties, the majority of these are not computationally feasible at finite temperatures, even for atomic systems. We shall use KS-DFT as it has the best balance between accuracy and speed, particularly in the low-temperature part of the WDM regime. In its original formulation [37], DFT establishes a one-to-one mapping between the ground-state electronic density and external potential v ext (r); this means the ground-state density is (in principle) sufficient to compute all observables. Mermin [47] extended the DFT formalism to ensembles at finite temperatures; in this case, there is a mapping between the equilibrium ensemble density n 0 (r) and the external potential minus the chemical potential, v ext (r) − µ, and hence n 0 (r) can be used to compute observable averages. The grand canonical potential Ω = Ω[n] is thus a functional of the density [58], where T τ [n] denotes the electronic kinetic energy at temperature τ , V τ ee [n] the electron-electron interaction energy, and v ext (r) the external potential (which in our case is just the electron-nuclear attractive field). F τ [n] is denoted the universal functional because it has no dependence on the external potential; minimizing Ω[n] with respect to the density yields the equilibrium free energy.
Like in the ground-state case, it is convenient to introduce an auxiliary system of non-interacting electrons which has the same density and temperature as the fullyinteracting system. This is known as the KS system [38,55]. The definition of the grand free energy within KS theory is given by where v ext (r) is the external potential experienced by the electrons (in our model, the electron-nuclear potential); and F τ xc is the exchange-correlation (XC) free energy functional, which is defined as and U [n] is the usual Hartree energy, given by The exact XC free energy functional contains all the information about electron-electron interactions. The KS orbitals, and hence the density, are obtained by solving the finite-temperature KS equations, where v H (r) is the Hartree potential, given by the functional derivative of the Hartree energy (34) with respect to the density. In the above, we have moved from a pure density-functional formalism to a spin-density-functional formalism [133], which allows for a formal treatment of magnetic systems, and often yields more accurate results even in the absence of an external magnetic field [134]. The variable σ =↑, ↓ represents the spin-channel. Strictly speaking, we note that the KS potential, and thus by extension the KS eigenvalues, should depend on the chemical potential due to the aforementioned mapping between the density and the external potential minus the chemical potential. However, since this makes no practical difference in the static case aside from a shifting of the orbital energies, we adopt the more common convention seen in the equations above in order to better connect with existing AA and KS-DFT literature. However, if one wanted to extend the model to compute linear response quantities (for example), the formally correct definition for the KS potential and orbital energies should be used. The notation ψ σ i, (r) denotes respectively bound orbitals ψ σ i (r), with discrete energy levels τ iσ , and unbound orbitals ψ σ (r), with continuum energies τ σ . The spindensities are given by where f iσ denotes the occupation number of the i th energy level τ iσ , f σ ( ) the continuum distribution function, and g σ ( ) the density of states. The total density is simply equal to the sum over the spin-densities, n(r) = n ↑ (r) + n ↓ (r). Since the KS system is a noninteracting system of electrons, the occupation numbers are determined by the Fermi-Dirac (FD) distribution, with f σ ( ) defined in a similar way. As discussed, finite-temperature KS-DFT is formally an exact method; however, the XC functional F τ xc [n] has to be approximated in practice. Formally, F τ xc [n] should depend explicitly on the temperature (a dependence which is absent in many approximations), and satisfy various exact conditions [58]. Furthermore, even XC approximations which satisfy these requirements can suffer from various errors. We outline two such errors, particularly relevant in the context of this work, below: (i) Self-interaction error : The classical Hartree energy (34) in ground-state KS-DFT introduces a spurious selfinteraction (SI) -repulsion of the electron from its own charge density [135] -which must be compensated by the XC term. However, for many approximate XC functionals the self-interaction is not compensated fully, and one is left with a self-interaction error (SIE). A prototypical example is the H atom, for which approximate-XC calculations can be compared against an analytic result, quantifying the SIE. Self-interaction causes an array of problems, among which is also the under-prediction of ionization potentials [136][137][138][139]. Development of methods to mitigate the SIE is a very active area of research [138,[140][141][142][143][144][145][146][147][148]. The SIE is also present in finitetemperature KS-DFT and also there it must be compensated by the XC functional; finite temperature makes this task even more difficult.
(ii) Ghost-interaction error : In finite-temperature KS-DFT, the density is constructed as a weighted ensemble of densities of KS determinants (see Eq. (26)). Formally, the KS states which form this ensemble should not interact with each other; however, the Hartree energy is a functional of the total density and thus a ghostinteraction error is present due to repulsion between electrons in different KS states [149]. This error is less studied than the SIE, but some correction schemes have been proposed when the ensemble weightings are defined a priori [150,151].

D. Reduction to average-atom model in finite-temperature KS-DFT
We shall now apply the finite-temperature KS scheme detailed above to the reduced Hamiltonian in Eq. (14). Using this Hamiltonian, the computational cost of solving the KS equations is significantly reduced, because the finite-temperature KS equations (35) need only to be solved for a single atom. In the single-atom picture, the electron-nuclear potential is spherically symmetric. The same is not necessarily true for the KS potential defined in Eq. (36); however, we make the spherical approximation (which is usually made in atomic KS calculations), in which the KS potential is assumed equal to its spherically-averaged value. This means the bound KS orbitals ψ σ i (r) can be decomposed into a product of radial and angular components, where the angular component Y m l (θ, φ) is a spherical Harmonic function. The eigenstates are now characterised by the quantum numbers n, l, m, and σ, where n denotes the energy level, l and m characterise the orbital angular momentum, and σ refers to the spin channel.
The radial component X σ nl (r) is determined by solving the one-dimensional differential equation given by where r > (x) = max(r, x). All integrals are performed within the Voronoi sphere, whose radius R VS [152] is determined from the average density of the nuclei n i , In the WDM regime, i.e. at high pressures and temperatures, the unbound electrons play an important role. Due to the nature of continuum states [153], it is difficult to compute the unbound density exactly and thus it is determined (in our model) in an approximate manner; we discuss this further in the following section. For now, we note that, without further approximation, the total (spin) electron density is split into bound and unbound components, with the bound component given by where f σ nl are the Fermi Dirac occupations, as defined by Eq. (38), and with the KS orbitals normalized within the Voronoi sphere, Since we adopt the formalism of spin-DFT, the number of electrons in each spin-channel is fixed to some integer value N σ e , which also fixes the total number of electrons in the Voronoi sphere equal to its average value N e = σ N σ e . The (spin-dependent) chemical potential µ σ in the FD distribution is thus determined from the condition that the number of electrons in each spin channel is fixed. In theory, the values of N σ e should be determined by whichever configuration minimizes the grand free energy [154].
We make a clarification here regarding the fact that there are two chemical potentials, µ ↑ and µ ↓ . This is a necessary consequence of using spin-DFT; however, one could say that in the absence of a magnetic field, there should be only a single chemical potential for the electrons. In fact, this is not a problem, because only one of the chemical potentials is physically relevant in this scenario: when adding an electron, the lower chemical potential is relevant, and when removing an electron, the higher chemical potential is relevant. If non-spindependent DFT were used, as is quite common in AA models, then only the total electron number is conserved and a single chemical potential returned. However, it is well known in DFT literature that using spin-dependent DFT is advantageous, because it produces more accurate results for systems whose spin is not zero [134,155,156], so we use the spin-dependent formulation in our model.
In this KS-AA model, we impose the condition that the electron density is smooth at the boundary between neighbouring spheres. Physically, this is motivated by the fact that the real electron density (which is formally equal to the KS density) should be smooth everywhere. Of course, the true system cannot be split into identical spheres, so this condition is an approximation designed to mimic the real electron density. Mathematically, this means the following boundary condition should be imposed on the (spin) density, The above condition on the density is physically intuitive, and leads to boundary conditions on the radial KS orbitals at the Voronoi sphere's edge that are widely used in AA models. However, it is not a necessary condition, and alternative choices can be made to model the concept of an atom immersed in a plasma. For example, in Ref. 120, there is no boundary condition applied to the density or KS orbitals at the sphere's edge; instead, the potential is fixed to a constant value outside the sphere ) and the KS orbitals are solved out to some radius r max R VS . A similar approach (with the KS potential instead modified inside the Voronoi sphere) is used in the MUZE code [119], and we shall compare results obtained with that approach with the above condition on the density in Section IV C.

E. Boundary conditions and treatment of unbound electrons
As mentioned in the previous sub-section, we do not explicitly solve the KS equations for the continuum orbitals in our model. Instead, we treat the unbound electrons as being completely free (the ideal approximation), i.e., having uniform density, when we compute their contribution to the total density. For this distinction between bound and continuum states, we assume the continuous part of the energy spectrum starts at v τ,σ s (R VS ). This is equivalent to shifting the KS potential everywhere by the constant v τ,σ wherev τ,σ s (r) signifies the shifted potential, and then assuming the continuum starts at energies above zero (sincē v τ,σ s (R VS ) = 0). We emphasize that shifting the potential by a constant has no effect on the KS orbitals, only their energy eigenvalues which shift by the same constant v τ,σ s (R VS ). Since it makes the notation easier, we assume the potential has been shifted in this way for the rest of the paper, unless specified otherwise.
We note that, in modern AA codes, the unbound electrons are usually treated in a more sophisticated manner, either semi-classically (TF) or (more typically) in a fully quantum manner [102,110,111], for example by expanding the continuum states in a discrete set of normalizable states [157,158]. The TF approximation for the unbound density is known to have certain limitations, for example systematically over-estimating the chemical potential [119]; these limitations are likely to be exacerbated using the even simpler ideal approximation. Nevertheless, our model will yield important insights into the comparison of different XC functionals and boundary conditions. We also compare this ideal treatment with TF and quantum unbound electrons in Section IV. Furthermore, we clarify that the density used to construct the KS free energy functional (32) and potential (36) is the full density (bound and unbound), with the kinetic energy of the unbound electrons being given by the ideal expression.
The value of the unbound electron (spin)-density n σ ub is determined from the number of unbound electrons, VS is the volume of the sphere. The total numbers of bound and unbound electrons are determined according to the FD distribution, We recall that the chemical potentials µ σ are chosen such that the sum N σ b + N σ ub equals a preset value, N σ e . We now proceed to discuss the question of boundary conditions. The boundary condition for the density has been specified by us in Eq. (46). However, to solve the KS equations (40), one needs to specify the boundary conditions for the radial orbitals X σ nl (r). Notably, Eq. (46) does not determine the orbital boundary conditions uniquely. For the unbound density, Eq. (46) is satisfied automatically. For the bound states, it implies the following equality for the radial orbitals X σ nl (r) and their derivatives: There is no unique way to satisfy the above relation. The two most simple choices are either to require the radial wavefunctions to be zero at the boundary, or their derivatives to be zero, i.e.
Mathematically, many other choices are also possible. Both of the above boundary conditions have been used in AA models (see for example Ref. 114 for an example of the former or Ref. 116 for an example of the latter), and it has been observed [116] that using the former choice (51) yields markedly different results than the latter (52) for the average ionization in Aluminum. We also observe that the choice of boundary conditions has a major impact on results in Section IV of this paper. We additionally note that we will sometimes refer to the boundary conditions in Eqs. (51) and (52) as b.c. (i) and b.c. (ii) respectively, in order to simplify notation. In spite of this, there has been limited analysis of the impact of choosing one of the above conditions over the other (perhaps because many models, instead or additionally, enforce boundary conditions on the potential as earlier discussed). However, Rozsnyai has conceptually identified [101,106] these boundary conditions as corresponding to upper (51) and lower (52) limits of a bandstructure, due to their association with bonding (52) and anti-bonding (51) molecular orbitals (MOs); Massacrier and co-workers have explored this band-structure interpretation further, interpolating between the bandstructure limits via a Hubbard functional form for the density of states [125,159,160].
In Fig. 2, we compare the density distribution resulting from the two different boundary conditions for Beryllium with density 0.266 g cm −3 and temperature 4 eV. We find that near the origin (r < 1 a 0 ) the densities are the same, whereas close to the boundary, they differ significantly. Notably, using condition (51) leads to a higher unbound density: the requirement of the wavefunction vanishing at R VS "pushes" the electrons to the unbound states.
This completes our first-principles reduction from the intractable minimization of the grand canonical potential (23) for a fully-interacting system of electrons and nuclei (6) into a finite-temperature KS model for a single atom. In the remainder of this paper, we briefly discuss some numerical aspects and present some results to explore the behaviour of this AA model under different approximations.

III. NUMERICAL ASPECTS
The majority of calculations were performed with the ORCHID code, an atomic DFT code which has been extended to include temperature and the boundary conditions described in Section II E; ORCHID has been used previously in Refs. [161][162][163][164][165]. An open-source successor to ORCHID, named atoMEC, is under development and can be used to reproduce many of the results in this paper [166]. For the bound KS orbitals, the second-order differential equation (40) is solved on a grid using the Numerov method in combination with a two-side shooting method [167][168][169]. The grid uses a logarithmic scale for the radial co-ordinate r, x = ln(Zr), where Z is the atomic number, to ensure a higher density of grid points near the nucleus. To solve for the radial wave-functions X σ nl (r), a transformation P σ nl (x) = X σ nl (x)e x/2 is made, so that the differential equation (40) becomes where the notation¯ τ ,σ nl denotes the KS eigenvalues corresponding to the shifted KS potentialv τ,σ s (47). Additionally, some calculations (involving different treatments of unbound electrons, to be later discussed) were performed using the MUZE code [119,170,171]. Below, we discuss details of the implementation of our model in ORCHID, and refer readers to the aforementioned references for further details of the MUZE code.
The KS orbitals and their eigenvalues from the differential equation (40) are used to determine the total density via Eqs. (44), (48) and (49). The whole process proceeds self-consistently until the total energy is converged to within 10 −6 Hartree, the spatially-averaged bound spin-densities are converged enough to satisfy dr(n σ (r) − n f σ (r)) < 10 −6 a.u. (the superscript f distinguishes the density of the previous iteration from the density of the last iteration) and the KS potentials (for both spin channels) satisfy dr(v S,σ (r) − v f S,σ (r)) < 10 −6 a.u. Different guesses were trialled for initializing the KS potential in the self-consistent cycle; it was found the bare Coulomb potential is always a suitable initial choice for the KS potential.
For each atom, XC functional, density and temperature range [172], convergence was checked with respect to the number of grid points N grid , the leftmost gridpoint x 0 = ln(Zr 0 ), and which bound states (characterised by the quantum numbers n, l) to account for. The required grid size depends strongly on the specific calculation, ranging between 1,000 to 40,000 points. The number of bound states is also highly dependent on the particular calculation. However, a lower grid bound of r 0 = e −13 /Z a 0 was found to be sufficient in almost all cases.

IV. RESULTS
In the following, we use our KS-DFT AA model as a surrogate of AA models to probe their sensitivity to various choices of approximation. We focus mostly on the impact of XC functional and boundary conditions, but also briefly explore the influence of more advanced treatments of unbound electrons in Section IV C.
We probe a range of temperatures from 0 to ∼ 25 eV, and densities from roughly 10 −3 to 20 g cm −3 , corresponding to to Voronoi sphere radii in the range 2 a 0 ≤ R VS ≤ 10 a 0 ; this roughly covers the lower-density and temperature region of the WDM phase space, as denoted in Fig. 1. For the temperature-dependent plots which are shown at fixed values of R VS , we connect these Voronoi sphere radii with their mass densities ρ m , and the corresponding dimensionless parameters r s , Γ i , and Θ e (the latter two computed at the mid-point of the temperature range we consider, τ = 13.1 eV) in Table I.
We note that the definitions of r s , Γ i , and Θ e in Eqs. (2)-(5) use as input the free electron density, and furthermore that the definition for the Fermi energy is strictly valid for an ideal Fermi gas only. However, we of course deal with interacting electrons, some fraction of which are bound by the nuclei. For the numbers quoted in Table I, we assume for simplicity that the free electron density is a fixed quantity for a given ionic density (independent of temperature) and equal to the valence electron density; we also use the definition in Eq. (4) for the Fermi energy regardless. These numbers should therefore be seen as rough indicators only rather than well-defined parameters. Based on these values, we note that the lowest density Hydrogen example (R VS = 10.0 a 0 ) is at the limit of what is typically considered WDM conditions; however, it is a useful test case for functional comparison as the AA approximation (meaning the neglect of intercell interactions) should be very accurate in this case.
From our KS-DFT AA model, we directly obtain the free energy, KS orbitals and their energies, occupation numbers, and the number of unbound electrons. To access the electronic pressure P e , we use the following relationship between the electronic free energy F and the volume V for a fixed temperature τ , which we compute numerically via finite differences. The free energy F is defined from the grand free potential Ω of Eq. (32) via the relationship Please see Appendix B for details of the construction of the free energy F [n] in our AA model. We focus on the electronic pressure because KS-DFT does not give access to the ionic pressure. We have explored adding c. X nl (R VS ) = 0 (51) and on the right are for the b.c. X nl (R VS ) = 0. Insets: percentage difference in pressure (using a log scale for the y-axis) between the approximate functionals and the exact functional. The x-axes of the inset plots are identical to the (temperature) x-axes of the main plots. Smaller plots below main ones: actual difference in pressure between approximate and exact XC functionals, ∆P e = P app e − P exact e . an approximate ionic pressure using the ideal gas law, pV = nRT , and observed this results in a noticeable increase in pressure. However, since we want to explore the impact of approximations which have no effect on the value of the ionic pressure at a given temperature and density, we only present results for the electronic pressure.
Furthermore, the number of unbound electrons, or equivalently (in our model) the mean ionization state (MIS), is an important property in dense plasmas [119], as are the KS orbital energies, which are used (for example) in the computation of thermal and electrical conductivities (Section I, Ref. 49). We therefore focus on the aforementioned quantities.
We compare results for three XC functionals: firstly, the zero temperature local spin-density approximation (LSDA) [38,173], which is widely used in finitetemperature KS-DFT and AA models; secondly, the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation [174], which is used extensively in groundstate DFT calculations [175]. Thirdly, we consider the temperature-parameterized LDA by Groth et al. (GDSMFB) [73]; this functional retains the computational advantages of the ground-state LDA, and can easily be integrated into other AA or KS-DFT codes via the LIBXC package [176], as has been done here in ORCHID. We also investigated the temperatureparameterized LDA by Karasiev et al. [56], and found that results were always in very close agreement with the GDSMFB functional (echoing similar observations in Refs. 177 and 178), and therefore we present results for LSDA, PBE and GDSMFB only.  In the following sections, we explore results for Hydrogen (Section IV A) and Beryllium (Section IV B); although the model we have derived is valid for systems with a macroscopic net charge, we consider only charge neutral examples (N e = Z). We drop the spin-dependent notation for all spin-dependent quantities such as the KS orbitals and their eigenvalues. Beryllium has an even number of electrons and thus the different spin channels give the same results. For Hydrogen, one of the spin channels is completely devoid of electrons (amounting to a constraint of N σ e = 0): since this spin channel does not contribute to many physically meaningful quantities, where we present results for spin-dependent quantities, these are for the occupied spin channel only.

A. Hydrogen
We first apply our AA model to Hydrogen. Besides being an element of high interest in the WDM regime in its own right [178], for Hydrogen we can solve our AA model exactly. Since we consider interactions only within the Voronoi sphere, which itself contains only one electron, there are no inter-electron interactions and thus the KS potential is given exactly by the electron-nuclear attraction, The XC functional also cancels the Hartree energy, which is the exact XC functional in this case. We stress here that this is the exact result within the limits of the model we define, in which various approximations have already been made, such as the neglect of explicit intra-cell interactions; it therefore does not represent the truly exact limit for the Hydrogen plasma in general. Another assumption that was mentioned but not discussed at length is that we take N ↑ e = 1, N ↓ e = 0. We fix this for convenience, but in principle one could search over all fractional N ↑,↓ e , with N ↑ e +N ↓ e = 1, and choose the configuration which minimizes the grand free energy. However, the exact XC functional is a reference from which we can isolate the errors that result from approximations for the XC functional (as opposed to other approximations in the model). Henceforth, we use the term "exact" when referring to this choice given the fact that we mean only to the exactness of the XC functional and potential.
In Fig. 3, we plot the P e -τ curve for Hydrogen at different values of R VS . The most notable observation in this series of the plots is that the choice of boundary condition has a significant impact, particularly at lower temperatures and higher densities. Indeed, for R VS = 2.0 a 0 , the EOS data is qualitatively different, with the boundary condition X nl (R VS ) = 0 showing large negative pressures for all the functionals (including the exact one) at lower temperatures. The approximate functionals show good agreement with the exact result, rarely differing by more than 10% and tending towards < 1% for higher temperatures. An additional observation is that the LDA and PBE functionals yield almost identical results throughout, with the temperature-dependent GDSMFB functional tending to deviate slightly more from the exact result.
Next, in Fig. 4, we compare the number of unbound electrons N ub as a function of temperature for two choices of density. For these (relatively low) densities, the two boundary conditions seem to be in quite good agreement. What is more interesting is that the approximate functionals tend to systematically over-predict the mean ionization state relative to the exact functional; furthermore, the approximate functionals differ minimally relative to each other. This is indicative of some kind of common error pertaining to semi-local XC functionals in general.
In Fig. 5, we consider the 1s energy level¯ τ 1s as a function of temperature for R VS = 4.0 a 0 . It is apparent that the approximate XC functionals systematically overpredict the 1s energy level (for both b.c.s), with once again minimal differences between the functionals themselves. The energy levels for the exact (bare Coulomb) XC functional are of course independent of temperature since v en (r) has no temperature dependence; interestingly, the energy levels from the approximate functionals also do not vary significantly across this temperature range (with an exception for very low temperatures for the b.c. X nl (R VS ) = 0). An important observation from these curves is that the temperature-dependent GDSMFB functional does not seem to improve much the prediction of the¯ τ 1s level: therefore, new XC functionals (most likely going beyond semi-local approximations) are required for finite-temperature KS-DFT.
Having explored the behaviour of various quantities as a function of temperature for fixed mass density, we now investigate the dependence of the same quantities on the mass density, for a fixed temperature 10 eV. We only compare results from the LSDA and exact XC functionals, since we already observed that all the approximate functionals gave very similar results; we also compare the  two boundary conditions directly in the same plots.
In Fig. 6, we see the striking impact of the boundary condition on the electronic pressure as the density increases. As would be expected, for lower densities, when there are fewer interactions between neighbouring atoms, the two boundary conditions agree relatively well, with the difference between them (shown in the middle panel of Fig. 6) rarely exceeding 10% up to around 0.04 g cm −3 . However, at higher densities, the two boundary conditions diverge strongly; the boundary condition X nl (R VS ) = 0 actually has a turning point at which the pressure starts to decrease with increasing density and becomes strongly negative. The bottom panel of Fig. 6 shows the difference between the functionals as The main message from this figure is the huge impact of the boundary condition, and the limitations of choosing a single boundary condition as is frequently done in AA models. As discussed earlier, one possible solution is to consider a band-structure picture, either via some sensible approximation such as that employed by Massacrier and co-workers [125,159,160], or even better using a first-principles approach which preserves the smoothness of the density at the sphere's edge (46). Another quite striking feature of Fig. 6 (also seen in Figs. 4 and 7) are sharp discontinuities, which appear due to the ionization of the 2s energy level at around 10 −2 g cm −3 for the b.c. X nl (R VS ) = 0, and the ionization of the 1s energy level at around 10 −1 g cm −3 for the b.c. X RVS nl (r) = 0. We shall later see similar discontinuities for Beryllium in Section IV B, when we will discuss them in greater detail.
The divergence of the pressure towards negative infin- ity for the X nl (R VS ) = 0 boundary condition is related to the ionization degree. In Fig. 7, we plot the number of unbound electrons (top) and the 1s energy level¯ τ 1s as a function of the mass density, again for τ = 10 eV. As was observed for the electronic pressure, results from the two boundary conditions diverge with increasing density. In particular, the¯ τ 1s level for the boundary condition X nl (R VS ) = 0, for both LSDA and the exact XC functionals, has a turning point and starts to decrease for densities above about 0.1 g cm −3 . This effect was also observed in Ref. 125 for Aluminium and Carbon. As a result, the number of unbound electrons with this boundary condition actually decreases (and seems to be approaching a value of zero) with increasing density. This results in decreasing pressure, because (for example) the kinetic energy decreases with a lower ionization degree, and thus the free energy decreases as the mass density increases. It should therefore be noted that the extreme divergences in pressure between the two boundary conditions would likely be suppressed to some degree if a different definition of pressure is used, which is not so directly influenced by the ionization degree. Regardless, decreasing pressures and ionization degrees with increasing densities is a counter-intuitive and seemingly unphysical result, since we expect greater ionization and pressures at higher densities. Furthermore, the behaviour of the 1s energy level raises important questions related to the concept of ionization potential depression (continuum lowering), a critical effect in materials under WDM conditions [179].
It is well documented in both experiments [180][181][182] and theoretical models [183][184][185] that ionization potentials -defined as the energy required to excite a given bound electron into the continuum -are lower for atoms immersed in a plasma relative to the isolated atom case, though there remains uncertainty regarding the precise nature of this effect [179,186]. In KS-AA models, it is typical to associate the KS orbital energies with the actual electronic energy levels, and by extension (as the continuum levels are usually defined as those with positive energy, > 0) the orbital energies define the ionization potentials. As an aside, we note that there is in fact no formal relationship between the KS orbital energies, which belong to a fictitious system of non-interacting electrons, and the actual electronic energy levels, with the exception of the HOMO level whose (negative) value is equal to the ionization potential in ground-state KS-DFT only [187][188][189][190]. Nevertheless, it has been postulated that the KS orbitals are a reasonable surrogate for the molecular orbitals of the real interacting system [191,192], which justifies to some extent the association of the KS orbital energies with the real electronic energy levels.
In light of the above, the density dependence of the 1s energy for the boundary condition X nl (R VS ) = 0 seems to be a strange result. This, and likewise the behaviour of the number of unbound electrons, points to the limitations of models such as our own which only take into account screening effects from the surrounding plasma in a coarse manner through boundary conditions on the density or potential, neglecting all explicit interactions between charge densities in the central sphere and its neighbours. In this sense, it is possible that the approximate XC functionals actually benefit from error cancellation within the model we define, as their errors relative to the exact XC functional may be partially compensated for by the opposing error induced by the neglect of intercell interactions. However, this error cancellation is not omnipresent (it will not occur in the low-density limit for example); moreover, when more advanced models are constructed which include to a greater extent the effects of inter-cell interactions, the approximate XC functionals will no longer benefit from this error cancellation.
Although it is clear that the choice of boundary condition typically far outweighs any error than the choice of XC functional, especially at higher densities, it is worth analysing in more detail the functional errors. The systematic under-binding of the electron density and overprediction of the 1s energy level are both related to the SIE discussed in Section II C. We further analyse the SIE in Fig. 8. Since there is only one electron in Hydrogen, the exact XC energy should exactly cancel the Hartree energy and therefore the ratio −F xc /U H should equal 1.
Interestingly, it appears that the GDSMFB functional is contaminated by a larger SI error (which increases with temperature and decreases with density) than the zerotemperature LSDA functional. This may explain why the GDSMFB functional seems to yield slightly larger errors, relative to the exact XC result, for the electron pressure in Fig. 3. Of course, the SI error is particularly important for Hydrogen in the AA model, and may be overwhelmed by other factors in different examples. Nevertheless, this figure does not explain the failings of the LSDA functional relative to the exact result for smaller values of R VS , since the ratio −F xc /E xc does not significantly deviate from unity.
One of the ramifications of the SI error is the incorrect asymptotic behaviour of the KS potential, which (among other factors) contributes to the electron density being too delocalized; this error, known as the delocalization error, is ubiquitous to (semi)-local XC functionals in DFT [193,194]. In Fig. 9 we show an example of how both the LSDA and GDSMFB potentials differ from the exact result,with most notable differences in the asymptotic region in which they decay incorrectly. In Section V, we consider some possibilities to mitigate against SI and delocalization errors. Inset: potentials plotted against 1/r to emphasize different asymptotic behaviours (1/r → 1/R VS ).

B. Beryllium
We now apply our AA model to the Beryllium atom. Beryllium is used in ICF capsules [4] and relevant to astrophysical processes [195], and thus accurate simulations of Beryllium under WDM conditions are of high interest [196,197]. Although there is no benchmark for the XC functional as in the case of Hydrogen, it is interesting nevertheless to compare choices of boundary conditions and XC functionals.   In Fig. 10, we plot the pressure P e as a function of temperature τ , for two values of the Voronoi sphere radius which correspond to mass densities of 3.04 and 0.379 g cm −3 respectively (for reference, the ambient solid density of Beryllium is 1.85 g cm −3 ). Under these conditions, it seems the pressure for the two different boundary conditions is in relatively good agreement, though more significant for the lower density with R VS = 4.0 a 0 ; furthermore, the LDA and GDSMFB functionals also agree very closely for both boundary conditions, with no observable differences between the functionals for the lower density. Next, in Fig. 11, we plot the mean ionization state N ub for Beryllium with density 0.379 g cm −3 as a function of temperature, again showing both boundary conditions and the LSDA and GDSMFB functionals in the same plot for comparison. Here, particularly at low temperatures, we see significant differences due to the boundary conditions, but the choice of functional has very little impact. The large deviation in the mean ionization state due to the choice of boundary conditions is explained by the eigenvalue spectrum. In Table II, we see that the 2s energy level for the X nl (R VS ) = 0 boundary condition is consistently in the discrete part of the energy spectrum; by contrast, it is unbound up to around τ = 25 eV for the X nl (R VS ) = 0 boundary condition. The eigenvalues in Table II also explain the discontinuities in the P e vs τ and N ub vs τ curves, since they arise when the 2s or 2p level (depending on the boundary condition) transitions from the continuum to the discrete part of the energy spectrum.
The discontinuities observed in the pressure and number of unbound electrons relate to two limitations of the model. Firstly, the fact that we treat the unbound electron density as a constant means the physical problem being solved changes significantly when a new bound level emerges. In Section IV C, we shall explore the impact of treating the unbound electron density in a quantum manner, which should alleviate this problem. Secondly, our definition of "unbound" orbitals -namely those orbitals with energies above the value of the KS potential at the sphere boundary -is an oversimplification. In partially ionized plasmas like the ones we study, there may be some core states which are clearly bound to the nuclei, and likewise some clearly free electron density, but states with energies τ,σ nl ∼ v τ,σ s (R VS ) probably exhibit both bound and free characteristics and therefore cannot be neatly categorised as one or the other. More meaningful definitions for the mean ionization state than a simple energy threshold (which is commonly used in AA models) could make use of quantities such as the electron localization function [198,199], the inverse participation ratio [200,201], or electrical conductivity data [202], but such an analysis is beyond the scope of this paper.
Next, we analyse the electronic pressure as a function of the mass density for fixed temperature τ = 20 eV in Fig. 12. The densities considered range from about 0.01 to 10 times the ambient solid density of Beryllium (indicated by ρ sol in the figure). The results in this plot are strongly reminiscent of what we observed for Hydrogen: namely, the two boundary conditions usually yield pressures within < 20% of each other up to around twice the ambient density (at this temperature); after that, they diverge significantly, with very large negative pressures again observed for the X nl (R VS ) = 0 condition. The reasons for this unphysical behaviour were discussed already for Hydrogen. Comparison between the LSDA and GDSMFB functionals (shown in the lower panel of Fig. 12) indicate that the inclusion of temperature in the XC approximation is most important for Beryllium compressed to around 5 times its ambient density, with this being consistent for both boundary conditions. However, Interestingly, almost the opposite effect is observed for the eigenvalues, shown in Fig. 13. Here, we observe the functional has a small but non-trivial impact on the¯ τ 1s level at densities lower than and including the ambient density. The boundary condition has very little impact on the¯ τ 1s level up to that point, indicating these orbitals do not feel any effect from neighbouring spheres up to that density. The valence energy levels¯ τ 2s and τ 2p separate at much lower densities, demonstrating that the choice of boundary condition is important for predicting ionization energies even for quite diffuse plasmas. The choice of functional has essentially no effect on these energy eigenvalues. Note the different scales on the x-axis, since the orbitals move into the continuum at different densities.

C. Connection with other average-atom models
In this section we make a connection to other existing AA models, by analysing how more sophisticated treatments of unbound electrons and an alternative boundary condition affect properties we discussed so far. We therefore hope to gain an understanding of how our analysis in previous sections might affect the development and usage of AA models more broadly. We use the MUZE code for this comparison.
For the analysis of different treatments of unbound electrons, we consider the transition where ideal refers to the uniform approximation to unbound electrons we used so far, the TF unbound electron density [203] is given by [101,204] n ub (r) = and the quantum unbound density is constructed in a similar way to the bound density (37), by solving explic-itly the radial KS equations (40) for continuum states X l (r) discretized on the energy scale [102,205]. The alternative boundary condition with which we compare does not impose any constraints on the bound KS orbitals at the edge of the Voronoi sphere. Instead, the KS potential is modified as follows, where v τ,M s (r) denotes that the KS potential is modified from its pure form v τ s (r) as defined by Eq. (36). Using the above form for the KS potential, the radial KS orbitals are allowed to "leak out" of the Voronoi sphere and are thus computed up to an infinite radius, since they rapidly decay naturally to zero outside of the sphere. They are still normalized within the sphere according to (45). In MUZE, the unbound orbitals in the quantum treatment satisfy different boundary conditions which amount to continuity of the orbitals at the sphere's edge. Although the above modification to the potential is not strictly a boundary condition in the mathematical solution of the KS differential equations, we henceforth call it the "potential condition" because its role is essentially that of a boundary condition.
For simplicity, we compare using only the LDA functional. Additionally, to make direct comparisons between the MUZE and ORCHID codes more straightforward, we switch from computing pressure with finite differences and instead use the following definition [206]: In Fig. 14, we compare the effect of both the boundary condition used and the treatment of unbound electrons for Hydrogen, on the pressure and MIS. Qualitatively, all results are in quite good agreement for this example, so we also plot the differences relative to the result obtained with ideal unbound electrons and the b.c. X nl (R VS ) = 0, which we call the reference result. Specifically, ∆P e =P ref e −P e , and ∆N ub = N ref ub − N ub , where the superscript ref denotes the reference result. We observe a few common trends from these plots. Firstly, the reference result is largely an upper bound for both the pressure and MIS. Additionally, the MIS results obtained via the potential condition seem to mostly lie in between the two orbital boundary conditions (regardless of the treatment of unbound electrons). Secondly, particularly with increasing temperature, the TF and quantum treatment of unbound electrons yields lower pressures than the ideal results, when compared with the same boundary condition on the potential or the orbital boundary conditions. Finally, especially at higher temperatures (> 10 eV), all approximations yield very similar predictions for the MIS.
Next, in Fig. 15, we plot an analogous set of results for Beryllium, this time with R VS = 4.7 a 0 (ρ m = The second from top panel shows the differences inP e relative to the result obtained from the b.c. X nl (R VS ) = 0, and the bottom panel the equivalent difference for N ub (with differences between data points obtained via linear interpolation).
0.232 g cm −3 ). In this example, the differences between the various approximations are clearer. It is again the case that the boundary condition X nl (R VS ) = 0 seems to consistently predict the highest pressures; also as was observed in Hydrogen, using the ideal approximation for the unbound electron density seems to yield higher pressures than the more advanced TF and quantum models. Of note here is that the choice of boundary condition, particularly for the MIS, seems to have more impact than the choice of treatment of the unbound electrons.
Furthermore, it appears to be the case in this example that the MUZE solutions (using the potential boundary condition) jump from being roughly equal to solutions from our model using the b.c. X nl (R VS ) = 0 to the b.c. X nl (R VS ) = 0 at a temperature around 10 ∼ 15 eV. The comparison of the 2s and 2p energy levels [207] for Beryllium with R VS = 4.7 a 0 in Table III helps explain the MIS results seen in Fig. 15 and also (less directly) Fig. 14. In this table we observe that the eigenvalues given by the potential boundary condition always lie in between those given by the two orbital boundary conditions. This fits with the observation that the X nl (R VS ) = 0 and X nl (R VS ) = 0 boundary conditions are respectively approximate upper and lower bounds for the MIS compared to the results obtained via the potential condition. Also in this table, it appears that changing the treatment of unbound electrons from ideal to quantum has only a small effect on the bound energy eigenvalues.
Based on the above results, the main conclusions that can be drawn are (i) that the ideal approximation has a tendency to overestimate pressure, particularly as the ionization degree increases, and (ii) that the potential boundary condition (60) seems to yield results somewhere between the two orbital conditions (51), (52). In future works, it would be interesting to explore the nature of this relationship at higher densities, when the two ideal quantum  Eventually, evaluating the accuracy of different boundary conditions should be guided by comparisons with more advanced theoretical models, such as KS-DFT-MD models, neutral pseudo-atom AA models that constrain the potential within a correlation rather than Voronoi sphere [122,124], and experimental results.

V. DISCUSSION AND SUMMARY
In this paper, we presented a fully first-principles derivation of a KS-AA model -starting with the fully-interacting, many-body Hamiltonian of electrons and nuclei -and ending with finite-temperature radial KS equations. We methodically considered the underlying assumptions and the interactions that are neglected in this model, yielding insight into the density and temperature limits under which the AA approximation is expected to be accurate. This analysis already yields some ideas regarding future directions for improving AA models: for example, one could go beyond the Born-Oppenheimer approximation and include nonadiabatic effects using the exact factorization method [130,208,209]. Furthermore, through the inclusion (exact or approximate) of higher order terms in the perturbative expansion of the coupling termsŴ J (15), there are possibilities to systematically improve AA models.
In our model, we impose the intuitive criterion that the KS density (which is formally equal to the real electronic density) must be smooth at the Voronoi sphere boundary, Eq. (46). Imposing this criterion leads to the relation (50) for the KS orbitals, which has no unique solution. We considered two options for satisfying this criterion -Eqs. (51) and (52) -in our model, because they are the most simple options and are frequently imposed in AA models.
We observed that these different boundary conditions had a significant impact on results, particularly for higher densities and lower temperatures (as expected), echoing the observation in Ref. 116. This implies that AA models should carefully consider the choice of boundary condition under these limits; even better, further developing the analysis by Rozsynai [101,106] (that these conditions represent band-structure limits) might remove the need to choose one particular condition, and instead one could envisage a scheme (such as the one applied in Refs. 125, 159, and 160) that interpolates between the two. The ultimate goal is to deduce a more accurate boundary condition (or set of boundary conditions) from first-principles, perhaps by considering the effect of the terms neglected during the reduction of the Hamiltonian in Section II A, such as the inhomogenous component H in , or the higher order terms in the perturbative expansion of the average componentĤ av .
We also compared results given by these boundary conditions (51) with an alternative condition on the potential (60), which is used in the MUZE code; there were also significant differences between these results, further emphasizing the importance of boundary conditions in AA models and hence the benefit of further investigation on this subject.
Furthermore, we also investigated the impact of using different approximations for the XC functional. In the case of Hydrogen, we compared with an 'exact' benchmark (F xc [n] = −U [n]). These comparisons indicated that some well-known errors of ground-state DFT, namely the self-interaction and delocalization errors, also affect the prediction of properties in our finitetemperature AA model, particularly the frontier energy levels. To mitigate the impact of these errors, one can borrow from the abundance of solutions suggested in ground-state KS-DFT, such as using hybrid functionals [194,210]. Whilst recognizing that the functional choice is generally much less significant than the choice of boundary condition in our AA model, existing temperature-dependent functionals (such as the GDSMFB functional [73] which we have tested) do not offer serious improvement compared to the standard LSDA and PBE functionals. This motivates the development and use of more advanced temperature-dependent functionals in finite-temperature KS-DFT. Furthermore, the results in this paper offer some support for the observation in Ref. 57, namely that going beyond semi-local approximations may be more important than including explicit temperature-dependence in XC functionals.
To conclude, AA models are a crucial tool in the simulation of materials under WDM conditions. Their computational efficiency not only facilitates calculations over large temperature and density ranges, but also offers an avenue to incorporate advanced features such as non-adiabatic or non-equilibrium effects; such effects are likely to be important in the WDM regime but are too complex to be included in full KS-DFT codes. The firstprinciples derivation and results presented in this paper yield insights regarding potential limitations of KS-DFT AA models: by understanding these limitations, and systematically improving the underlying approximations, there is scope to even further increase the usefulness of AA models.

DATA AVAILABLITY
The data for all the figures in the paper can be downloaded from Ref. [211].

ACKNOWLEDGMENTS
We are grateful to Hardy Gross, Gérard Massacrier and Sang-Kil Son for useful and detailed discussions. T.C. also thanks Zhandos Moldabekov and Nikitas Gidopoulos for helpful comments. Part of A.C.'s initial work on this manuscript was supported by Sandia's Laboratory Directed Research and Development Project No. 200202. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This work was partly funded 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 Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament.

A. DERIVATION OF INTER-CELL COUPLING TERMSŴJ IN MANY-BODY HAMILTONIAN
We wish to expand the terms composingŴ J , W a iJ = 2Z 1 − 2R J · y iJ + |y iJ | 2 in powers of |y iJ |. We recall first the binomial expansion for 1/ √ 1 + , we shall expandŴ J up to second order only in |y iJ |. We henceforth use the notation y i = y iJ for simplicity. The expansions (ignoring higher order terms) forŴ J are thuŝ =Z 2 + 2R J · y i + 3(R J · y i ) 2 − |y i | 2 (67) where we have adopted a form that will be more convenient for expansions in going from Eq. (67) to Eq. (68), using the fact that Ne j=1 = N e . We can expandŴ b ijJ in a very similar manner, We now group terms of the same order in |y iJ | together. We start with the zeroth-order termŴ This term vanishes for charge neutral systems, Z = N e . Next, we consider the first-order termŴ where we used the notationŶ J = Ne i=1 y iJ .
Finally, we consider the second-order termŴ The first term in square brackets vanishes for Z = N e . This completes our derivation of the coupling terms up to second-order inŴ J .

B. CONSTRUCTION OF FREE ENERGY IN OUR AA MODEL
In finite-temperature KS-DFT, the free energy is equal to In our AA model, the unbound electron density is given by the ideal approximation and thus the usual orbitalbased expressions for the KS kinetic energy and entropy cannot be applied. The kinetic energy and entropy are therefore split into bound and unbound components as follows, where the superscripts b and ub denote bound and unbound terms respectively. In our AA model, these components are computed as 1 + e β( −µ σ ) , The remaining terms in the internal energy E[n] take as input the full density (i.e. the sum of the bound and unbound components). In the AA model, these are given by where e xc [n ↑ , n ↓ ](r) is the XC energy density.  A 55, 191 (1997).