Current issues in finite-$T$ density-functional theory and Warm-Correlated Matter

Finite-temperature DFT has become of topical interest, partly due to the increasing ability to create novel states of warm-correlated matter (WCM). Subclasses of WCM are Warm-dense matter (WDM), ultra-fast matter (UFM), and high-energy density matter (HEDM), containing electyrons (e) and ions (i). Strong e-e, i-i and e-i correlation effects and partial degeneracies are found in these systems where the electron temperature $T_e $ is comparable to the electron Fermi energy. The ion subsystem may be solid, liquid or plasma, with many states of ionization with ionic charge $Z_j$. Quasi-equilibria with the ion temperature $T_i\ne T_e$ are common. The ion subsystem in WCM can no longer be treated as a passive"external potential", as is customary in $T=0$ density functional theory (DFT) dominated by solid-state theory or quantum chemistry. Hohenberg-Kohn-Mermin theory can be used for WCMs if finite-$T$ exchange-correlation (XC) functionals are available. They are functionals of both the one-body electron density $n_e$ and the one-body ion densities $\rho_j$. A method of approximately but accurately mapping the quantum electrons to a classical Coulomb gas enables one to treat electron-ion systems entirely classically at any temperature and arbitrary spin polarization, using exchange-correlation effects calculated {\it in situ}, directly from the pair-distribution functions. This eliminates the need for any XC-functionals, or the use of a Born-Oppenheimer approximation. This classical map has been used to calculate the equation of state of WDM systems, and construct a finite-$T$ XC functional that is found to be in close agreement with recent quantum path-integral simulation data. In this review current developments and concerns in finite-$T$ DFT, especially in the context of non-relativistic warm-dense matter and ultra-fast matter will be presented.


INTRODUCTION
Although there are no systems at zero temperature available to us, it is the quantum mechanics of the simpler T =0 systems that has engaged the attention of theorists. Thermal ensembles usually require the study of extended systems attached to a "heat bath", and within some statistical ensemble. Even perturbation-theory approaches to model systems like the electron gas at finite-T were full of surprises [1].
Condensed matter physics and chemistry could get by with T = 0 quantum mechanics as the input to some sort of thermal theory which is not integrated into the many-body problem. Much of plasma physics and astrophysics could manage with simple extensions of hydrogenic models, Thomas-Fermi theory, extended-Debye theory, and classical 'one-component-plasma' models as long as the accuracy of observations, experiments and theoretical models did not demand anything more from quantum mechanics. On the other hand, at the level of foundations of quantum mechanics, the whole issue of a quantized thermo-field dynamics has been an open problem [2]. Similarly, the theory of 'mixed' systems with classical and quantum components is also a topic of discussion [3]. It is in this context that we need to look at the advent of density-functional theory (DFT) as a great step forward in the quantum many-body problem. The Hohenberg-Kohn theorem published in 1964 was soon followed by its finite-T generalization by Mermin, providing a 'thermal' density-functional theory (th-DFT) in 1965 [4,5], which also saw the advent of Kohn-Sham theory. Hence in 2015 we are celebrating the fiftieth anniversary of both Kohn-Sham theory, and Mermin's extension of Hohenberg-Kohn theory to finite-T .
While DFT provided chemistry and condensed-matter physics an escape from the intractable 'n-electron' problem, in addition to its computational implications, DFT has deep epistemological implications in regard to the foundational ideas of physics. DFT claims that the many-body wavefunction can be dispensed with, and that the physics of a given system can be discussed as a functional of the one-body density. Thus even entanglement can be discussed in terms of density functionals [6,7]. However, it is the computational power of DFT that has been universally exploited in many fields of physics.
The interest in thermonuclear fusion via laser compression and related techniques, and the advent of ultra-fast lasers have created novel states of matter where the electron temperature T e is usually of the order of the Fermi energy E F , under conditions where they are identified as warm dense matter (WDM) [8]. When WDM is created using a fast laser within femto-second time-scales, the photons couple strongly to the electrons which are heated very rapidly to many thousands of degrees, while the ions remain essentially at the initial 'ambient' temperature [9,10]. In addition to highly non-equilibrium systems, this often leads to two-temperature systems with the ion temperature T i = T e , with T e ≫ T i . Alternatively, if shock waves are used to generated a WDM we may have T i > T e . Such ultra-fast matter (UFM) systems can be studied using a fs-probe laser within timescales t such that t ≪ τ ei , where τ ei is the electronion temperature relaxation time [11] of the UFM system. These WCM systems are of interest in astrophysics and planetary science [12], inertial fusion [8], materials ablation and machining [13], in the hot-carrier physics of field-effect transistors and other nano-devices [14].
Early attempts to apply thermal-DFT (also called finite-T DFT, th-DFT) to WDM-like systems were undertaken by the present author and François Perrot in the early 1980s as reviewed in Ref. [17]. This involved a reformulation of the neutral-pseudoatom (NPA) model that had been formulated by Dagens [15] for zero-T problems, as it has the versatility to treat solids, liquids and plasmas.
Originally it was J. M. Ziman [16] (and possibly others) who had proposed the NPA model as an intuitive physical idea in the context of solid-state physics. The electronic structure of matter is regarded as a superposition of charge densities n j ( r − R j ) centered on each nuclear center at R j . In other words, if the total charge density in momentum space was n T (k), then this is considered as being made up of the individual charge distributions n j ( k) put together using the ionic structure factor S( k). This was more explicitly implemented in muffin-tin models of solids, or 'atoms-in molecules' models of chemical bonds that were actively pursued in the 1960s, with the increasing availability of fast computers. The NPA model was formulated rigorously within T = 0 DFT by Dagens who showed that it was capable of the same level of accuracy, at least for 'simple metals' as the LMTO, APW or Korringa-Kohn-Rostoker codes that were becoming available in the 1970s [15]. Wigner's T = 0 exchange-correlation (XC) 'functional' in the local-density approximation (LDA) was used by Dagens. In the finite-T NPA that we have used as our "workhorse", we solve the Kohn-Sham Mermin equation for a single nucleus placed at the center of a large 'correlation sphere' of radius R c which is of the order of 10r ws , where r ws is the Wigner-Seitz radius per ion. Here r ws = {3/(4πρ)} 1/3 whereρ is the ion density given as the number of ions per unit atomic volume. For WDM aluminum at normal compression r ws ≃ 3 a.u. All types of particle correlations induced by the nucleus at the center of the 'correlation sphere' would have died down to bulk-values when r → R c . The ion distribution ρ(r) =ρg ii (r) is approximated as a spherical cavity of radius r ws surrounding the nucleus, and then becoming a uniform positive background [19,20]. This is simpler to implement than the full method implemented in Ref. [21]. The latter involved a self-consistent iteration of the ion density ρ(r) and the electron density n(r) obtained from the Kohn-Sham procedure coupled to a classical integral equation or even molecular dynamics; the simpler NPA procedure is sufficient in most cases.
There have also been several practical formulations of NPA-like models in more recent times. Some of these [22] are extensions of the INFERNO cell-model of Lieberman [23], while others [24] use a mixture of NPA ideas as well as elements of Chihara's "quantum-HNC" models [25]. We have discussed Chihara's model to some extent in Ref. [26]. In true DFT models the electrons are mapped to a non-interacting Kohn-Sham electron gas having the same interacting density but at the noninteracting chemical potential. This feature is absent in INFERNO-like cell-models where the chemical potential is determined via an integration within the ion-sphere or by some such consideration. Thus different physical results may arise (e.g., for the conductivity) depending on how the chemical potential is fixed. Chihara's models use an ion susbsytem and an electron subsystem coupled via a 'quantal Ornstein-Zernike' equation. However, if a one-component electron-gas calculation were attempted via the 'quantal HNC', the known g ee (r) are not recovered. In the two component case, as far as we can ascertain, the S(k → 0) limit is not correctly related to the compressibility.
Thus the Kohn-Sham NPA calculation provides the free-electron charge density pile-up n f (r) around the nucleus. This is sufficient to calculate an electron-ion pseudopotential U ei , and hence an ion-ion pair potential V ii (r) as discussed in, say, Ref. [20]. Once the pair-potential is available, the Hyper-Netted Chain equation (and its modified form incorporating a bridge function) can be used to calculate an accurate g ii (r) if desired, rather than via the direct iterative procedure used in Ref. [21]. This finite-T NPA approach is capable of accurate prediction of phonons (i.e., milli-volt energies) in WDM systems, as shown explicitly by Harbour et al. [27] using comparisons with results reported by Recoules et al. [28] who used the Vienna Ab-initio Simulation Package (VASP).
. Since the XC-functional of DFT is directly connected with the pair-distribution function (PDF), or equivalently with the two-particle density matrix [18], we sought to formulate the many-body problem of ionelectron systems directly in terms of the pair distribution functions g α,β of the system, where α and β count over types of particles (ions and electrons, with two types of electrons with spin up, or down) [19][20][21]. The ionic species may be regard as classical particles without spin as their thermal de Broglie length is in the femto-meter regime at WDM temperatures. This approach led to the formulation of the Classical-map Hyper-Netted-chain (CHNC) method that will be briefly described in sec. .
The attempt to use thermal DFT for actual calculations naturally required an effort towards the development of finite-T XC-functionals [29][30][31][32][33][34][35] [41]) as well as T = 0 DFT-based pseudopotentials are implemented. Currently, these codes also included versions where the single-particle states could be chosen as a Fermi distribution [42] at a given temperature, while they do not include the finite-T XC functionals that are needed for a proper implementation of thermal DFT. These codes are meant to be used at T = 0 or small T since finite-T calculations require a very rapid increase in the basis sets needed for such calculations. It should also be mentioned that Karasiev et al.l [43] have recently implemented finite-T XC within the 'Quantum Espresso' code, as well as given an "orbital-free" implementation, although, as far as we can see, the non-locality problem in the kinetic-energy functional has not been resolved.
However, the availability of DFT-electronic structure codes have opened up the possibility of using them even in the WDM regime. We give several references to such work which contain additional citations to other calculations [28,[44][45][46][47][48]. This renewed interest has re-kindled an interest in the theory of thermal DFT in the context of current concerns [49]. In the following we discuss some of the typical issues that arise in applying thermal-DFT to current problems, as these may range from basic issues to the simple question of "if one can get away with" just using the T = 0 XC functional.
The use of a functional, augmented with gradient approximations etc, is satisfactory as long as the 'external potential' can be considered fixed, as is the usual case in quantum chemistry and solid-state physics. In situations where the external potential arises from a dynamic ion distribution ρ(r), since ρ(r) as well as the electron distribution n(r) depend self-consistently on each other, it is clear that the XC-contribution is a functional of both ρ and n, i.e., the XC-functional is of the form F [n(r), ρ(r)]. Under such circumstances, a direct in situ calculation of the electron g(r) in the presence of the ion distribution has to be carried out, and an 'on-the-fly' coupling constant integration is needed for each self-consistent loop determining n(r) and ρ(r). We presented examples of such calculations for a system of electrons and protons at finite temperatures, in Refs. [50,51], using the classicalmap Hyper-netted Chain technique (CHNC) that enables an easy in situ calculation of the g ee (r), g ei (r) and g ii (r). This approach is at once non-local and hence avoids the need for gradient approximations. Furthermore, the ionion correlations are highly non-local and the LDA or its extensions are totally inadequate since they are described by the HNC approximation.

EXCHANGE-CORRELATION AT FINITE-T
It may be useful to present this section as an 'FAQ' (frequently asked questions) rather than a formal discussion on thermal-XC functionals.

Do we have reliable thermal-XC functionals?
The finite-T XC-functional in the random-phase approximation (RPA) [30][31][32] has been available since 1982, while formulations and parametrizations that go beyond RPA have been available since the late 1980s [33][34][35]. Finite-T XC-data from quantum simulations for the uniform finite-T electron fluid were provided in 2013 by Brown et al. [52], while an analytical fit to their data is found in Karasiev et al. [56]. The XCparametrization of Perrot and Dharma-wardana given in 2000, Ref. [35], was based on a coupling-constant evaluation of the finite-T electron-fluid PDF calculated via the Classical-map Hyper-Netted-Chain (CHNC) [59] method. It closely agrees with the recent quantumsimulation results (Fig 1). Finite-T CHNC-based results are available for the 2D-and 3D-electron gas, as well as other electron-layer systems [61][62][63]. They are in good agreement with path-integral and other Monte Carlo (PIMC) calculations where available.
We consider the data for the 3D system that have been conveniently parametrized by Karasiev et al., (labeled KSDT in Fig. 1). The CHNC f xc (T ) at high temperatures (beyond what is displayed in the figure) show somewhat less correlation than given by PIMC, but correctly approaches the Debye-Hückel limit at high temperatures. In the high-density regime (r s < 1), the RPA-functionals become increasingly accurate as r s → 0. The small-r s regime has also been recently treated by Schoof et al. [55]. It should be stated that when the CHNC mapping was constructed, Franois Perrot and the present author did not attempt to map the r s < 1 regime in detail as it is fairly well treated by RPA methods. Recent simulations by Malone et al. [53] find some differences between their work, and that of Brown at al for r s in the neighbourhood of unity. Similarly, the CHNC data show differences for the r s = 1 curve, as shown in Fig 1. However it is too early to re-examine the small r s regime and review the data of Ref. [53] which are given as the internal energy and not converted to a free energy.
However, it is clear that there is no shortage of reliable finite-T XC-functionals for those who wish to use them.
Can we ignore thermal corrections and use the T = 0 implementations?
While finite-T XC functionals can be easily incorporated into the NPA model or average-atom cell models etc. [23], this is much more difficult in the context of large DFT codes like VASP or ABINIT. Hence the already installed T = 0 XC-functionals have been used as a part of the 'package' for a significant number of calculations for WDM materials, ranging from equation of state (EOS), X-ray Thomson scattering, conductivity etc. Hence the question has been raised as to whether the thermal corrections to the T = 0 XC-functional may be conveniently disregarded.
The push for accurate XC-functionals in quantum chemistry came from the need for 'chemical accuracy' in predicting molecular interactions in the milli-Rydberg range. The current level of accuracy in WDM experiments is nowhere near that. Furthermore, many properties (e.g., the EOS and the total energy) are insensitive to details since total energies are usually very large compared to XC-energies, even at T = 0, unless one is dealing with unusually contrived few-particle systems. However, one can give a number of counter examples which are designed to show that there are many situations where the thermal modification of the T = 0 XC-energy and XC-potential are important.
As a model system we may consider the uniform electron fluid with a density of n electrons per atomic volume, and thus having an electron-sphere radius r s = {3/(4πn)} 1/3 . Since the Fermi momentum k F = 1/(αr s ), where α = (4/9π) 1/3 , the kinetic energy at T = 0 scales as 1/r 2 s , while the Coulomb energy scales as 1/r s . Hence the ratio of the Coulomb-interaction energy to the kinetic energy scales as r s . Thus, the electron-sphere radius r s is also the 'coupling constant' that indicates the deviation of the system from the non-interacting independent particle model. The RPA is valid when r s < 1 for T = 0 systems, for Coulomb fluids. On the other hand, at very high temperatures, the kinetic energy becomes T (or k B T where k B = 1 in our units), while the Coulomb energy is Z 2 /r s , where Z = e = −1 for the electron fluid. Hence the ratio of the Coulomb energy to the kinetic energy, viz., Γ = Z 2 /(r s T ) for Classical Coulomb systems. Here the role of r s is reversed to that at T = 0, and the system behaves as an "ideal gas" for large r s in systems where T ≫ E F . The equivalent of the RPA-theory in the high-T limit is the Debye-Hückel theory which is valid for Γ < 1. A generalized coupling constant that 'switches over' correctly from its T = 0 behavior to the classical-fluid behaviour at high T can be given as: The equivalent kinetic temperature T kin referred to in the above equation can be constructed from the mean kinetic energy as in Eq. (A2) given in the appendix to Ref. [35]. However, the main point here is that there are two noninteracting limits for studying Coulomb fluids. We can start from the T = 0 non-interacting limit and carry out perturbation theory, or coupling-constant integrations to include the effect of the Coulomb interaction λZ 2 /r, with λ moving from 0 to unity (e.g., see Eq. 71 of Ref. [49] for a discussion and references). Alternatively, we can start from the T → ∞ non-interacting limit. This high temperature limit is the 'classical limit' where the system is a non-interacting Boltzmann gas. One can do perturbation theory as well as coupling constant integrations over Γ ′ going from 0 to its required value Γ. The latter approach is well known in the theory of classical fluids. Such results provide standard 'benchmarks' in the context of the classical one-component plasma [34,65], just as the electron gas does for the quantum many-electron problem. However, there is no clear way of evolving from a classical Boltzmann gas at Γ = 0 into a quantum fluid by increasing the Coulomb coupling to its full value, as the anti-symmetry of the underlying wavefunction needs to be included. This problem does not arise if we start from a non-interacting Fermi gas at T = 0. How this problem is solved within a classical scheme is discussed below, in the context of the CHNC method. The 'temperature connection formula' referred to recently by Burke et al. [64] in a thermal-DFT context may be closely related to this discussion.
Although XC effects are important, it is a small fraction of the total energy. They become negligible as T becomes very large, when the total energy itself becomes very large. Thus it is easy to understand that finite-T XC effects are most important, for any given r s , in the WDM range where 0 ≤ T /E F ≤ 1, with E F = 0.5/(αr s ) 2 . Furthermore, in any electron-ion system containing even one bound state, the electron density n(r) becomes large as one approaches the atomic core, and hence there are spatial regions r where T /E F (r) ≤ 1, when finite-T XC comes into play. Since the 'free-electron' states are orthogonal to the core states, the free-electron density pileup n f (r) near a nucleus immersed in a hot-electron fluid is also equally affected, directly and via the core. Furthermore, n f (r) is a property that directly enters into the calculation of the X-Ray Thomson scattering signal as well as the electron-ion pseudopotential U ei (r). Hence the effect of finite-T XC, and the need to include thermal-XC functionals in such calculations can be experimentally ascertained.
In Fig. 2 we present the n f (r) near an Aluminum nucleus in an electron fluid of density 1.81 × 10 23 electrons/cm 3 , i.e., at r s = 2.07 and at T = 10 eV, calculated using the neutral-pseudo-atom method. This temperature corresponds to T /E F ≃ 0.84. Calculations using VASP code for an actual experiment covering this regime has also been given by Plageman et al. [45]. Although the difference in charge densities that arises from the difference between the T = 0 XC and the finite-T XC shown in Fig. 2 may seem small, such charge-density differences translate into significant energy differences as well as into significant X-ray scattering features.
Although Kohn-Sham energies are not to be interpreted as the one-particle excitation energies of the system, they can be regarded as the one-particle energies of the non-interacting electron fluid (at the interacting density) that appears in Kohn-Sham theory. These eigenenergies are also sensitive to whether we use the T = 0 XC-functional, or even to different finite-T functionals. For instance, in Sec. 6 of Ref. [35] we give the Kohn-Sham energy spectrum of warm-dense Aluminum at 15 eV calculated using the PDW-finite-T XC-functional [35], as well as the finite-T Iyatomi-Ichimaru (YI) functional. In summary, the KS-bound states obtained by the two methods (with YI given second) are: at energies (in Rydbergs) of 2115.044 and 2110.199 for the 1s level, 27.86214 and 27.53968 for the 2s level. The outermost level, the 2p-state, has an energy of 25.05646 and 24.81116 from PWD and YI, respectively. Similar proportionate changes are seen in the phase shifts of the continuum states. Thus it is clear that the XC-potentials should have a significant impact, especially in determining the regimes of plasma phase transitions [20,54], finite-T magnetic transitions, as well as in the theory of ionization processes [47] and transport properties.
Another example of the need for finite-T XC functionals is given by Sjostrom and Daligault [57] in their discussion of gradient-corrected thermal functionals. They conclude that "finite-temperature functionals show improvement over zero-temperature functionals, as compared to path-integral Monte Carlo calculations for deuterium equations of state, and perform without computational cost increase compared to zero-temperature functionals and so should be used for finite-temperature calculations". Karasiev et al. [43] have recently implemented the PDW-finite-T XC functional as well as their new fit to the PIMC data in the 'Quantum Espresso' code. They have made calculations of the bandstructure and electrical conductivity of WDM Aluminum. They find that the use of finite-T XC is necessary if significant errors (upto 15% at T /Ef ≃ 0.11 in the case of Al) are to be avoided [58].
Can we define free and bound electrons in an 'unambiguous' manner?
In a 'fully-ionized' plasma all the electrons are in delocalized states. Thus, in stark contrast to quantum chemistry, most of plasma physics deals with continuum processes. WDM systems usually contain some partially occupied bound states as well as continuum states. Thus, if the Hamiltonian is bounded, and if there is no frequency dependent external field acting on the system, there is no difficulty in identifying the bound states and continuum states of the non-interacting electron system used in Kohn-Sham theory. If a strong frequency-dependent external field is acting on the system, the concept of 'bound' electrons as distinct from 'free' electrons becomes much more hazy, and will not be discussed here.
Depending on the nature of the 'external potential', a system at T = 0 may be such that all electrons are in 'bound states'. The latter are usually eigenstates ψ j whose square ψ j (r) become rapidly negligible as r goes beyond a region of localization. The spectrum contains occupied and unoccupied 'bound states' as well as positive-energy states which are not localized within a given region. All states become partially occupied in finite-T systems, and treatments that restrict themselves to a small basis set of functions localized over a finite region of space become too restrictive. Most DFT codes use a simulation cell of linear dimension L with periodic boundary conditions. In such a model the smallest value of k in momentum space is ∼ π/L and this prevents the direct evaluation of various properties (e.g., S(k)) as k → 0. In the NPA model a large sphere of radius R such that all particle correlations have died out is used, and phase shifts of continuum states, taken as plane waves, are calculated. This procedure allows an essentially direct access to k → 0 properties as well as the bound and continuum spectrum of the ion in the plasma. However, the difficulty arises when the electronic boundstates spread beyond the Wigner-Seitz radius of the ion.
The question of determining the number of free electrons per ion, viz.,Z is usually posed in the context of the mean-ionic chargeZ used in metal physics and plasma theory. If the nuclear charge is Z n , and if the total number of bound electrons attributed to that nucleus is n b , then clearlyZ = Z n −n b if the charge distribution n b (r) is fully contained within the Wigner-Seitz sphere of the ion. While n b is well-defined in that sense for many elements under standard conditions, giving, for exampleZ = 3 for Al at normal compression and up to about T =20 eV, this simple picture breaks down for many elements even under normal conditions. If the electronic charge density cannot be accurately represented as a superposition of individual atomic charge densities, the definition of n b becomes more complicated since a bound electron may be shared between two or more neighbouring atoms that form bonds. Transition-metal solids and WDMs have delectron states which extend outside the atomic Wigner-Seitz sphere. Hence assigning them to a particular nuclear center becomes a delicate exercise. However, even in such situations there are meaningful ways to define n b andZ that lead to consistence with experiment. In such situations the proper value ofZ may differ from one physical property to another as the averaging involved in constructing the mean valueZ may change. A similar situation applies to the effective electron mass m * e which deviates from the ideal value of unity (in atomic units), and takes on different values according to whether we are discussing a thermal mass, an optical effective mass, or a band mass that we may use in a Luttinger-Kohn k · p calculation.
Experimentally,Z is a measure of the number of free electrons released per atom. This can be measured from the ω → 0 limit of the optical conductivity σ(ω). Thus, although transition metals like gold have delocalized delectrons, the static conductivity upto about 2 eV is found to indicate thatZ = 1, with the optical mass m * e = 1. Another property which measuresZ is the electronic specific heat. Here again the specific heat evaluated from DFT calculations that use aZ = 1 pseudopotential for Au agrees with experimental data up to 2 eV, while those that use the density of states from all 11 electrons as free-electron states will obtain significantly different answers [66,67] that need to be used with circumspection. That is, such a calculation will be valid only if the d-electrons are fully delocalized and partake in the heating process by being coupled with the pump laser creating the WDM.
The argument thatZ is not a valid concept or a quantum property because there is no 'operator' corresponding to it has no merit. The temperature also does not correspond to the mean value of a quantum operator. In fact, T is a Lagrange multiplier ensuring the constancy of the Hamiltonian within the relevant times scales, whileZ is the Lagrange multiplier that sets the charge neutrality conditionn =Zρ relating the average electron density to the average ion density [21].
Additional discussions regardingZ may be found in Refs. [17,51] and in Ref. [20] where the case of a WDM mixture of ions with different ionization, viz., Al Zj + is treated within a first-principles DFT scheme.

FUTURE CHALLENGES IN FORMULATING FINITE-T XC FUNCTIONALS
In considering a system of ions with a distribution ρ( r) = j δ( r − R j ), and an electron distribution interacting with it, the free energy F has to be regarded as a functional of both ρ(r) and n(r). Hence the ground state has to be determined by a coupled variational problem involving a constrained-search minimization with respect to all physically possible electron charge distributions n(r), and ion distributions ρ(r), subject to the usual formal constraints of n-representability etc. The Euler-Lagrange variational equation from the derivative of F with respect to n(r), for a fixed ρ(r) would yield the usual Kohn-Sham procedure with the rigid electrostatic potential of ρ(r) providing the external potential. However, if no static approximation or Born-Oppenheimer approximation is made, we can obtain another Euler-Lagrange variational equation from the derivative of F with respect to ρ. This coupled pair of equations treated via density-functional theory involves not only the f ee xc , but also f ei xc and f ii c , the latter involving correlations (but no exchange) as it arises from ion-ion interactions beyond the self-consistent-field approximation. In effect, just as the electron many-body problem can be reduced to an effective one-body problem in the Kohn-Sham sense, we can thus reduced the many-ion problem into a "single-ion problem". Such an analysis was given by us long ago [21].
The ion-ion correlations cannot be approximated by any type of local-density approximation, or even with a sophisticated gradient approximation. However, Perrot and the present author were able to show that a fully nonlocal approximation where an ion-ion pair-distribution can be constructed in situ using the HNC equation provides a very satisfactory solution. This is equivalent to positing that the ion-ion correlation functional is made up of the hyper-netted-chain diagrams. However, significant insights are needed in regard to the electron-ion correlation functionals which involve the coupling between a quantum subsystem and a classical subsystem. This is largely an open problem that we have attempted to deal with via the classical-map HNC approach, to be discussed below.
The advent of WDM and ultra-fast matter has thrown out a number of new challenges to the implementation of thermal DFT. A simple but at the moment unsolved problem in UFM may be briefly described as follows. A metallic solid like Al at room temperature (T r ) is subject to a short-pulse laser which heats the conduction electrons to a temperature T e which may be 6 eV. The core electrons (which occupy energy bands deep down in energy and hence not excitable by the laser) remain essentially unperturbed in the core region and at the core temperature, i.e., at T r ≃ 0.026 eV. The temperature relaxation by electron-ion processes is 'slow', i.e., it occurs in pico-second times scales. On the other hand, electronelectron processes are 'fast', and hence one would expect that the conduction-band electrons at T e to undergo exchange as well as Coulomb scattering within femtosecond time scales, consistent with electron-electron interactions timescales. Thus, while we have a quasiequilibrium of a two-temperature system holding for up to pico-second timescales, the question arises if one can meaningfully calculate an exchange and correlation potential between the bound electrons in the core at the temperature T r , and the conduction-band electrons at T e , with T e ≫ T i . While we believe on physical grounds that a thermal DFT is applicable at least in an approximate sense, an unambiguous method for calculating the two-temperature XC-energies and potentials is as yet unavailable.

Classical-map Hyper-Netted Chain Method
Once the pair-distribution function of a classical or quantum Coulomb system is known, all the thermodynamic functions of the system can be calculated from g(r). The XC-information is also in the g(r). Only the ground-state correlations are needed in calculating the linear transport properties of the system. Hence, most properties of the system become available. It is well known that correlations among classical charges (i.e, ions) can be treated with good accuracy via the the hyper-netted-chain equation, but dealing with the quantum equivalent of hyper-netted-chain diagrams for quantum systems is difficult, even at T = 0 [68].
When we have an electron subsystem interacting with the ion subsystem, obtaining the PDFs becomes a difficult quantum problem even via more standard methods. We need to solve for a many-particle wavefunction which rapidly becomes intractable as the number of electrons is increased beyond a small number. The message of DFT is that the many-body wavefunction is not needed, and that the one-particle charge distribution n(r) is sufficient. While the charge distribution at T = 0 involves a sum over the squares of the occupied Kohn-Sham wavefunctions, at very high T the classical charge distribution is given by a Boltzmann distribution containing an effective potential felt by a single 'field' particle, and characterized by the temperature which is directly proportional to the classical kinetic energy.
In CHNC we attempt to replace the quantum-electron problem by a classical Coulomb problem where we can use a simple method like the ordinary HNC equation to directly obtain the needed PDFs, at some effective 'classical fluid' temperature T cf having the same density distribution as the quantum fluid. The electron PDF g 0 (r) of the non-interacting quantum electron fluid is known at any temperature and embodies the effect of quantum statistics (Pauli principle). Hence we can ask for the effective potential βV P au (r) which, when used in the HNC, gives us the g 0 (r), an idea dating back to a publication by F. Lado [69]. This ensures that the noninteracting density has the required "n-representable" form of a Slater determinant. Of course, only the product P (r) = βV P au (r) can be determined by this method, and it exists even at T = 0. Then the total pair potential to be used in the equivalent classical fluid is taken as βφ(r) = P (r) + βV Cou (r). How does one choose β = 1/T cf since the Pauli term is independent of it?
To a very good approximation, if T cf is chosen such that the classical fluid has the same Coulomb correlation energy E c as the quantum electron fluid, then it is found that the PDF of the classical Coulomb fluid is a very close approximation to the PDF of the quantum electron fluid at T = 0. There is of course no mathematical proof of this. However, from DFT we know that only the 'correct'  [71] and HNC. In (c) the paramagnetic g(r) at rs=1 and rs=10, T=0 are compared with DMC. (d)Finite temperature PDFs (CHNC) for T /EF =2, rs=5 would correspond to a WDM at ≃ 3.6 eV (∼ 42,000 K).
ground state distribution will give us the correct energy, and perhaps it is not surprising that this choice is found to work. The T cf that works for the T = 0 quantum electron gas is called the "quantum temperature" T q . More details of the method are given in Ref. [35]. There it is argued that, to a good approximation, for a finite-T electron gas at the physical temperature T , the effective classical fluid temperature T cf = T 2 q + T 2 . This has been confirmed independently by Datta and Dufty [70] in their study of classical approximations to the quantum electron fluid. Thus CHNC provides all the tools necessary for implementing a classical HNC calculation of the PDFs of the quantum electron gas at finite-T .
We display in Fig. 3 pair-distribution functions calculated using CHNC, and those available in the literature from quantum simulations at T = 0, as finite-T PDFs from quantum simulations are hard to find. In any case the classical map is expected to be better as T increases and the T = 0 comparison is important. In the figure, diffusion Monte Carlo (DMC) and variational Monte Carlo (VMC) data [71] are compared with CHNC results. In fig. 3 the parallel-spin PDF is marked g 11 (r), while the anti-parallel spin PDF is marked g 12 (r). The latter has a finite value as r → 0 as there is no Pauli exclusion principle operating on them. Furthermore, the the mean value of the operator of the Coulomb potential, i.e., e 2 /r, is of the form {1 − exp(−k dB r)}/r, where k dB is the thermal de Broglie wavelength of the electron pair, as discussed in Ref. [59]. This 'quantum-diffraction' correction ensures that g 12 (r → 0) has a finite value, as seen in the figure. It is in good agreement with Quantum Monte Carlo results. Thus the CHNC is capable of providing a good interpretation of the physics underlying the results of quantum simulations. Needless to say, unlike Quantum Monte Carlo or Path-Integral simulation methods, the CHNC integral equations can be implemented on a laptop and the computational times are imperceptible.
Using the PDFs g(r, T, λ) calculated with a scaled Coulomb potential λV Cou (r), a coupling constant integration over λ can be carried out to obtain the XC-free energy F xc (r s , T ) as described in detail in Ref. [35]. As seen from Fig. 1, this procedure leads to good agreement with the thermal-XC results from the PIMC method, while only the T = 0 spin-polarized E c data were used in constructing T q . Furthermore, since T cf tends to the physical temperature at high T , and since the HNC provides an excellent approximation to the PDFs of the high-T electron system, the method naturally recovers the high-T limit of the classical one-component plasma. Note that we could NOT have started from the high-T limit of an ideal classical gas and used the well-known classical coupling constant (i.e., Γ integration method, e.g., see Baus and Hansen or Ichimaru [34,65]) to determine f xc from an integration that ranges from Γ = 0, T = ∞ to the needed temperature (i.e, the needed Γ). This is because there is no clear method of capturing the physics contained in T q , and ensuring that Fermi statistics are obeyed (e.g., via the introduction of a βV P au (r)), as there is only Boltzmann statistics at Γ = 0.
The ability of the CHNC to correctly capture the thermal-DFT properties of the finite-T quantum fluid suggests its use for electron-ion systems like compressed hydrogen (electron-proton gas), or complex plasmas with many different classical ions interacting with electrons [50], without having to solve the Kohn-Sham equations as in Bredow et al. [72]. The extensive calculations of Bredow et al. establish the ease and rapidity provided by CHNC, without sacrificing accuracy. CHNC has potential applications for electron-positron systems or electron-hole systems where both quantum components can be treated via the classical map. It also provides a partial solution to the still unresolved problem of formulating a fully-nonlocal 'orbital-free' approach that directly exploits the Hohenberg-Kohn-Mermin theory, without the need to go via the Kohn-Sham orbital formulation. CONCLUSION We have argued that our current knowledge of the thermal XC-functionals is satisfactory and the stage is set for their implementation in practical DFT codes. Noting the complexity of warm-dense matter, we have emphasized simplifications as well as extensions which do not sacrifice accuracy. In this respect the neutral-pseudo atom model can, in most cases do the work of the ab initio codes like VASP, and handle high-temperature problems that are beyond their scope. Orbital-free approaches [43] will also become increasingly useful, especially at intermediate and high T /E F . Nevertheless, the ab initio codes are needed at low-temperature low-density situations involving molecular formation, where the NPA breaks down as it is a "single-center" approach. However, in many WDM cases, we need to go beyond the picture where the ion subsystem is held static, and the electrons only feel them as an 'external potential'. Hence we have emphasized the need for calculating not just the XC-functionals for electrons, but also the classical correlation functionals for ions, as well as the ion-electron correlations directly, in situ, via direct coupling-constant integrations of all the pair-distribution functions of the system, ensuring a fully non-local formulation where gradient expansions are not needed. In fact, there is no need for any XCfunctionals in such a scheme. To do this efficiently and accurately, we have proposed a classical map of the quantum electrons and implemented it in the CHNC scheme which depends on DFT ideas. This capacity is not found in any of the currently available methods. CHNC has been used to construct a finite-T XC functional for electrons more than a decade before PIMC results became available, and it turns out that the CHNC results are accurate. The CHNC scheme has been successfully used for calculating the equation of state and other properties of warm dense matter as well as multi-component T = 0 electron-layer systems, thick layers etc., that are expensive to treat by quantum simulation methods, but relevant for nanostructure physics.