Reconciling ionization energies and band gaps of warm dense matter derived with ab initio simulations and average atom models

Average atom (AA) models allow one to efficiently compute electronic and optical properties of materials over a wide range of conditions and are often employed to interpret experimental data. However, at high pressure, predictions from AA models have been shown to disagree with results from ab initio computer simulations. Here we reconcile these deviations by developing an innovative type of AA model, AVION, that computes the electronic eigenstates with novel boundary conditions within the ion sphere. Bound and free states are derived consistently. We drop the common AA image that the free-particle spectrum starts at the potential threshold, which we found to be incompatible with ab initio calculations. We perform ab initio simulations of crystalline and liquid carbon and aluminum over a wide range of densities and show that the computed band structure is in very good agreement with predictions from AVION.


I. INTRODUCTION
The electronic and optical properties of materials change drastically with increasing pressure or density, which affects their equation of state, opacity, and transport properties. These changes have profound implications for the interior structure and evolution of stars and planets [1,2]. The modification of electronic states with temperature and pressure is a very active research subject [3,4]. The characterization of electronic properties is especially challenging in the regime of warm dense matter (WDM) because pressure and temperature are high, particles are strongly interacting, and the system is partially degenerate. One relevant phase transition is the dissociation and metallization of hydrogen [5][6][7][8] that changes from an insulating, molecular substance to a monoatomic, electrically conducting fluid. This phase change is a manifestation of a Mott transition [9,10].
Due to the development of high-power lasers and other facilities, the capabilities to create and diagnose WDM have increased substantially and a variety of new experimental results have been obtained [11][12][13][14][15][16][17]. Many characterize modifications of the atomic structure due to the plasma environment. A key quantity is the degree of ionization, which may be enhanced by the continuum lowering phenomenon, or ionization potential depression (IPD), that describe the reduction of the ionization energy compared to the value of an isolated ion. X-ray Thompson scattering measurements (XRTS) provide direct access to the IPD [18][19][20]. The analysis of the scattering signal is based on the Chihara formalism [21,22], which relies on the separation of bound and free electronic states as well as well-defined direct transitions between the two. There are experimental results both supporting and contradicting the two classic IPD models of Stewart-Pyatt (SP) [23] and Ecker-Kröll (EK) [24]. SP models are widely used to predict plasma properties like ionization balance, equation of state, opacities, electrical and thermal conductivities.
There are a variety of theoretical approaches to study WDM and the IPD [25][26][27]. With ab initio simulations, one follows the evolution of many ions in a periodic sim-ulation cell by coupling density functional theory (DFT) [28,29] to molecular dynamics (MD) [30]. In comparison, average atom (AA) models are computationally far less demanding because they place only a single nucleus at the center of a cloud of electrons and determine the resulting states. Many variations of the AA approach exist [25,[31][32][33]. However, when these methods are applied to guide or interpret experiments at extreme densities or highly degenerate plasma conditions, the results appear to deviate from predictions of ab initio methods [26,[34][35][36], which provided the initial motivation for this study.
In this paper, we use DFT calculations to study the evolution of electronic states of dense C and Al over three orders of magnitude in density and then reproduce them with a novel, general-purpose type of AA model. For the latter, we adopt an alternate view of the continuum. Usually, AA models embed the ion in some kind of jellium and the continuum of states starts at the threshold of the potential. We determine that it is this property that leads to a disagreement with ab initio results. As a matter of fact, valence and conduction bands as derived from DFT-MD simulations exhibit no clear link to such a potential threshold.

II. AB INITIO SIMULATIONS
We performed ab initio calculations for dense carbon and aluminum with the Abinit code [37][38][39][40][41][42]. For both materials, we performed static DFT calculation using perfect lattices (T = 0 eV) and DFT-MD simulations at a finite temperature (T = 12.5 eV). For aluminum, we employed the local density approximation (LDA) to incorporate exchange-correlation (XC) effects [43] and the PBE-GGA functional was used for carbon [44]. For both materials, we generated new pseudopotentials with very small core radii using the Opium code [45]. The 1s core electrons were frozen in the norm conserving pseudopotential for aluminium. For carbon, we even performed all-electron calculations.
As the predicted crystal structures of aluminium and carbon change with density, we change the structure at T = 0 according to Ref. [46] for carbon and Ref. [47] for aluminium. For carbon, the order of the phases is diamond (3.52 to 7.66 g cm −3 ), BC8 (7.87 g cm −3 ), simple cubic (sc, 15.9 g cm −3 ), simple hexagonal (sh, 23.29 g cm −3 ) and face-centered cubic (fcc, 46.36 to 640.83 g cm −3 ). A 12 3 k-point grid with a plane wave cutoff energy of 600 Hartrees (Ha) was used for diamond. The BC8, sc, and sh calculations used a k-point grid of 32 3 points and a cutoff energy of 700 Ha. The static T = 0 calculations for aluminium were carried out for both the fcc and body-centered cubic (bcc) [48] phases to get a better understanding on how the phase affects the valence band gap. For the fcc unit cell a 64 3 k-point grid with 12 bands and a plane wave cutoff energy of 400 Ha has been used. The bcc calculations used a 12 3 k-point grid with 20 bands and a plane wave cutoff at 400 Ha. The band gap reported for these solid phases is the indirect gap over the sampled Brillouin zone. The MD simulations at T = 12.5 eV were run using a Nosé-Hoover thermostat and with a time step of δt = 0.09 fs and with 8 atoms in the periodic supercell which was sampled only at the Γ-point. The number of bands considered was 450 and the plane wave cutoff was 500 Ha for carbon. The aluminum DFT-MD simulations considered 630 to 640 bands and a cutoff energy of 400 Ha. Electronic density of states (DOS) and band gap results at finite temperatures were derived by averaging over the DOS of 10 sample configurations obtained via DFT-MD and which where postprocessed with a static run using 8 3 k-points. These DOS were averaged after aligning them first at the chemical potential µ (a few eV below the Fermi energy at these low temperatures). Figure 1 shows the DOS, weighted by the Fermi-Dirac factor 1/(1 + e (ε−µ)/kT ), at T = 12.5 eV and for densities a few times the solid density. For illustration purposes the DOS were smoothed using a Gaussian with a HWHM of 0.8 eV. For carbon, the peaks below −250 eV repre-sent the 1s states while the 2s and 2p states have already merged with the continuum at these conditions. With increasing density, we find that the continuum edge shifts towards lower energies when compared with the Fermi energy, consistent with earlier findings [49]. For the following discussion, we define the valence band gap (VBG) of carbon to be the gap between the 1s states and the bottom of the continuum band.
For aluminum at the lowest density of 8.03 g cm −3 , the 2s and 2p states form two well-separated peaks. At 22.0 g cm −3 , these peaks have broadened. At yet higher densities, they will merge (not shown). The VBG of aluminum is defined to be the gap between the uppermost 2s/2p state and the continuum band.
Our simulations show consistently that the continuum edge shifts towards lower energies, well below the Fermi energy, as density increases. Eventually gaps between valence and continuum bands close. In traditional AA models, the bottom of the continuum begins at the maximum of the potential and a state is pressure-ionized as soon as it crosses this threshold. This leads to inconsistent results between the DFT and AA approaches because in DFT, one derives the free-particle spectrum without introducing such a threshold. In order to resolve this discrepancy and to reproduce the fully self-consistent many-body DFT calculations, we developed a novel type of AA approach, AvIon, that is much more efficient than DFT-MD and can be readily applied to arbitrary materials and thermodynamic conditions. As we will show below, AvIon reproduces the DFT-MD predictions very well. It also explains apparent discrepancies between the usual AA models and DFT-MD simulations.

III. THE AVION MODEL
In ab initio many-body simulations, the band structure is obtained from eigensolutions of the Schrödinger equation in a periodic box for an effective potential resulting from a collection of electrons and nuclei, while taking advantage of the Kohn-Sham scheme [50]. In our AvIon approach, we intend to solve the same equations to obtain the band structure and wavefunctions around a single nucleus in a limited volume. The only connection with the plasma environment arises from the boundary conditions of the wavefunctions at the surface of that volume. In an actual many particle system these boundary conditions are highly complex and vary from ion to ion depending on the environment. Still these variations can be represented well enough with approximate methods as long as the predictions are validated through comparisons with many-body simulations.
AvIon is a spherical model with a neutrality radius, R. To obtain agreement with many-body simulations, we introduce novel boundary conditions at R that allows us to derive bound and free states within a common scheme. As we will show, this has deep implications for the band structure of the continuum states and for the magnitude of band gaps. Conversely, traditional AA models set the boundary conditions at infinity. This leads to a clear distinction between a bound (discrete) and a free (continuous) spectrum, which is not compatible with predictions of many-body simulations at high density.
For an element with nuclear charge, Z, and at given temperature, T , the AvIon approach self-consistently determines the electronic density profile, n e (r), and an effective potential, V (r), within the Kohn-Sham scheme [50]. The spherical potential is the sum of the Coulombic nuclear attraction, the Hartree potential v H resulting from the electronic density through Poisson equation, and an exchange-correlation part in the LDA approximation. (We use [51][52][53].) The potential v H inside R can be written, up to a constant, as a function of the sole charge inside this same volume: where Q(r) = r 0 4πr 2 n e (r ) dr is the number of electrons inside radius r. We do not need to assume any definite value for V , and leave v H (0) unspecified.
The electronic density (including all electrons) is in turn determined from the spherical eigenstates φ ε m = 1 r P ε (r)Y m (r) in the potential V . Their radial parts obey (with m the electron mass, the Planck constant), are normalized to unity inside the ion sphere, and populated according to a Fermi-Dirac distribution. Accordingly, where the g (ε) are partial DOS for the angular momentum . The chemical potential µ is adjusted so that Z electrons exactly pile up inside of radius R, i.e. Q(R) = Z. From electro-neutrality one has dV dr (R) = 0 (neglecting v xc ). The value V (R) defines the potential threshold on the energy scale.
The distinctive feature of AvIon is in the way the DOS is constructed. The continuum in traditional AA models starts above that threshold (setting V (r ≥ R) = 0), and includes all ε > V (R) solutions, matching them to a combination of Bessel functions (free solutions) at the R boundary. We show in this paper that the band structure inherent to condensed matter, which persists for crystals under WDM conditions as shown recently by Bekx et al. [54] (see in particular their Fig. 7), or in DFT simulations, challenges this simple view, yet can be recovered with a modified AA model.
For bound states, the concept of bands instead of isolated states (which give a δ-function contribution to some g (ε)) has already been used in a number of AA models [55][56][57][58]. Regarding states above threshold, some preliminary works have explored the modification of the DOS due to neighbors using multiple scattering theory [59][60][61]. While this is in the spirit of our work, our aim is to avoid the huge complexity it adds.
In AvIon for every angular momentum, , the partial DOS g (ε) = k g k (ε) is a sum over successive energy bands k. For given , all wavefunctions within the k-th band have k − 1 nodes in the interval ]0, R[. In addition, the lower ε − k and upper ε + k energy limits of this k-th band result from the eigensolutions of the radial Schrödinger equation (Eq. (3)) with the respective boundary conditions d(P ε /r) dr (R) = 0 and P ε (R) = 0. With these two conditions, we intend to mimic the bonding and antibonding character of typical molecular wavefunctions. Note that the P ε (R) = 0 condition for the upper band limits does not imply that there is a hard wall at radius R, but that the total (over the whole plasma) wavefunction, a part of which we solve for in the ion sphere, has a node at R. We emphasize that we no longer make a distinction between bound (below threshold) and free (above threshold) states, as is done in traditional AA models. Still we need to make an assumption for the DOS inside each band.
Between the respective limits ε ± k , we choose in this work the Hubbard functional form g k (ε) = 2 . This is undoubtedly oversimplified. The DOS of the bands is actually controlled by a complex interplay between the ion and its neighbours, and its computation requires involved techniques. For instance the recent work of Starrett & Shaffer [61] uses multiple scattering approach, but requires around one hour of CPU time. Our goal is precisely to bypass this computationally highly demanding step: AvIon necessitates a few seconds per point [63]. With this DOS at hand and wavefunctions computed for a set of energies inside the bands, the electronic density can be obtained from Eq. (4).
As a summary, the three parameters (Z, T, R) entirely determine an AvIon calculation. The iterative procedure is initiated with a trial electronic density as input. The potential is determined through Eqs. (1) and (2). Bands and wavefunctions are obtained from the solutions of Eq. (3). The resulting electronic density is constructed from Eq.(4) with µ adjusted for electroneutrality. The result is compared with the input. If different, a new loop is initiated with a mix of the two densities, and the procedure is continued until full convergence is achieved.
To compare with results at given ion density n i (or mass density ρ = Am u n i , where A is the atomic mass of the element), we set in this paper R equal to the Wigner-Seitz radius, R WS = (3/(4πn i )) 1/3 . Note that this is only one plausible choice among other possibilities that we don't explore in this work [64].
As an illustration of our approach, Fig. 2 shows the evolution of band energies over a large density range for carbon and aluminum at temperatures T = 12.5 eV and T = 0.1 eV. The reference energy in panel (a) and (b) is the potential threshold V (R). In order to better visualize the VBG, energies are counted from the top of the 1s band for carbon in (c), and the top of the 2p band for aluminum in (d). At low density, bands below the threshold appear as sharp, localized bound states. With increasing density, they widen as they approach the potential threshold, which they eventually cross in a continuous way. Above a density of a few g cm −3 the bottom of the continuum is associated for carbon first with the 2s and later with the 2p band, while for Al, it is first the 3s and later the 3d band. The void region between the top of the valence band (1s for C, 2p for Al) and the continuum defines the VBG, as shown with arrows in the panels.
Our AvIon model does not distinguish bound and free states. States rather evolve from being localized to delocalized. To compare with traditional AA and IPD models, we introduce a pseudo-ionization potential (PIP) in AvIon, which we define to be the difference between a state's energy and the threshold potential, V (R), even though the PIP is not linked to any physical quantity in our approach. We show the PIP for the 1s state of C in panels (a) and (c) of Fig. 2, and for the 2p state of Al in (b) and (d). Both decrease rapidly with density. The VBG predicted by our AvIon approach remains much larger because the region void of states extends well above the potential threshold. This explains discrepancies reported in Ref. [34] between very small gaps in AA models and large VBGs in DFT-MD results, as is demonstrated by further analysis in the following section.
In AvIon, we can predict the K-edge by subtracting the 1s state's energy from the chemical potential (or Fermi energy), as is illustrated for carbon in Fig. 2(c). In the same panel, we have included as blue dots the values computed from DFT simulations in Ref. [26] for T = 1.3 eV. The agreement with the K-edge inferred from our AvIon approach is excellent. The increase of the K-edge energy emerges as a 'competition between continuum lowering and Fermi surface rising' [26,35,36] when including the potential threshold as an intermediate step in the analysis. Still in AvIon and in DFT, the K-edge energy can be determined directly.

IV. COMPARISON BETWEEN METHODS
The DOS as obtained from AvIon are compared in Fig. 1 to DFT-MD results at T = 12.5 eV and for the lowest densities. To ease comparison, the DOS have been convolved using a Gaussian distribution for the neutrality radii R with a HWHM equal to R WS /40, which slightly broadens the valence peaks. As AvIon is based on a quite simple prescription for the DOS shape inside a band, a perfect match is not expected. Nevertheless the model reproduces remarkably the DFT-MD results and their trends, especially since it requires a few seconds against hundreds of hours for DFT. The continuum edges and the valence band positions agree quite well for carbon, The position of the potential thresholds V (R) are also drawn as vertical bars in the top of the panels. For carbon, the continuum starts just below the threshold at 3.51 g cm −3 , and just above at 21.2 g cm −3 . For aluminum, it is just below the threshold at 8.03 g cm −3 . These situations would be hard to distinguish experimentally from the usual equality in traditional AA models. In contrast, for aluminum at 22.0 g cm −3 , the continuum edge starts some 40 eV above the threshold. These re-sults are easily visualized through inspection of the upper panels of Fig. 2. In Fig. 3, we compare results of DFT, AvIon, and SP methods over three orders of magnitude in density. For both T = 0 conditions and fluid states at T = 12.5 eV, we show results from DFT(-MD) simulations and AvIon calculations. The latter use T = 0.1 eV instead of T = 0 to avoid numerical issues, but as a matter of fact the AvIon results are fairly similar in these low temperature range except for a slight difference for Al at low density.
For carbon, AvIon reproduces the DFT band gaps very well but a small deviation remains for the highest densities. AvIon predicts a gap closure at a density of ≈ 500 g cm −3 while the gap in static DFT calcula- For aluminum, the agreement is very good up to a density of around 20 g cm −3 . Then the AvIon band gap increases to higher values than is predicted with DFT methods. However, the following decrease and eventual band gap closure at ≈ 190 g cm −3 occurs in a range compatible with DFT predictions. The difference in the latter may be traced back to the different crystal or fluid state in the DFT simulations, while the environment in AvIon is always with spherical symmetry. Note that AvIon shows a marked change of the band gap slope around 32 g cm −3 , which is linked to the switch from the 3s to the 3d band used to define the continuum edge (see Fig. 2(d)). This change in slope is also seen in our DFT simulations. A similar slope change, though less dramatic, occurs for carbon at around 35 g cm −3 due to the switch of the lowest continuum band from 2s to 2p.
In Fig. 3, we also plot the PIP computed with AvIon for T = 12.5 eV (it is almost the same for T = 0.1 eV on this density range.) Its value goes to zero at already ≈ 170 g cm −3 for C and ≈ 30 g cm −3 for Al. From the point of view of traditional AA models, this would imply that the involved bands have already merged with the continuum at these conditions. However, both our AvIon and DFT results predict there is still a sizable VBG present under these conditions.
It is useful to compare when AvIon and SP models predict a given state to become pressure ionized. From around 1 g cm −3 up to the density where the valence band crosses the threshold, carbon has two localized electrons and aluminum has ten. We hence plot in Fig. 3 the modified IP, E 0 − ∆ SP , for the 1s state of C 4+ and C 5+ , and for the 2p state of Al 3+ and Al 4+ , where E 0 is the isolated ion spectroscopic value and ∆ SP the IPD predicted by the SP model. It shows little resemblance to the VBG curves computed with AvIon and DFT. SP models predict pressure ionization to occur at lower densities, at values in between the point where AvIon's PIP goes to zero and AvIon and DFT predict the VBG to close.
It should be noted that AA models and DFT simulations involve approximations which do not allow reproduction of the measured energies of isolated ions, which are input parameters in SP models. By shifting SP curves for C 4+ and Al 3+ , i.e. adjusting E 0 , it is possible to match the evolution of AvIon's PIP over most of the density range (see lines labelled "shifted" in Fig.3.) In these degenerate conditions, AvIon and SP models make similar predictions for screening effects and ∆ SP reduces to the ion sphere expression, 3ze 2 /(2R), where z = 4 for C and z = 3 for Al. The residual discrepancy for the densities where both curves go to zero may be explained by noticing that the PIP is counted from the top of the valence band; the SP model rather refers to an eigenstate that would be computed with the condition P ε (∞) = 0. Hence it is situated in the band below its top and is ionized at a higher density.

V. CONCLUSION
In this paper, we have developed a novel type of AA approach, AvIon, to obtain consistency with much more computationally demanding many-body DFT simulations. We studied solid and liquid carbon and aluminum over three orders of density. With our AvIon method, we find excellent agreement with DFT predictions for the continuum edge shift below the Fermi energy, the carbon K-edge increase and the valence band gaps. Traditional AA models on the other hand predict valence bands to enter the continuum at far too low densities. Still, our AvIon results agree with a continuum lowering predicted by SP model up to a few times solid density. But at higher densities, drastic errors in the SP models become apparent as the continuum band shifts away from the potential threshold that is incorrectly equated with the beginning of the continuum band in the SP model.
Such predictions may be tested in experiment at high energy laser facilities like the National Ignition Facility or Laser Mégajoule (LMJ) [65] that are capable of producing compressed matter of Gbar pressures [66] and multiple times solid density [67] and may be combined with x-ray diffraction [68] and x-ray Thomson scattering [69].
Since our AvIon approach is as efficient as traditional AA models and can deliver DFT compatible predictions in only a few CPU seconds, it allows to fit and interpret experimental data and infer the temperature and density without resorting to SP models. The capability of AvIon will make joint experimental-theoretical explorations of WDM more efficient.