Transition to metallization in warm dense helium-hydrogen mixtures using stochastic density functional theory within the Kubo-Greenwood formalism

The Kubo-Greenwood (KG) formula is often used in conjunction with Kohn-Sham (KS) density functional theory (DFT) to compute the optical conductivity, particularly for warm dense matter. For applying the KG formula, all KS eigenstates and eigenvalues up to an energy cutoff are required and thus the approach becomes expensive, especially for high temperatures and large systems, scaling cubically with both system size and temperature. Here, we develop an approach to calculate the KS conductivity within the stochastic DFT framework, which requires knowledge only of the KS Hamiltonian but not its eigenstates and values. We show that the computational effort associated with the method scales linearly with system size and reduces in proportion to the temperature, unlike the cubic increase with traditional deterministic approaches. In addition, we find that the method allows an accurate description of the entire spectrum, including the high-frequency range, unlike the deterministic method which is compelled to introduce a high-frequency cutoff due to memory and computational time constraints. We apply the method to helium-hydrogen mixtures in the warm dense matter regime at temperatures of $\sim 60\mathrm{kK}$ and find that the system displays two conductivity phases, where a transition from nonmetal to metal occurs when hydrogen atoms constitute $\sim 0.3$ of the total atoms in the system.

Figure 1

The conductivity $\sigma \left(\omega \right)$ and DOS $\rho \left(\varepsilon \right)$ (in the insets) using stochastic versus deterministic (Quantum Espresso) calculations. We show two examples, each calculated in a periodic simulation cell of length $L$ at the $\mathrm{\Gamma }$ point: an insulator He atom at 3200 K with $L=5.3\text{Å}$ (upper panel) and a metallic system, ${\text{H}}_{256}$ at 4000 K, with $L=8\text{Å}$ (lower panel). The deterministic QE results used 200 KS eigenstates for the first system and 1700 for the second. For the stochastic DFT calculation we used 960 stochastic orbitals for the insulator and 480 orbitals for the conductor. The conductivity of both systems was calculated using 120 stochastic orbitals. A kinetic energy cutoff of $762\text{eV}$ for He and 525 eV for ${\text{H}}_{256}$ was used. Each peak was Gaussian broadened, deploying a width parameter equal to $\eta =1.2$ eV [in sKG this parameter affects Eq. (2)].

Figure 2

The conductivity spectrum for a ${\text{He}}_{128}$ system at $T=27\text{kK}$ and density of $0.83\frac{\text{g}}{{\text{cm}}^3}$. The sKG-sDFT LDA conductivity with error bars (with $I_{\sigma }=I_{\mathrm{H}}=160$) is compared to the deterministic results by vasp based on PBE, as described in Ref. [12]. The deterministic calculations were done with an energy cutoff of 800 eV using 570 KS states. Inset: The spectrum decay. The red line is proportional to ${\omega }^{-5/2}$.

Figure 3

The conductivity spectrum of ${\text{He}}_{128}$ at 27 kK and density of $0.83{\mathrm{g}/\mathrm{cm}}^3$. Top panel: The conductivity based on one sDFT Hamiltonian (using $I_{\mathrm{H}}=150$) calculated with an increasing number $I_{\sigma }$ of sKG-o's. Inset: The standard deviation (stdev) of the conductivity, averaged over all frequencies as a function of $I_{\sigma }$. The dashed line is proportional to $I_{\sigma }^{-1/2}$. Bottom panel: The conductivity based on three sDFT Hamiltonians, each obtained using $I_{\mathrm{H}}$ sDFT-o's. In order not to clutter the plot, error bars are given only for the $I_{\mathrm{H}}=150$ sDFT-o's calculation. The sKG calculations were all done using $I_{\sigma }=150$ sKG-o's.

Figure 4

The CPU time (blue circled markers) and $N_C\left(n\mathrm{\Delta }t\right)/n$ (purple square markers) as a function of the number $n$ of time step propagation operators used in the calculation. The calculation was done for ${\text{He}}_{128}$ using $N_g={60}^3$ grid points, at 57 kK on 8 processors, where $N_{\mathrm{ts}}=128$ and $\mathrm{\Delta }t=0.25$ a.u., for one dipole direction.

Figure 5

Wall time for stochastic $\mathrm{DFT}+\mathrm{KG}$ calculations with 600 time steps with $\mathrm{\Delta }t=0.25\hslash E_h^{-1}$ as well as a deterministic calculation done using Quantum Espresso (QE), as a function of the number of He atoms, at $0.75\frac{\text{g}}{{\text{cm}}^3}$ and 57 kK. The orange curve is proportional to $N^3$, and the blue curve is linear with $N$, the number of atoms.

Figure 6

Comparison of the calculated conductivity $\sigma \left(\omega \right)$ (a), the DOS $\rho \left(\varepsilon \right)$ (b), and the radial pair correlation $g\left(r\right)$ (c) for ${\text{He}}_{128}$ at 57 kK and density of $0.75\frac{\text{g}}{{\text{cm}}^3}$ based on configurations generated by AIMD [12] versus eFF dynamics. The dashed line in panel (b) shows the Fermi-Dirac level occupation, transiting from a value of 1 at low energies $\varepsilon$ to 0 at high energies.

Figure 7

The conductivity (upper panel) and the DOS (lower panel) of different systems containing 1024 atoms with different hydrogen percentages at $T=57\text{kK}$ and an average atomic volume of $60a_0^3$ per atom. The conductivity was normalized according to the number of electrons in the system. The DOS is shifted so that the chemical potential is zero.

Figure 8

The maximal conductivity frequency ${\omega }_{\max }$ (top panel), maximal conductivity ${\sigma }_{\max }$ (middle panel), and the DC conductivity of the actual (round blue markers) and the linear-mixing model ${\sigma }_{\text{LM}}$ (square orange markers), as a function of the hydrogen ratio $X_{\mathrm{H}}$ in He-H mixtures.