Equation of state and shock compression of warm dense sodium - a first-principles study

TL;DR

This study combines first-principles methods to accurately determine sodium's equation of state across wide conditions, revealing ionization effects and benchmarking existing models for shock compression behavior.

Contribution

It presents a coherent EOS for sodium using PIMC and DFT-MD, identifying ionization-induced compression maxima and providing benchmark data for popular EOS models.

Findings

Identified two compression maxima due to K-shell and L-shell ionization.

Provided benchmark EOS data for sodium over wide density-temperature ranges.

Showed radiation effects dominate at very high temperatures, affecting shock compression.

Abstract

As one of the simple alkali metals, sodium has been of fundamental interest for shock physics experiments, but knowledge of its equation of state (EOS) in hot, dense regimes is not well known. By combining path integral Monte Carlo (PIMC) results for partially-ionized states [B. Militzer and K. P. Driver, Phys. Rev. Lett. 115, 176403 (2015)] at high temperatures and density functional theory molecular dynamics (DFT-MD) results at lower temperatures, we have constructed a coherent equation of state for sodium over a wide density-temperature range of $1.93-11.60$ g/cm$^{3}$ and $10^3-1.29\times10^8$ K. We find that a localized, Hartree-Fock nodal structure in PIMC yields pressures and internal energies that are consistent with DFT-MD at intermediate temperatures of $2\times10^6$ K. Since PIMC and DFT-MD provide a first-principles treatment of electron shell and excitation effects, we are able to identify two compression maxima in the shock Hugoniot curve corresponding to $K$-shell and $L$-shell ionization. Our Hugoniot curves provide a benchmark for widely-used EOS models, SESAME, LEOS, and Purgatorio. Due to the low ambient density, sodium has an unusually high first compression maximum along the shock Hugoniot curve. At beyond 10$^7$ K, we show that the radiation effect leads to very high compression along the Hugoniot curve, surpassing relativistic corrections, and observe an increasing deviation of the shock and particle velocities from a linear relation. We also compute the temperature-density dependence of thermal and pressure ionization processes.

Tables

Table 1. Table 1: Comparison of energies and pressures calculated using different methods. For the PIMC calculations, we also specify the node type, the basis of localized orbitals and their numbers n s subscript 𝑛 𝑠 n_{s} . The internal energies are in units of Ha/atom, pressures are in GPa.

Method	Basis	$n_s$	$E$	${\sigma }_E$	$E$-$E_{\text{L}DA}$	$P$	${\sigma }_P$	$P$-$P_{\text{L}DA}$
$T$=2020958 K
LDA			-60.61			11354.6
PBE			-61.32		-0.71	11349.6		-4.9
PIMC-FP		0	-58.67	0.46	1.93	11356.4	68.1	1.8
PIMC-HF	6-31G	6	-59.06	0.37	1.55	11357.7	56.0	3.1
PIMC-HF	6-31G	13	-59.45	0.19	1.16	11313.6	29.7	-41.0
PIMC-HF	6-31+G	13	-59.15	0.27	1.46	11362.1	42.6	7.5
PIMC-HF	6-31+G	17	-59.13	0.31	1.48	11360.8	49.3	6.2
PIMC-LDA	6-31G	13	-59.77	0.43	0.83	11283.4	66.3	-71.2
PIMC-PBE	6-31G	13	-59.28	0.34	1.32	11340.0	53.0	-14.6
PIMC-PBEX	6-31G	13	-58.92	0.39	1.69	11405.4	59.9	50.8
$T$=1010479 K
LDA			-112.87	0.00	0.00	4528.0	0.0	0.0
PBE			-113.54	0.00	-0.67	4531.8	0.0	3.8
PIMC-FP		0	-111.14	0.80	1.73	4474.0	125.2	-54.0
PIMC-HF	6-31G	6	-113.12	0.33	-0.24	4544.7	51.5	16.7
PIMC-HF	6-31G	13	-113.44	0.29	-0.57	4508.9	45.8	-19.1
PIMC-HF	6-31+G	13	-113.47	0.23	-0.60	4508.2	37.0	-19.8
PIMC-HF	6-31+G	17	-112.73	0.31	0.15	4627.1	49.9	99.1
PIMC-LDA	6-31G	13	-113.53	0.36	-0.65	4473.6	55.7	-54.4
PIMC-PBE	6-31G	13	-113.67	0.32	-0.80	4466.6	52.6	-61.4
PIMC-PBEX	6-31G	13	-113.38	0.32	-0.51	4507.5	50.8	-20.5

Equations

$${\rho }_F(\text{R},{\text{R}}^{\prime };\beta )=\frac{1}{N!}\sum _{\mathcal{P}}(-1{)}^{\mathcal{P}}{\int }_{\text{R}\to \mathcal{P}{\text{R}}^{\prime },{\rho }_T>0}d{\text{R}}_t{e}^{-S[{\text{R}}_t]},$$

$${\rho }_T(\text{R},{\text{R}}^{\prime };\beta )={\left|\left|{\rho }^{[1]}({\text{r}}_i,{\text{r}}_j^{\prime };\beta )\right|\right|}_{ij},$$

$${\rho }_0^{[1]}(\text{r},{\text{r}}^{\prime };\beta )=\sum _{\text{k}}{e}^{-\beta E_{\text{k}}}{\Psi }_{\text{k}}(\text{r}){\Psi }_{\text{k}}^{\ast }({\text{r}}^{\prime }),$$

$$\begin{array}{rl}{\rho }^{[1]}(\text{r},{\text{r}}^{\prime },\beta )& ={\rho }_0^{[1]}(\text{r},{\text{r}}^{\prime };\beta )+\\ & \sum _{I=1}^N\sum _{s=0}^{n_s}{e}^{-\beta E_s}{\Psi }_s(\text{r}-{\text{R}}_I){\Psi }_s^{\ast }({\text{r}}^{\prime }-{\text{R}}_I).\end{array}$$

$$\mathcal{H}(T,\rho )=(E-E_0)+\frac{1}{2}(P+P_0)(V-V_0)=0,$$

$$u_s^2=\xi /\eta \text{and}{u}_p^2=\sqrt{\xi \eta },$$

$$P_{\text{radiation}}=\frac{4\sigma }{3c}{T}^4\text{and}{E}_{\text{radiation}}=3P_{\text{radiation}}V,$$

$$N(r)=⟨\frac{1}{N_I}\sum _{e,I}\theta (r-\left|{\vec{r}}_e-{\vec{r}}_I\right|)⟩,$$

Full text

I Introduction

The behavior of sodium and other alkali metals at extreme conditions has generated intense scientific interest over many decades as experimental and theoretical technology has evolved to facilitate studies of increasingly atypical states. At ambient conditions, sodium has long been known as a prototypical, simple, free-electron metal Wigner and Seitz (1933) with a high-symmetry, bulk-centered cubic (bcc) structure. Recent interest in sodium has been driven by the the ability to carefully probe exotic high-pressure states Neaton and Ashcroft (2001), made possible by improvements in static-compression Eremets (1996); Dubrovinsky et al. (2012) and single-crystal experimental techniques Gregoryanz et al. (2008). A number of experimental and theoretical works have revealed that sodium exhibits anomalous structural Aleksandrov et al. (1982); Fritz and Olinger (1984); Hanfland, Loa, and Syassen (2002); Raty, Schwegler, and Bonev (2007); McMahon et al. (2007); Gregoryanz et al. (2008); Rousseau et al. (2011), electronic and optical Raty, Schwegler, and Bonev (2007); Ma et al. (2009); Lazicki et al. (2009); Gatti, Tokatly, and Rubio (2010); Mao et al. (2011); Marqués et al. (2011); Loa et al. (2011); Pozzo, Desjarlais, and Alfè (2011); Ibañez Azpiroz et al. (2014); Miao and Hoffman (2015); Naumov and Hemley (2015), and melting Gregoryanz et al. (2005); Koči et al. (2008); Yamane, Shimojo, and Hoshino (2008); Maksimov, Lepeshkin, and Magnitskaya (2011); Marqués et al. (2011); Eshet et al. (2012); González and González (2015) behavior at high pressure. Among the structural studies, high-pressure X-ray diffraction experiments have explored sodium up to pressures and temperatures of 150 GPa and 1000 K, observing a complex series of phase transitions Hanfland, Loa, and Syassen (2002); McMahon et al. (2007); Gregoryanz et al. (2008) from the bcc phase to at least seven lower-symmetry phases. Sodium also displays an anomolous melting curve maximum Gregoryanz et al. (2005) in this pressure range as the liquid undergoes structural electronic changes, becoming more dense than the solid over a large pressure range of $\sim$60 GPa. At even higher pressures, beyond 200 GPa, sodium undergoes a metal-insulator transition Ma et al. (2009) by forming an electride Miao and Hoffman (2015).

Ab initio random structure searches predict that sodium will ultimately become a reentrant metal at a pressure of 15.5 TPa at 0 K Li et al. (2015).

Analogous to the impact that improved static-compression technology has had for alkali metals in the condensed matter regime, we can expect future dynamic compression research to routinely probe the gigabar (Gbar) Swift et al. (2012); Kritcher et al. (2016) pressures and provide new data on the behavior of warm dense matter (WDM). Shock measurements of hot, dense alkali metals are expected to be an important component of solving problems in astrophysics, planetary physics, stockpile stewardship, and inertial fusion. Thus far, only relatively low pressure ($P<228$ GPa Golyshev et al. (2011), $T<8168$ K) shock experiments have been performed on sodium. Early dynamic compression experiments were motivated by compressibility of metals and the simplicity of ambient sodium as a material for the study of fundamental shock physics behavior Rice (1965); Bakanova, Dudoladov, and Trunin (1965); Mar . More recent shock experiments have been motivated by the prospect of observing the metal-insulator transition Golyshev et al. (2011).

Complementary to the shock experiments, there have been a number of theoretical works aiming at calculating the equation of state (EOS) and shock Hugoniot curve using either analytic approaches Ehrenfeld, Krimsky, and Selvitella (1966); Hickey (1967) or semi-classical free-energy models Pastine (1967, 1968); Khishchenko (2008), fluid-variational theory Ross (1980), local pseudopotential theory Young and Ross (1984), or EAM potentials Belashchenko (2009, 2013a, 2013b). While each of these methods agrees reasonably well with the available low-pressure shock Hugoniot data ($P<100$ GPa), none possess the rigor that is necessary to treat the strong coupling, quantum degeneracy, and ionization physics associated with hot, dense plasma regimes, that Gbar-shock experiments will access.

Development of rigorous, first-principles frameworks, such as path integral Monte Carlo (PIMC) Pierleoni et al. (1994); Driver and Militzer (2012); Militzer and Driver (2015) or density functional theory molecular dynamics (DFT-MD) Hu et al. (2016); Lambert and Recoules (2012); Gao et al. (2016); Zhang et al. (2016a), for computing the EOS and behavior of materials in extreme pressure and temperature conditions is needed as Gbar-range shock experiments are emerging and data need to be interpreted. Recently, we have been developing PIMC for simulating heavier elements Militzer and Ceperley (2001); Militzer (2009); Driver and Militzer (2012); Benedict et al. (2014); Driver and Militzer (2015); Militzer and Driver (2015); Driver et al. (2015); Driver and Militzer (2016); Hu et al. (2016); Zhang et al. (2016b), which is efficient at high temperatures, and can be used to complement DFT-MD calculations that are efficient at comparatively low temperatures. Combined data from PIMC and DFT-MD provide a coherent EOS over a wide density-temperature range that spans the condensed matter, WDM, and plasma regimes. In the current work, we use PIMC and DFT-MD to compute the EOS, shock compression behavior, and plasma structure evolution of sodium across a much larger density-temperature range than has ever been studied in experiments or in theory. We provide our EOS table and shock Hugoniot curve as a theoretical benchmark in comparison with widely-used EOS database models and experiments.

The paper is organized as follows: Section II discusses details of our PIMC and DFT-MD simulation methods. In Sec. III, we discuss several results: (1) comparison of the performance of different choices for PIMC nodal surfaces, (2) the internal energies and pressures, (3) shock Hugoniot curves, and (4) the density-temperature evolution of ionization processes. Finally, in Sec. IV, we conclude.

II Simulation methods

A reliable theoretical scheme for simulating materials across a wide range of density and temperature conditions naturally involves treating physics appropriately in different regimes. In the limit of high temperature (above $10^{8}$ K), materials tend to behave as weakly-interacting plasmas because of strong thermal ionization. Under such conditions, the ideal Fermi gas model and the Debye-Hückel theory Debye and Huckel (1923) are good approximations. At low to intermediate temperatures ($T<10^{6}$ K), standard, orbital-based, Kohn-Sham DFT-MD is a suitable option to derive the EOS because it fully accounts the bonding effects and bound states within certain approximations of the exchange-correlation functional. However, the need to explicitly compute all partially occupied electronic orbitals causes DFT-MD to become computationally intractable beyond temperatures of roughly $10^{6}$ K. Recent work on orbital-free DFT Sjostrom and Daligault (2014); Hu et al. (2016); Karasiev, Calderín, and Trickey (2016); Gao et al. (2016) and extended DFT-MD with free electron approximation for high energies Zhang et al. (2016a) have made progress toward overcoming this difficulty, but their general applicability and accuracy needs to be further examined.

PIMC simulates all quantum, many-body exchange and correlations effects and provides the most natural formulation to compute accurate EOSs at high temperatures ($T>10^{6}$ K). In PIMC, the thermal density matrix is expressed as Feynman path integral Feynman (1953), which treats electrons and nuclei as quantum paths that are cyclic in the imaginary time $0\leq t\leq\beta=1/k_{\text{B}}T$, where $k_{\text{B}}$ is the Boltzmann constant. Therefore, PIMC becomes increasingly efficient at higher temperatures as paths become shorter and more classical in nature. However, application of PIMC to study real materials other than hydrogen Pierleoni et al. (1994); Magro et al. (1996); Militzer and Ceperley (2001); Militzer et al. (2001); Hu et al. (2010); Militzer and Ceperley (2000); Hu et al. (2011); Militzer, Magro, and Ceperley (1999); Militzer and Graham (2006), helium Militzer (2009, 2006), hydrogen-helium mixtures Militzer (2005), and one-component plasmas Militzer and Pollock (2005); Pollock and Militzer (2004), is difficult because of the complex fermion sign problem, nonlocal pseudopotentials, and complex nodal structures Ceperley (1991).

The sign problem in fermionic PIMC simulations is usually addressed with the fixed-node approximation Ceperley (1991) that restricts paths to positive regions of a trial density matrix, $\rho_{T}(\textbf{R},\textbf{R}_{t};t)>0$. The restricted path integral reads,

[TABLE]

where $\mathcal{P}$ denotes permutations of identical particles, and $S$ is the action that weights the paths. The most common approximation to the trial density matrix is a Slater determinant of single-particle density matrices,

[TABLE]

in combination with the free-particle (FP) density matrix,

[TABLE]

which represents a sum over all plane waves $\Psi_{\textbf{{k}}}(\textbf{{r}})$ with energy $E_{\textbf{{k}}}$.

PIMC with FP nodes gives exact results in the limit of high temperature Militzer (2009). Developments of this method have allowed for remarkable progress. Using FP nodes, we have successfully obtained the EOS of several heavy elements (C Driver and Militzer (2012), N Driver and Militzer (2016), O Driver et al. (2015), Ne Driver and Militzer (2015)) and compounds (H2O Driver and Militzer (2012)). The combined results of PIMC and DFT-MD have bridged the WDM gap between DFT-MD and the high temperature limit and provided consistent sets of coherent EOS from first principles.

Several attempts have been made to go beyond FP nodes in PIMC simulations Kha ; Militzer and Ceperley (2000). In recent work Militzer and Driver (2015), we show that the applicability range of PIMC simulations can extend to lower temperatures when a number of $n_{s}$ atomic orbitals at each ion $I$ are added to the FP nodes,

[TABLE]

Our results for a single silicon atom in periodic boundary condition showed that nodes derived from Hartree-Fock (HF) orbitals yield highly accurate predictions for the pressure and the internal energy at much lower temperature than is accessible with FP nodes. The combined results using this PIMC method and DFT-MD provided a coherent EOS for dense silicon plasmas over a wide density-temperature grid (2.3-18.6 g/cm3, $5\times 10^{5}$-$1.3\times 10^{8}$ K). In this work, we will also investigate the effects of various nodal surfaces in PIMC calculations and show that the localized, Hartree-Fock orbitals yield accurate pressures and internal energies for sodium.

Our PIMC simulations are performed within the fixed node approximation Ceperley (1991) and based on the CUPID code Militzer (2000). Similar to the PIMC simulations of Si Militzer and Driver (2015) and recently of Na at one density Zhang et al. (2016b), we treat the nuclei classically because of the high temperatures considered here. Electronic Coulomb interactions are introduced via pair density matrices Militzer (2016).

DFT-MD simulations use the Vienna Ab initio Simulation Package (VASP) Kresse and Furthmüller (1996) and implement exchange-correlation functionals within the local density approximation (LDA) Perdew and Zunger (1981); Ceperley and Alder (1980). The effect of using other forms of exchange-correlation functionals is small, in comparison to the effects of thermal excitations, as will be disscussed in Sec. III.1. We use a projected augmented wave (PAW) pseudopotential Blöchl, Jepsen, and Andersen (1994) with a frozen, 1s2 electrons core and a small core radius of 1.45 Bohr, which is the hardest one available for Na in VASP com .

Since the width of the electronic bands is small compared to the thermal excitation energies, we use the $\Gamma$ point for sampling the Brillouin zone. We choose plane wave basis with a large cutoff of 4000 eV to account for the electrons that are excited to the very high-energy states. We use a time step of $0.2$ fs and perform DFT-MD simulations typically for over 2000 steps. A large number of temperature, energy, and pressure fluctuations are observed, which indicate the simulations have reached equilibrium. We discard 20% of the trajectories at the beginning and average the internal energy and pressure to derive the EOS. In order to put the PAW-LDA pseudopotential energies on the same scale as all-electron calculations, we shifted all of our VASP DFT-MD energies by -161.3386 Ha/atom. This shift was determined by performing isolated, all-electron atomic calculations with the OPIUM code opi and corresponding isolated-atom calculations in VASP.

Most of our simulations start from a bcc cell with 8 atoms. In comparison with using a 54-atom cell, DFT-MD simulations with the 8-atom cell is sufficient in providing well converged internal energies and pressures down to temperatures below 5.0$\times 10^{4}$ K. For example, the difference is 0.1 eV/atom in internal energy and 0.3 GPa in pressure between simulations using an 8- and a 54-atom cells, at 3.87 g/cm3 and 5.05$\times 10^{4}$ K. At lower temperatures, the system could tend to freeze and is subject to larger finite-size errors. In order to maximally eliminate this effect, we use a larger, 54-atom cell at these temperatures up to 2.0$\times 10^{5}$ K, as long as the computational cost is affordable. We simulate data along nine isochores corresponding to 2-, 3-, 4-, 5-, 6-, 7-, 8-, 10-, and 12-fold compression of ambient density $\rho_{\text{ambient}}$=0.96663 g/cm3. For each density, we study the temperature range from 103 to 1.29$\times 10^{8}$ K, relevant to the regime of WDM and stellar interiors (Figs. 1-2).

III Results

III.1 Comparison of PIMC methods

In Ref. Militzer and Driver, 2015, we have shown that accuracy of PIMC simulations can be improved by introducing atomic orbitals, derived with the HF method, to the fermion nodes. Here we investigate whether any further improvement can be made by representing the orbitals with more accurate basis sets, including a large number of localized orbitals or by deriving the orbitals with LDA or generalized gradient approximation in the Perdew-Burke-Ernzerhof (PBE) type Perdew, Burke, and Ernzerhof (1996), instead of HF. Before performing many-atom simulations, we tested various methods of introducing atomic orbitals to study one stationary sodium atom in a periodic 5-Bohr cubic box, and compared the results in Fig. 3 and Table 1.

Figure 3 shows the difference in internal energy and pressure between PIMC and DFT-MD calculations. At temperatures of $T>2\times 10^{6}$ K, all PIMC energies and pressures are systematically higher than our DFT-MD results within LDA or PBE. This can be attributed to the excitations of 1s electrons, which are treated explicitly as the other outer-shell electrons in PIMC, but are frozen in the pseudopotential of DFT-MD calculations.

At temperatures of $T\leq 2\times 10^{6}$ K, the energies from PIMC computations with FP nodes are systematically too high, while the pressures remain quite reasonable. In contrast, using PIMC with HF nodes, the results are in much better agreement with PBE predictions. This suggests that FP nodes do not lead to the correct $K$ and $L$ shell structures of the atom, which is primarily a local property, so the error in the energy does not depend significantly on density. At $T<10^{6}$ K, pressures from our PIMC+13 HF simulations are high, so are the uncertainties. This is due to the low sampling efficiency of PIMC at these temperatures.

Table 1 lists additional PIMC results for two temperatures of 1 and 2 million Kelvin. In all cases, the orbitals were derived with spin-restricted GAMESS gam calculations of the Na*+* ion so that we can use the same orbitals for both spin channels in PIMC calculations. This is a reasonable approximation because the spin state is of minor importance at high temperatures. This is confirmed by DFT-MD calculations: within LDA and PBE, the spin-polarized (5+6) and spin-unpolarized calculations yield similar results for the temperature range under consideration.

In our PIMC calculations with localized nodal surfaces, we need a minimum of 6 atomic orbitals to provide at least one bound state for every electron. Table 1 shows that we found no statistically significant difference between using 6 and 13 HF orbitals in PIMC calculations, for both temperatures. At 1010479 K, the PIMC pressure and energy become too high when we increase the number of orbitals to 17. We attribute this deviation to our small simulation cell of 5.0 Bohr. The highest atomic orbitals are too delocalized to fit into this cell.

We find no statistically significant differences in the PIMC results between using a 6-31G and a slightly more accurate 6-31+G basis set. Both yield similar HF energies that are 0.176 Ha higher than basis set-converged HF calculations with a TZV basis. This difference is within the error bar of typical PIMC calculations.

We also test whether there would be an advantage in deriving the atomic orbitals with LDA or PBE methods rather than with HF theory, and find no significant difference. Furthermore, our studies with orbitals that are derived with just PBE exchange (PBEX) lead to similar PIMC results.

In the following many-atom PIMC calculations, we choose 13 HF orbitals that are generated with 6-31G basis in Eq. 4.

III.2 Equation of state

Figure 4 shows the calculated internal energies and pressures along nine isochores relative to the ideal Fermi electron gas model 111Here, we include the kinetic energy of the nuclei when calculating the total internal energy of the ideal Fermi electron gas model, but have not considered the charge and interaction of the nuclei.. Our results show excellent agreement between PIMC and DFT-MD at 2$\times$106 K for all densities. The difference between PIMC and DFT-MD is typically less than $3$ Ha/atom in internal energy and within 3% in pressure. We have thus succeeded in constructing a coherent first-principles EOS table for warm dense sodium, over a wide range of densities (2-12 $\rho_{\text{ambient}}$) and temperatures ($10^{3}$-$1.29\times 10^{8}$ K).

At high temperatures of $T>2\times 10^{7}$ K, our PIMC results agree with those from the Debye-Hückel model as well as the Fermi electron gas theory, because the temperatures are so high that the atoms are fully ionized.

PIMC successfully bridges DFT-MD at 2$\times$106 K and analytical models in the high-temperature limit. This cross-validates the EOS data from present PIMC simulations with fixed-node approximation, and the use of a frozen core and a zero-temperature exchange-correlation functional in DFT-MD up to $2\times 10^{6}$ K.

As temperature decreases, the ideal Fermi electron gas model significantly overestimates the energy and the pressure because it neglects all interactions. This is partially improved in the Debye-Hückel model, which treats weak interactions correctly within the screening approximation. However, because the Debye-Hückel model still does not treat bound states, it leads to unphysically low pressures and energies at low temperatures, as the screening approximation breaks down and electrons occupy bound states.

In addition, Fig. 5 shows the nuclear pair-correlation functions $g_{\text{Na-Na}}(r)$ computed in DFT-MD and PIMC simulations using an 8-atom simulation cell at $10^{6}$ K and 8-fold compression. The good agreement in the $g_{\text{Na-Na}}(r)$ between the two methods shows that PIMC and DFT-MD predict consistent ionic plasma structures. This is further confirmation that the fixed-node approximation in PIMC and the exchange-correlation and pseudopotential approximation in DFT-MD do not inhibit the accuracy of these methods in the WDM regime.

III.3 Shock compression

In dynamic compression experiments, the conservation of mass, momentum, and energy constrains a steady shock to follow the Rankine-Hugoniot equation Ahrens (2003)

[TABLE]

where ($E_{0},P_{0},V_{0}$) and ($E,P,V$) represents the internal energy, pressure, and volume of the initial and the shocked state, respectively. The EOS ($E$ and $P$ on a grid of $T$ and $V$) data in the previous section allows us to solve Eq. 5 with spline fitting.

Given that samples in shock experiments may be pre-compressed to reach higher-density, lower-temperature states off the principal Hugoniot, we consider four different initial conditions corresponding to 0.75, 1.0, 1.5, and 2.0 times $\rho_{\textrm{ambient}}$. DFT simulations are performed at each of these densities to determined the corresponding initial pressures and internal energies.

We thus calculate the Hugoniot curves and represent them with pressure-density plots in Fig. 6. Two compression peak maxima are predicted along each of the Hugoniot curves, one above $2\times 10^{6}$ K and the other below $10^{6}$ K. We tested different bivariate interpolation methods for the EOS grid in $\rho-T$ space and observed consistent compression curves. Small discrepancies ($<$0.1) in the compression minimum in between the two peaks is observed. We attribute this to the non-smoothness caused by the small differences between EOS values obtained by PIMC and DFT-MD at the corresponding pressure and temperature conditions.

For comparison, we consider two shock compression profiles predicted by SESAME Lyon and Johnson (1992), which is a tabular database for the thermodynamic properties of materials that is constructed by connecting available shock wave data with Thomas-Fermi-Dirac and Mie-Grüneisen theory at high densities, and some simple analytic forms at low densities. The high-temperature limit is sufficiently described by the Thomas-Fermi-Dirac theory, as is shown by the remarkable consistency between the SESAME database and our present first-principles calculations at $T>2\times 10^{7}$ K. However, at lower temperatures when the system includes bound states, Fermi-Dirac breaks down and the SESAME database becomes insufficient. SESAME therefore fails to capture shell effects and only exhibits a single density peak along each Hugoniot curve.

Additionally, Fig. 7 shows temperature and pressure along the Hugoniot curve as functions of the shock compression ratio $\rho/\rho_{0}$. When the initial density is the lowest ($\rho_{0}=0.75~{}\rho_{\text{ambient}}$), a maximum compression ratio of $\sim$5.6 is reached at 6$\times$105 K. The value of this peak decreases with increasing $\rho_{0}$, and reduces to 4.7 when $\rho_{0}=2.0\rho_{\text{ambient}}$. The higher-temperature compression peak reaches slightly larger compression ratio (4.9) than the lower peak as the initial density increases. In the high-temperature limit, the system is almost fully ionized and approaches the the non-relativistic ideal limit ($\rho/\rho_{0}=4.0$), regardless of the initial density.

Figure 8 compares the principal Hugoniot curve of several elements including He Militzer (2009), C Driver and Militzer (2012), N Driver and Militzer (2016), O Driver et al. (2015), Ne Driver and Militzer (2015), and Si Militzer and Driver (2015). Interestingly, sodium together with helium show higher maximum compression ratio than other elements. We attribute this to the unusual low ambient density of sodium. In addition, each element imprints a particular structural signature in its principal Hugoniot curve, exhibiting compression maxima with different values and at varied pressures and temperatures, due to the interplay of excitation of internal degrees of freedom, electronic interactions, and the interaction effects tied to the initial conditions Militzer (2006).

Figure 1 compares the principal Hugoniot curve of sodium from our first-principles calculations and those predicted by widely-used EOS models, SESAME Lyon and Johnson (1992), LEOS More et al. (1988); Young and Corey (1995), and average-atom Purgatorio (Lynx) Wilson et al. (2006); Whi . The compression maxima along the Hugoniot curves are closely related to ionization and interaction effects. SESAME and LEOS do not explicitly include information about electronic shell structure, and therefore do not show two distinct compression maxima. On the other hand, the DFT-based, average-atom Purgatorio (Lynx) model does compute the shell structure for an average of multiple ionization states. Therefore, Purgatorio (Lynx) agrees well with our first-principles results at above 200 GPa or 105 K. Below that, it is less reliable, because average-atom approaches cannot treat bonding and many body effects in a dense fluid properly. This is demonstrated by the deviation of the Hugoniot curve by Purgatorio (Lynx) from our calculations, which agrees with those predicted by SESAME and LEOS database, which were constructed by extrapolating experimental values.

We also compare our principal Hugoniot results with available experimental Rice (1965); Bakanova, Dudoladov, and Trunin (1965); Mar and theoretical Young and Ross (1984); Belashchenko (2013a) data available at low temperatures and pressures. In comparison with the experimental data, the local pseudopotential theory Young and Ross (1984) underestimates the pressure by up to 30 GPa. The EAM model Belashchenko (2013a) agrees well with experiments up to 110 GPa, above which the model becomes invalid. Notably, DFT-MD, LEOS, and SESAME lie 10-20 GPa below the higest-pressure experimental data, but all methods tend to agree with the experimental data at the lowest pressures.

The shock velocity $u_{s}$ and the particle velocity $u_{p}$ that are of interests in shock experiments can be calculated using Militzer et al. (2001)

[TABLE]

where $\xi=(P-P_{0})/\rho_{0}$, $\eta=1-\rho_{0}/\rho$. Figure 9 compares the relation between $u_{s}$ and $u_{p}$ of sodium from our first-principles calculations and those predicted by the SESAME database, low-pressure ($P<10^{2}$ GPa) shock experiments on sodium Bakanova, Dudoladov, and Trunin (1965); Mar , and an “universal Hugoniot” Ozaki et al. (2016) obtained from high-pressure (up to 0.2 Gbar) experimental Hugoniot data on Al, Fe, Cu, and Mo. We find the $u_{s}-u_{p}$ relation is well consistent with the linear relations predicted by experiments at $u_{s}<15$ km/s Bakanova, Dudoladov, and Trunin (1965); Mar , but increasingly deviates from them when the shock velocities exceeds $\sim$500 km/s. This high velocity corresponds to a shock pressure of $\sim$1.9 Gbar along the principal Hugoniot curve, and temperature of $\sim 10^{7}$ K, at which the atoms are nearly fully ionized (see the discussion in Sec. III.4). The deviation from the “universal Hugoniot” of fluid metals Ozaki et al. (2016) is more evident. This reflects a fundamental difference between the effects of ionization on shock compression in sodium and in heavier metals, such as Fe, Cu, and Mo. The $u_{s}-u_{p}$ profile by SESAME agrees remarkably well with our first-principals predictions, which is not unexpected because of the nearly linear relation between $u_{s}$ and $u_{p}$ and the consistency between SESAME and our Hugoniot data at low ($<10^{4}$ K) and high ($>10^{7}$ K) temperature and pressure conditions (Figs. 6-7). Starting from pre-compressed sodium at 2 times ambient density, the shock velocity is higher but the slope $du_{s}/du_{p}$ remains 1.31, according to linear interpolation.

When temperature exceeds 107 K, the radiation contribution becomes important. In order to evaluate the radiation effect on shock compression, we consider an ideal black body scenario and add the photon contribution to the EOS using

[TABLE]

where $\sigma$ is the Stefan-Boltzmann constant and $c$ is the speed of light in vacuum. We then re-construct the Hugoniot curves and show them Figs. 1 and 7. We find that the two compression peaks remain unchanged as we include radiative effects. However, the Hugoniot curves deviate significantly from the classical limit above $10^{7}$ K and tend towards a compression ratio asymptote of 7. These results imply that the radiation contribution plays a significant role in the shock compression of materials at extreme temperatures ($T>10^{7}$ K) and pressures ($P>1$ Gbar). In comparison, the pressure-contribution from relativistic effects, included in the Purgatorio (Lynx) model is much smaller (Fig. 1).

III.4 Ionization

The structure of the Hugoniot curves observed in the previous section can be understood from the density-temperature dependence of the ionization process. In order to examine this relation, Fig. 10 shows the number of electrons near the nucleus for a given temperature and density, $N_{\textrm{Na}-e}(r)$, given by the formula

[TABLE]

where the sum includes all electron-ion pairs and $\theta$ represents the Heaviside function.

The $N_{\textrm{Na}-e}(r)$ functions at two different densities and a series of temperatures are compared with ionization states of isolated sodium ions (calculated with GAMESS). The results show that $N_{\textrm{Na}-e}(r)$ decreases with $T$, indicating the gradual ionization of the atoms with temperature. The energy levels of the 1s electrons are much deeper than those of the outer-shell electrons, and therefore only become excited above $2\times 10^{6}$ K. This indicates that the upper peaks on the Hugoniot curve in Figs. 6 and 7 are related to the excitation of the $K$-shell electrons, while the lower peaks correspond to the excitation of $L$-shell electrons.

IV Conclusions

In this work, we construct a thermal density matrix with localized, HF orbitals to construct fermion nodal surfaces and perform PIMC simulations of the second-row element sodium. We obtain an accurate EOS from temperatures of 129 million to 1 million K, at which the results are consistent with DFT-MD. The excellent agreement between the PIMC and DFT-MD validates the use of the pseudopotential that freezes 1s electrons and the use of zero-temperature exchange-correlation functionals up to temperateures of $2\times 10^{6}$ K.

By investigating the shock compression curves using the obtained EOS data, we find two compression maxima along the Hugoniot curves. This is in contrast to the single-peak Hugoniot curve predicted by the SESAME and the LEOS database, which are based on models that neglect bonding and many-body effects. The higher-temperature compression maxima occurs at $6\times 10^{6}$ K due to the significant excitation of 1s electrons, while the lower compression maximum is due to the thermal excitation of outer-shell electrons. We predict a maximum compression ratio of 5.3 along the principal Hugoniot, and compression ratios of greater than 5 can be reached when the initial density is less than 1.5 $\rho_{\text{ambient}}$. The value decreases at higher initial densities due to stronger particle interaction, and vice versa. Including a radiation contribution, shock Hugonits are modified significantly above 107 K and 1 Gbar and a much higher compression ratio can be reached.

V Supplementary Material

See the supplementary material for the tables of EOS data.

Acknowledgements.

This research is supported by the U. S. Department of Energy, grant DE-SC0010517 and DE-SC0016248. Shuai Zhang is partially supported by the PLS-Postdoctoral Grant of the Lawrence Livermore National Laboratory. Computational support was provided by the Janus supercomputer, which is supported by the National Science Foundation (Grant No. CNS-0821794), the University of Colorado, and the National Center for Atmospheric Research, and NERSC. This research is part of the Blue Waters sustained-petascale computing project (NSF ACI 1640776), which is supported by the National Science Foundation (awards OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. S.Z. appreciates the helpful discussion with Mu Li.