1 Introduction↩︎

The advent of high-energy lasers and femto- and even atto-second technologies [1] have made it possible to create and study matter under extreme conditions [2], [3]. Prof. Michael Bonitz and his collaborators have made seminal contributions to the ab initio study of such systems [4]–[6], as reflected in many contributions to this Festschrift.

Material systems under extreme densities \(\bar{\rho}\), and temperatures \(T\), with extreme pressures \(P\) occur naturally in planetary interiors and astrophysical objects. They also occur as transient states that have to be probed on sub-nanosecond timescales [3], [7]–[11] via cutting-edge experiments on high-energy-density materials relevant to experimental astrophysics [8], fusion physics [9], [10] and even for nuclear-stockpile stewardship requirements. The Fermi energies \(E_F\) of compressed matter are large, and although \(T\) (where we use energy units) may be nominally high, \(\theta=T/E_F\) is small. Hence the name “warm-dense matter” (WDM) has been used for such highly-correlated energy-dense matter.

The number of free-electrons per ion, viz., \(\bar{Z}\), interaction potentials \(V_{ab}(r)\), static structure factors \(S_{ab}(k)\) where \(a, b\) = \(e,i\) (i.e., electrons ,ions) as well as their dynamic analogues, e.g., \(S_{ab}(k,\omega)\) [12] are fundamental to the evaluation of the equation of state [10], [13], and transport properties [13], [14] of WDM. They appear in the analysis of experimental results from X-ray Thomson scattering (XRTS) [15]–[19] and X-ray diffraction [7], [11], [20], [21]. The theoretical methods available for the study of WDM hydrogen and the Uniform Electron fluid (UEF). include:

1. \(N\)-atom density-functional calculations, coupled with molecular dynamics (DFT-MD) where \(N\) nuclei and the corresponding number of electrons are treated using Kohn-Sham-Mermin density functional theory theory (DFT). Here the ion configurations are evolved using classical molecular dynamics (MD), This method is referred to as DFT-MD, or more frequently as QMD (quantum MD).

2. Path-Integral-Monte-Carlo (PIMC) methods. Here the partition function is evaluated using Feynman trajectories, and works most efficiently at high-\(T\) where the trajectories are close to classical trajectories.

3. average atom (AA) methods. Here we use a reduction of the \(N\)-nuclei problem to a one-body (i.e., one-ion) problem using an ion-ion XC functional to construct a neutral-pseudo-atom (NPA) model which is an exact DFT procedure, to the extent that the e-e and ion-ion XC functionals are accurate.

4. Methods based on using classical potentials that include quantum effects. Here we use the Classical-Map (CMap) procedure of Dharma-wardana and Perrot to accurately map the quantum electron system at temperature \(T\) to a classical Coulomb fluid at an appropriate temperature \(T_{cf}\).

.

We present comparisons of \(S_{ab}(k)\) and other available quantities from these methods, for the case of the UEF, and for hydrogen plasmas.

2 Review of theoretical methods↩︎

We review the features of the theoretical methods enumerated above in greater detail below. We use atomic units with \(\hbar=m_e=|e|=1\), and \(T\) will usually be given in energy units of eV (1 eV = 11,604K). We also define the electron- and ion- Wigner-Seitz radii \(r_s = [3/(4\pi\bar{n}]^{1/3}, r_{ws} = [3/(4\pi\bar{\rho}]^{1/3}\), where \(\bar{n}\) is the average free-electron density in the plasma, while \(\bar{\rho}\) is the average ion density. The number of free electrons per ion (experimentally measured using a Langmuir probe in low-\(T\) plasmas) is \(\bar{Z}\), with \(\bar{Z} = \bar{n}/\bar{\rho}\). The Fermi wavevector \(k_F = (2\pi^2\bar{n})^{1/3}\) and the Fermi energy \(E_F = k_F^2/2\) are important scales of energy and momentum for the electron subsystem.

2.1 \(N\)-DFT simulations with MD evolution of ion positions↩︎

\(N\)-atom Density functional theory (DFT) uses \(N\) nuclei, located at \(\vec{R}_I, I=1,..,N\) and the corresponding number of electrons. usually with \(N\sim 64-256\). The multi-center electron density \(n(\vec{r};\vec{R}_I)\) is calculated using Kohn-Sham theory in the fixed potential of the ions treated as in a periodic (crystalline) potential. In effect, the Hohenberg-Kohn stationary condition on the one-body electron density leads to the Kohn-Sham equation for a set of non-interacting electrons at the interacting density. \[\label{eeEulerLag46eqn} \frac{\delta F([n(r)];\vec{R}_I)}{\delta n(r)}=0\tag{1}\] Thus the interacting plasma is mapped into a Lorentz plasma in QMD with the ions held fixed. Only the electron many-body problem is simplified in this manner, with the ion-ion manybody problem is left intact. So, the ion configuration \(\vec{R}_I\) is evolved to equilibrium via many static ionic configurations, using classical molecular-dynamics (MD) to obtain a statistical thermal average for the multi-centered one-body density. This procedure is implemented in many well-known computational packages [22], [23], where a choice of \(T=0\) XC-functionals (but currently no finite-\(T\) functionals) is available to reduce the electron many-body problem. The many-ion problem in QMD has to be treated “head-on” via the \(N\)-atom simulation. The MD procedure yields the ion-ion pair distribution function (PDF), viz \(g_{ii}(r)\), and the corresponding \(S_{ii}(k)\), but further computationally costly manipulations are needed to obtain the Kubo-Greenwood response functions required for the determination of the electrical and thermal conductivities \(\sigma\) and \(\kappa\) [24]. Similarly, further manipulations are needed [25] to determine \(S_{ee}(k)\) that enters into XRT spectra  [26]–[28].

The mean number of free electrons per ion, viz., the mean ionization \(\bar{Z}\) is not easily obtained from such \(N\)-atom QMD calculations, leading to several differing estimates of \(\bar{Z}\) [29]. This method, variously referred to as DFT-MD, KS-MD (Kohn-Sham MD) and QMD (quantum molecular-dynamics)  [10] will be referred to here as QMD due to its brevity and wider usage. QMD is an established “workhorse” although it has several major limitations:

1. its accuracy decisively depends on the electronic exchange-correlation (XC) functional. Thus, estimates of the conductivities (\(\sigma,\kappa\)) using QMD may differ by a factor of two depending on the selected XC-functional [14].

2. Given an exact XC-functional, the Hohenberg-Kohn-Mermin theorem of DFT promises to give the total free energy \(F\) (or the ground-state energy \(E\) at \(T=0\)) and the one-body distributions \(n(r),\rho(r)\) of the electrons and ions. But the theory does not provide valid excitation energies or electron eigenfunctions. The Kohn-Sham method maps the interacting many-electron system to a non-interacting system and provides the one-electron energies and eigenfunctions of that fictitious electron system. However, they are still a good approximation in many cases; but the bands gaps etc., of the non-interacting electron system obtained by the Kohn-Sham method need further corrections if they are to represent the actual physical system. For instance, DFT calculations of solid Ge return a metallic solid with no band gap.

Most implementations of the Kohn-Sham method do not correct the Kohn-Sham states for self-interaction effects [30], or for the discontinuities in the XC-functional across band gaps [31]. Thus, when the Kohn-Sham eigenstates and eigenfunctions are directly used for, say, the calculation of the mean ionization \(\bar{Z}\), results with varying degrees of inaccuracy may be obtained [29], [32], [33]. These effects may be substantially corrected if the DFT-MD calculation is augmented with a GW-quasiparticle calculation [34], [35] but the computational costs becomes prohibitive, especially for finite-\(T\) applications.

3. \(N\)-atom DFT coupled with MD (QMD), with \(N\sim 100-500\) at finite-\(T\) is computationally extremely expensive even without GW-quasiparticle corrections because the basis sets needed for the finite-\(T\) calculations increase rapidly for finite \(T\). For instance, the recent analysis of XRTS spectra of highly compressed Be [17]–[19] using QMD would take months of computation, and hence many laboratories choose to remain outside such warm-dense matter studies. Accurate orbital-free DFT methods have been a long-sought goal for speeding up these \(N-atom\) DFT calculations. However, while considerable progress has been made in developing kinetic-energy functionals [36], their accuracies remain well below Kohn-Sham calculations; they also do not provide any one-electron eigenfunctions that are convenient intermediaries in calculating many physical properties of WDM systems.

4. The structure factors \(S_{ab}(k)\) are available only for \(k>k_0\) where \(k_0\) depends inversely on the linear dimension of the simulation box. For instance, Moldabekov et al [25] do not report data for \(k<1.15\) Å\(^-1\) in their DFT-MD studies of hydrogen at 0.33 g/cc.

2.2 Path-integral Monte-Carlo methods↩︎

A quasi-exact method for quantum simulations is afforded by path-integral Monte Carlo (PIMC) simulations, if not for uncontrolled errors arising from the Fermion-sign problem (FSP)  [37]–[40]. Approximate methods for dealing with the FSA have been developed and applications to WDM studies, e.g., [41], [42] and for the calculation of finite-\(T\) XC-functionals [43], [44] are well known. Bonitz and coworkers have used PIMC to provide a data base for the uniform electron fluid and for warm-dense hydrogen [5].

PIMC is even more demanding than QMD for computer resources, and should be considered as a method of providing benchmark results, especially for systems where \(\theta = T/EF > 1\). PIMC becomes inaccurate for low-\(T\) (i.e., \(\theta < 1\)) applications.

We use in this study the PIMC results at \(r_s =\) 2, and 3.32 (where \(r_s\) is the electron Wigner-Seitz parameter) for warm dense hydrogen provided by Moldabekov et al [25] to compare with similar results from an approach based on a classical-map of the quantum electron system where a statistically equivalent classical Coulomb gas is studied.

2.3 Limitations of PIMC and QMD↩︎

Given that both PIMC and \(N\)-atom DFT-MD (i.e., QMD) an computationally very demanding, and may take months of calculation for analyzing current state-of-the-art experiments, accurate first-principles methods that are also computationally very efficient are badly need Such a capability would enable the experimentalist to optimize the diagnostics and settings “on-the-fly", i.e., while the experiment is in progress. For instance, we may consider the prediction of the XRT spectra, or, at the more fundamental level, the prediction of \(S_{ee}(k)\) within computational times of a few seconds as an objective for theory development.

2.4 Average-atom models↩︎

A serious bottleneck in \(N\)-atom DFT is the need to calculate Kohn-Sham states for millions of ionic configurations in resolving the ion-ion many-body problem. The ion-ion many-body problem can also be rigorously reduced to a one-body problem in DFT using an appropriate ion-ion XC functional [45], just as the electron-electron problem is reduced to an effective one-electron problem using an e-e XC functional. In effects, instead of using just Eq. 1 dealing with only electrons, we consider two coupled equations. The first of these is for the electron subsystem with the one-body density \(n(r)\), and the second deals with the ion one-body density \(\rho(r)\). We are considering a uniform system with the mean electron density \(\bar{n}\) and a mean ion density \(\bar{\rho}\). The resulting analysis leads to the construction of a single representative atom around a single nucleus placed in the appropriate plasma medium. It is referred to here as the neutral-pesudo-atom (NPA) model [45]. It is a rigorous DFT model limited mainly by the need to replace e-e and ion-ion many body effects by suitable XC-functionals, and the neglect of certain electron-ion correlations [46]. The latter was the main criticism leveled by Chihara [47] against the NPA model of Dharma-wardana and Perrot [45]. Details of computationally convenient implementations of the NPA have been given in many past publications [13], [48], [49].

While the NPA begins as an all electron calculation, it constructs, internally, many intermediate quantities to simplify the calculation. These involve defining an average ionization \(\bar{Z}\), ion-electron pseudopotentials, T-matrices, as well as ion-ion pair potentials and pair distribution functions. Such pseudopotentials and intermediate quantities are also used in QMD calculations, and they are supplied externally to the QMD calculation as a part of the package. In the NPA, all such intermediate quantities are constructed internally as they are all functionals of the one-body distributions, That is, the NPA uses no external empirical parameters, or elements from intuitive physical or “chemical" models. The inclusion of the XC function for classical ions involves, in some cases, the evaluation of a bridge function using a model hard-sphere fluid, somewhat as the electron XC-functional uses the UEF as the model fluid. Thus, all intermediate quantities other than XC functionals are directly computed from the results of the Kohn-Sham calculations. So, given the material density \(\bar{\rho}\), the nuclear charge \(Z_n\), and the temperature \(T\), the only quantities needed in the NPA from an external theory are the XC-potentials.

In the NPA, \(\bar{Z}\) is evaluated from the finite-\(T\) Friedel sum rule, using the phase shifts of the continuum electrons. These continuum states are not significantly subject to self-interaction errors or XC-discontinuities. Also, the evaluation is made at the minimum of the total free energy \(F([n],[\rho])\) for variations in the one-body densities \(n,\rho\). This means \(\bar{Z}\) corresponds to a full Saha-type calculation inclusive of all the many-body corrections automatically included via XC-functionals. This “automatic” free-energy minimization picture that arises in an NPA calculation is presented in some detail in Ref. [13].

That is, the NPA does not attempt to impose Yukawa- or Debey-Hükel models, Lennard-Jones models, bond-order models etc., taken from outside. However, the ion-electron pseudopotentials obtained from the NPA can be post facto further re-parametrized for convenience and transferabilty, e.g., as modified Heine-Abarankov pseudopotentials [50]. The ion-ion pair-potentials can be represented by Yukawa forms augmented by a Friedel-type oscillatory tail [51]. The neglect of these oscillatory tails prevents the correct prediction of subtle liquid-liquid phase transitions, or Martensitic-type transitions found in many materials.

It has sometimes been claimed that the NPA, or average-atom models cannot correctly simulate systems such as liquid carbon at low temperatures due to the existence of complex covalent bounds in them [52], [53]. We have demonstrated in previous publications that the ion-ion PDFs of complex fluids of such materials like C, Si, obtained using \(N\)-atom DFT-MD can also be obtained in close agreement with the QMD results, even close to the melting line. Since the NPA is based on static DFT, the calculated PDFs, \(S_{ab}(k)\) are long-time averages where transient bonding effects demonstrated in simulations, e.g., in Ref. [52] are of no consequence to thermodynamic and static quantities such as the PDFs of uniform fluids.

An advantage of the NPA over PIMC and QMD is that the structure factors \(S_{ab}(k)\) are available for the whole range of \(k\) values including the \(k\to 0\) limit since the calculation is not limited by the linear dimension \(L\) of the simulation cell used. For instance, a QMD calculation using 64 atoms corresponds to a linear dimension of only 4 atoms, i.e., only a single nearest-neighbour to any “central” atom is treated, and the accessible \(k\) is limited to values larger than \(\pi/L\). Unfortunately, Rayleigh-weight calculations of XRT spectra require results at low scattering wavevectors where many-body effects become more prominent; but such calculations for small-\(k\) are extremely expensive using PIMC or QMD.

2.5 The classical-map approach to WDM calculations↩︎

There is a long history of attempts to use classical “statistical” potentials that mimic quantum effects for use with classical simulations to obtain the properties of quantum systems [54]. The classical-map approach developed by the present author and Perrot [55] falls into this group, although motivated entirely by DFT ideas.

In our discussion of the NPA model, we noted that a fully non-local ion-ion XC-functional needed to reduce the ion-ion many-body problem could be internally generated within the calculation itself. This means that, given an accurate mapping of quantum electrons into a statistically equivalent classical Coulomb fluid, then the fully non-local XC-functionals needed for reducing the electron many-body problem could also be generated internally, within the calculation itself. Furthermore, the whole calculation can be done using only classical statistical mechanics, without recourse to the heavier machinery of quantum statistical mechanics and quantum kinetics. That is, CMap calculations for WDM are even more rapid and economical than NPA calculations.

The current state of the classical map (CMap) approach is that it is quite successful in dealing with the normal uniform electron fluid (UEF) at zero to finite \(T\), at any electron density \(\bar{n}\) and spin polarization \(\zeta\). The method can be justified using density functional theory. Briefly, what has been established [55]–[63] is that the PDFs of a classical fluid (consisting of classical electrons interacting via a modified Coulomb interaction) become equal to that of the quantum UEF at the given density and temperature, if the temperature \(T_{cf}\) of the classical fluid is chosen such that the XC-energy of the classical fluid is equal to the XC- energy of the quantum fluid at the corresponding \(\bar{n},T\) and \(\zeta\). At the time when the classical-map approach was developed, no data for the XC-free energy of the UEF at finite-\(T\) were available. Hence the classical-fluid temperature \(T_{cf}\) at \(T=0\) was evaluated using the \(T=0\) XC-energy. This was named \(T_q\), and the classical-fluid temperature \(T_{cf}\) at finite-\(T\) was estimated using the ansatz \[\label{tcf46eqn} T_{cf}=\surd\left[T^2_q +T^2\right]\tag{2}\] Here \(T_q\) depends only on the \(r_s\) of the electron fluid. Denoting the spin states as ‘\(u,d\)’, the classical Coulomb fluid at \(T_q\) recovers the PDFs \(g_{uu}(r), g_{ud}(r)\), and \(g_{dd}\) at \(T=0\) and for arbitrary spin polarizations \(\zeta\). Here we are concerned with the case \(\zeta=0\), and more details of CMap calculations may be found in Ref. [55], [59], [60].

In effect, since there is no Hartee contribution in uniform fluids, the total ground-state energy is simply the XC-energy. DFT asserts that when two systems modeled to have the same Hamiltonian have the same ground state energy, then their one-body distributions, viz., \(n(r)=\bar{n}g_{ee}(r)\) become equivalent. That is, the \(g_{ee}(r)\) of the classical fluid becomes equivalent to that of the quantum fluid.

While the classical map is successful for mapping the quantum UEF, it has only had limited success with partially ionized hydrogen plasmas [59]–[62] or for more complex systems like aluminum plasmas. First of all, if an ion forms bound states, the classical map needs the electron-ion pseudopotential (say, obtained from an NPA calculation) as an external input as it has no internal theory for constructing bound states and pseudopotentials. Even in the fully ionized case, while the effective classical temperatures for the e-e and i-i interactions are evidently \(T_{cf}\) and \(T\), there is no clear insight for the temperature to be associated with the electron-ion interaction. At this point, such a temperature can be constructed by appealing to the compressibility sum rule [59], or by demanding that \(g_{ei}(r)\) is consistent with what is provided by an NPA calculation.

However, as with PIMC calculations, for sufficiently high \(T\), the CMap method can be used. That is, at high \(T\), the differences among \(T\) and \(T_{cf}\), or what is obtained from the compressibility sum-rule model become small and the results are less sensitive to the choice the effective electron-ion temperature and the method begins to work well.

In this study we present CMap calculations of \(S_{ab}(k)\) for fully ionized hydrogen for comparison with PIMC calculations where available from the work of Bonitz and coworkers. Such comparison are a necessary step for improving the classical-map approach and related less microscopic methods.

3 The degree of ionization \(\bar{Z}\) for Warm dense Hydrogen↩︎

Although a system of \(N\) protons interacting with \(N\) electrons via the Coulomb potential within a given volume may seem a straight-forward problem in quantum statistical mechanics, it turns out to have a very rich phenomenology due to the possibility of formation of a variety of bound states of electrons and ions, phase transitions etc. Initially, when first-principles methods were not available, various phase-transitions involving normal electronic states and even superconducting states had been conjectured [64], [65]. Today, we have the essential tools for studying partially or fully ionized hydrogen plasmas using rigorous first-principles methods. However, even the very concept of “the degree of ionization” [66] has been challenged, leading to many recent discussions of the topic [67], [68].

Bonitz et al [5], in their discussion of the “degree of ionization" of hydrogen plasmas state that”the degree of ionization and the fractions of atoms and molecules are not physical observables. They cannot be rigorously computed, even by a first principles simulation. Any result depends on the used criterion,..., even though the sensitivity to the chosen criterion can be verified and minimized". They concede that “For completeness, we mention that there exist other approaches that use different criteria, in particular, dynamical quantities to estimate the degree of ionization".

More specifically, such computations cannot distinguish and differentiate the species H, (H\(_2)^-\), etc that are formed through bound-electron effects that manifest in partially ionized plasmas. The computed \(g_{\rm HH}(r)\) obtained from such methods, for such partially ionized plasmas may contain a featureless broad peak, or usually an extra “hump” prior to the main peak, with no species differentiation.

Other authors have claimed that there is no quantum operator corresponding to the “degree of ionization” in disputing the use of \(\bar{Z}\) in WDM calculations [69]–[71]. They ignore that temperature is a valid physical property although it has no quantum operator. In the study of WDM, we use quantum statistical mechanics, with a classical heat bath attached to the system rather than a purely quantum mechanical system. Hence, quantities such as the temperature \(T\), the chemical potential \(\mu\), and the mean ionization \(\bar{Z}\) can be introduced as Lagrange multipliers of the theory [45]. They can also be introduced into numerical simulations by introducing thermostats (for \(T\)), numberstats (for \(\mu\)), as well as chargestats (for \(\bar{Z}\)). While the use of thermostats is well developed, the other types of controls for the number of particles or charge neutrality have not been taken up extensively. This is partly why quantities like \(\bar{Z}\) remain currently unacknowledged in first-principles simulations.

Furthermore, as discussed in Ref. [19], there actually exists a quantum operator whose mean value gives us a precise value of \(\bar{Z}\) via the Friedel sumrule that arises naturally in the construction of the neutral-pseudo-atom model of DFT.

The PIMC method, which inputs only electrons and protons does not recognize or distinguish bound and free electrons at any stage of its calculation. Similarly, as we already noted, it fails to identify and distinguish species like H, (H\(_2)^+\), H\(_2\), H\(^-\) that exist in a partially ionized H-plasma. That does not mean that these species have no “physical existence". They all appear within a broadened "bump" in the \(g_{\rm HH}(r)\) calculated by these methods. This failure to recognize the speciation of bound atomic and molecular structures within those methods is no basis for claiming that such species do not exist, In reality, the PIMC, QMD or NPA calculations, all of which begin with an input of electrons and nuclei, should be regarded as a sequence of successive calculations that analyze the speciation in hydrogen and other partially ionized plasmas in an appropriate manner, as the primary plasma, i.e., electrons and nuclei, respond to changes in energy and length scales.

Another well-known feature of partially ionized plasmas is the existence of metastable states connected by liquid-liquid phase transitions. These transitions involve the cooperative action of a large numbers of atoms and distant neighbours. Such transitions exist even in solids (e.g., Martensitic transition) and are beyond the reach of simulations that employ a few hundred nuclei. Such metastable phases may have two different average ionizations at a given density and temperature, and in such cases \(\bar{Z}\) becomes a quick label for indicating such phases. For instance, hydrogen at a given density and temperature may exist as a fully ionized metastable phase with \(\bar{Z}=1\), or with a value less than unity.

Moldabekov et al [25] have studied hydrogen at \(T\)=4.8 eV and density 0.08 g/cc (\(r_s=3.23\)). Hydrogen plasmas within that range of densities is typical of those that show competing metastable states with \(\bar{Z}=1\) and \(\bar{Z}<1\). In fig. 1 we display the behaviour of \(\bar{Z}\) for a hydrogen plasma of density 0.10 g/cc (\(r_s=3)\) over a wide range of temperatures to illustrate the competition between a fully ionized phase and a partially ionized phase.

The region indicated by a box, within the temperature region 2 eV to 6 eV is characterized by two possible plasma phases that we could reach by starting our NPA calculations from either (a) a fully ionized trial solution or (b) from a partially ionized initial starting point, in iterating the Kohn-Sham equations to an equilibrium solution. This behaviour is found in a variety of H-plasmas for \(r_s \sim 3\), possibly including the plasma at a density of 0.08 g/cc studied in Ref.[25] using PIMC.

In Fig. 2, we display \(g_{\rm HH}(r), S_{\rm HH}(k)\) for two H-plasmas of density 0.1 g/cc. One of them is on the fully-ionized branch, while the other is a partially ionized system. The “physical reason’ for the existence of these two possibilities becomes clear from their PDFs. The \(g(r)\) of the fully-ionized system has a sharp peak near \(r\)=0.75 atomic units, and also a broad peak near \(r=3\) atomic units. The convergence in this region is difficult. Although the ion-ion potential is repulsive, by placing a nearest-neighbour proton at \(r < 1\) a.u., the system prevents the formation of bound electrons in atomic or molecular states; but it gains XC-energy in holding a high density of free electrons, with the plasma remaining in the fully ionized state at \(r_s = 3\). The alternative possibility is to allow atomic and quasi-molecular forms of hydrogen to form, gaining energy by reducing the ion-ion electrostatic repulsion, while sacrificing the XC-energy when the free-electron density decreases, with \(r_s\) dropping from 3 to 4.48.

To obtain a PDF on the partially ionized branch at 5 eV itself proved somewhat difficult. However, the results at \(T = 4.2\) eV, \(\bar{Z}\) = 0.3 are presented (curve marked with filled squares) as being typical of this regime. Here, although \(r_{ws} = 3\), \(r_s\) is 4.48 as 70% of the electrons are in bounds states. The interactions among the bound forms (quasi-molecules) are weak, and hence the \(g(r)\) is seen to be rather featureless, although screening is also weak at \(r_s = 4.48\).

As these results are preliminary, further calculations are needed to validate these results. More details of these metastable states, liquid-liquid phase transitions, and their Helmholtz free energies would be presented in a future study. .

4 Classical Map calculations for the UEF and H-plasmas↩︎

Although the classical-map method is not as well grounded as the NPA method, the CMap works well for fully ionized plasmas. These methods require only the electron density Wigner-Seitz parameter \(r_s\), the spin polarization \(zeta\) and the physical temperature as inputs. The spin polarization will be taken as zero, as WDM studies have currently addressed mainly the case \(\zeta=0\). However, the capacity to do spin-polarized calculations is important since many giant planets and other astrophysical objects have strong magnetic fields.

4.1 Calculations using the Classical Map for the UEF↩︎

There have been many calculations of the structure factors \(S_{\zeta, zeta'}(k)\) and PDFs of the uniform electron fluid at \(T=0\), using mainly fixed-node Monte Carlo, and diffusion Monte-Carlo methods. Comparisons of results from the CMap for \(T=0\), both for \(\zeta=0\) and \(\zeta\ne 0\), we given in previous publications [55]. In this study we limit ourselves to more recent finite-\(T\) UEF calculations. In this context we select the following densities \(\bar{n}\) and temperatures \(T\), as PIMC H-plasma calculations are also available for most of them. We use \(\theta=T/E_F\) in the following.

1. \(r_s = 2\), i.e., nominally \(\bar{n}\) = 0.33 g/cc at \(T\) = 4.8 eV.In Ref. [25], \(\theta\) = 1, 0.5,and 0.25 are equated with \(T\) = 12.5 eV, 6.27 eV and 3.13 eV. This is consistent with an assumed ionization of \(\bar{Z}=1\) in hydrogen. NPA calculations at \(r_s=2, T\)=12.5 eV, 6.25 eV, and 3.13 eV established that they can be treated as fully ionized plasmas with \(\bar{Z}=1\). No co-existence of partially ionized states were encountered.

However, Ref. [25] also uses a ‘chemical model’ with \(\bar{Z}\) = 0.72, 0.64 and 0.61. The case \(\bar{Z}\) = 0.72 should correspond to \(\theta\) = 1.26, \(r_s\) = 2.247, while \(\bar{Z}\) = 0.64 implies \(\theta\) =0.683, \(r_s\) = 2.337; and \(\bar{Z}\) = 0.51 implies that \(\theta\) = 0.352, \(r_s\) = 2.375. The authors have not indicated how these \(\bar{Z}\) have been obtained. We will do not examine them any further in this study.

In the case the UEF, we do not have to concern ourselves with \(\bar{Z}\), and we simply use \(r_s=2, \theta=1, 0.5 and 0.25\) for comparisons with the results given in the literature. In the case of H-plasmas, a first-principles determination of \(\bar{Z}\) becomes necessary. Ref. [25] has not reported PIMC calculations for \(\theta=0.5, 0.25\) as PIMC is not viable at low temperatures. They have provided results based on a new ansatz combining time-dependent DFT results for the dynamic structure factor \(See(q,\omega)\) with static DFT results for the density response. We defer the consideration of the cases \(\theta=0.5, 0.25\) to a future study within an improved form of the classical map.

2. \(r_s = 3.23\), i.e., nominally \(\bar{n}\) = 0.08 g/cc at \(T/E_F = 1\), i.e., taken as 4.8 eV in Ref. [25]. Furthermore, Ref. [25] refers to a “chemical model”’ (CM), where a mean ionization of \(\bar{Z}\) = 0.57 is used. If \(\bar{Z} = 1\) is assumed, then setting \(\theta = 1\) is consistent. However, the value of \(\bar{Z}\) for this hydrogen plasma (i.e., 0.08 g/cc at 4.8 eV), obtained using the Friedel sum rule implemented via the NPA code gives \(\bar{Z}\) = 0.365, significantly lower than the chemical models used by Moldabekov et al [25]. If the NPA value of \(\bar{Z}\) is adopted then, \(r_s\) = 4.519, while \(\theta\) = 1.956. Hence we are unable to make a proper comparison with the results given in Fig. 1 of the study by Moldabekov et al [25]. The classical-map approach, being established mainly for fully ionized systems, is also not easily applicable here if the plasma is in a partially ionized state.

However, as seen from Fig. 1 for hydrogen at 0.1 g/cc, a density of 0.08 g/cc at 4.8 eV is a good candidate for being in a fully ionized state and also in a partially ionized state. Unfortunately, even on starting with fully ionized trial inputs, we could not generate Kohn-Sham NPA solutions at 0.08g/cc and 4.8 eV that were fully ionized. As a somewhat close approximation to the density (0.08g/cc) studied by Moldabekob et al., Fig. 3 we display the \(S_{ee}(k)\) for fully ionized hydrogen, and the corresponding UEF, at \(r_s = 3\), \(\bar{\rho}\) = 0.1g/cc and at 5.0 eV.

5 Classical Map and NPA calculations for H-plasma at 0.33 g/cc↩︎

Moldabekov et al [25] provide PIMC calculations for \(S_{ee}(k)\) for hydrogen plasmas at the density 0.33 g/cc, and at \(\theta = 1, 0.5\) and 0.25. We use the nomenclature “uniform electron fluid” (UEF) for quantum electron systems with a uniform non responding neutralizing positive background. The name “uniform electron gas” is also used, most appropriately if \(r_s < 1\). It has been shown [55] with good accuracy that the quantum electron fluid at a physical temperature \(T\) is statistically equivalent to a classical Coulomb fluid at the temperature \(T_{cf}\), Eq. 2 . Now that results for \(F_{xc}(T)\) for the UEF are readily available, the ansatz used in Eq. 2 is not necessary and \(T_{cf}\) can indeed be calculated directly, but choosing \(T_{cf}\) to recover the known \(F_{xc}(r_s,T)\) from a classical-fluid calculation. However, in this study we retain the original ansatz whose accuracy is probably sufficient for astrophysical and WDM applications.

5.1 Results for \(r_s = 2, \theta = 1\)↩︎

Figure 4 is essentially similar to Fig. 4 of Ref. [25], and treats the case \(\theta = 1\), \(\bar{\rho}\) = 0.33 g/cc, at 12.5 eV which corresponds to \(r_s = 2\) since the system is fully ionized.

An approximate form of \(S_{ee}\) proposed by Chihara [72] can also be calculated using the so-called “ion feature" \(I(k)\) of XRTS theory, using the bound electron-density \(n_b(k)\) and the screening density \(n_f(k)\) that is formed around each ion of nuclear charge \(Z_n\). In the present case there are no bound electrons (\(\bar{Z}=1\) and \(n_b)k)=0\). \[\begin{align} \label{SeeChihara46eqn} I(k)=\{n_b(k)+n_f(k)\}^2S_{ii}(k)\\ S_{ee}(k)=S_{ee}(UEF,k)+I(k)/Z_n \end{align}\tag{3}\] Chihara’s approximation does not satisfy the compressibility sumrule for a two component fluid unless the compressibility of the UEF is negligible. We have calculated the relevant quantities, viz., \(n_f(k), S_{ii}(k)\) using the NPA and constructed \(S_{ee}\) according to Eq. 3 . This is displayed (dashed blue line) in Fig. 4. We find that the Chihara approach is only moderately successful here, unlike with more dense plasmas, e.g., highly compressed Be studied in Ref. [19]. In contrast, the results of the CMap approach seem to be in agreement.

6 Conclusion↩︎

Our study of hydrogen plasmas presented here has demonstrated the existence of competing metastable states of ionization in partially ionized hydrogen plasmas, for instance, for \(r_s\sim 3\) and \(0.35 < T/E_F \le 1\). We used a first-principles method, namely, density functional theory where both the e-e and ion-ion many-body problems have been reduced to one body problems. This one-body DFT approach, computationally realized as the NPA model, uses two XC-functionals (one for electrons, and another for ions), where the e-e functional is provided as an external input.

We have also used a classical-map scheme where the electron subsystem is replaced by a statistically equivalent classical Coulomb fluid whose temperature is selected to reproduce the correlation energy of the quantum fluid at \(T=0\). The exchange energy is exactly treated by the selection of the non-interacting PDFs of the finite-\(T\) UEF. This classical-map approach allows a convenient, numerically very fast treatment of fully ionized hydrogen, where all the structure factors (including the electron-electron structure factor) can be computed. The further characterization of competing metastable states of ionization in partially ionized H-plasmas (as revealed in this study) is important, both from the point of view of basic physics, as well as in regard to applications in astrophysics and high-energy density physics. Furthermore, the development of theoretical tools to expose and identify the atomic and quasi-molecular species (such as H, H\(^-_2\), H\(_2\), H\(^-\) etc.) that exist in partially ionized plasmas, is also an important task for future research.

References↩︎