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

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¹ with a high-symmetry, bulk-centered cubic (bcc) structure. Recent interest in sodium has been driven by the ability to carefully probe exotic high-pressure states,² made possible by improvements in static-compression^3,4 and single-crystal experimental techniques.⁵ A number of experimental and theoretical works have revealed that sodium exhibits anomalous structural,^5–11 electronic and optical,^9–21 and melting^16–27 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^5,8,10 from the bcc phase to at least seven lower-symmetry phases. Sodium also displays an anomalous melting curve maximum²² in this pressure range as the liquid undergoes structural electronic changes, becoming more dense than the solid over a large pressure range of ∼60 GPa. At even higher pressures, beyond 200 GPa, sodium undergoes a metal-insulator transition¹² by forming an electride.²⁰ Ab initio random structure searches predict that sodium will ultimately become a re-entrant metal at a pressure of 15.5 TPa at 0 K.²⁸ 

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)^29,30 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,³¹ $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.^32–34 More recent shock experiments have been motivated by the prospect of observing the metal-insulator transition.³¹ 

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^35,36 or semi-classical free-energy models,^37–39 fluid-variational theory,⁴⁰ local pseudopotential theory,⁴¹ or embedded atom model (EAM).^42–44 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)^45–47 or density functional theory molecular dynamics (DFT-MD),^48–51 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,^{46–48,52–58} 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.

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⁸ 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⁵⁹ are good approximations. At low to intermediate temperatures (⁠$T< 106$ 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⁶ K. Recent works on orbital-free DFT^48,50,60,61 and extended DFT-MD with free electron approximation for high energies⁵¹ have made progress toward overcoming this difficulty, but their general applicability and accuracy need 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> 106$ K). In PIMC, the thermal density matrix is expressed as the Feynman path integral,⁶² which treats electrons and nuclei as quantum paths that are cyclic in the imaginary time $0≤t≤β=1/kBT$⁠, where $kB$ is the Boltzmann constant. Therefore, PIMC becomes increasingly efficient at higher temperatures as paths become shorter and more classical in nature. However, the application of PIMC to study real materials other than hydrogen,^{45,52,63–69} helium,^53,70 hydrogen-helium mixtures,⁷¹ and one-component plasmas^72,73 is difficult because of the complex fermion sign problem, nonlocal pseudopotentials, and complex nodal structures.⁷⁴ 

The sign problem in fermionic PIMC simulations is usually addressed with the fixed-node approximation⁷⁴ that restricts paths to positive regions of a trial density matrix, $ρT(𝐑,𝐑t;t)>0$⁠. The restricted path integral reads

where $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,

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

which represents a sum over all plane waves $Ψ𝒌(𝒓)$ with energy Ek.

PIMC with FP nodes gives exact results in the limit of high temperature.⁵³ Developments of this method have allowed for remarkable progress. Using FP nodes, we have successfully obtained the EOS of several heavy elements (C,⁴⁶ N,⁵⁷ O,⁵⁶ and Ne⁵⁵) and compounds (H₂O⁴⁶). 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.^66,75 In a recent work,⁴⁷ 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,

Our results for a single silicon atom in the 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 temperatures than those which are 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/cm³, $5×105$-$1.3×108$ 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⁷⁴ and are based on the CUPID code.⁷⁶ Similar to the PIMC simulations of Si⁴⁷ and recently of Na at one density,⁵⁸ we treat the nuclei classically because of the high temperatures considered here. Electronic Coulomb interactions are introduced via pair density matrices.⁷⁷ 

DFT-MD simulations use the Vienna Ab initio Simulation Package (VASP)⁷⁸ and implement exchange-correlation functionals within the local density approximation (LDA).^79,80 The effect of using other forms of exchange-correlation functionals is small, in comparison to the effects of thermal excitations, as will be discussed in Sec. III A. We use a projected augmented wave (PAW) pseudopotential⁸¹ with a frozen, 1s² electrons’ core and a small core radius of 1.45 bohrs, which is the hardest one available for Na in the VASP.⁸² 

Since the width of the electronic bands is small compared to the thermal excitation energies, we use the Γ 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 that 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⁸³ and corresponding isolated-atom calculations in the 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 $× 104$ K. For example, the difference is 0.1 eV/atom in internal energy and 0.3 GPa in pressure between simulations using an 8-atom cell and a 54-atom cell, at 3.87 g/cm³ and 5.05 $× 104$ 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 $× 105$ 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 $ρambient$ = 0.966 63 g/cm³. For each density, we study the temperature range from 10³ to 1.29 $× 108$ K, relevant to the regime of WDM and stellar interiors (Figs. 1 and 2).

In Ref. 47, we have shown that the 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,⁸⁶ 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-bohrs cubic box and compared the results in Fig. 3 and Table I.

Figure 3 shows the difference in the internal energy and pressure between PIMC and DFT-MD calculations. At temperatures of $T>2×106$ 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≤2×106$ 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<106$ 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 I lists additional PIMC results for two temperatures of 1 and 2 × 10⁶ K. In all cases, the orbitals were derived with spin-restricted General Atomic and Molecular Electronic Structure System (GAMESS)⁸⁷ 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 I shows that we found no statistically significant difference between using 6 and 13 HF orbitals in PIMC calculations, for both temperatures. At 1 010 479 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 bohrs. 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).

Figure 4 shows the calculated internal energies and pressures along nine isochores relative to the ideal Fermi electron gas model.⁸⁸ Our results show excellent agreement between PIMC and DFT-MD at 2 × 10⁶ 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 $ρambient$⁠) and temperatures (10³-$1.29×108$ K).

At high temperatures of $T>2×107$ 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 × 10⁶ 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×106$ K.

As the 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 unphysical 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 $gNa–Na(r)$ computed in DFT-MD and PIMC simulations using an 8-atom simulation cell at 10⁶ K and 8-fold compression. The good agreement in the $gNa–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.

In dynamic compression experiments, the conservation of mass, momentum, and energy constrains a steady shock to follow the Rankine-Hugoniot equation⁸⁹ 

where (E₀, P₀, V₀) and (E, P, V) represent the internal energy, pressure, and volume of the initial and the shocked states, respectively. The EOS (E and P on a grid of T and V) data in Sec. III B allow 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 $ρambient$⁠. DFT simulations are performed at each of these densities to determine 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×106$ K and the other below 10⁶ K. We tested different bivariate interpolation methods for the EOS grid in $ρ – 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,⁹⁰ 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×107$ 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 the temperature and pressure along the Hugoniot curve as functions of the shock compression ratio $ρ/ρ0$⁠. When the initial density is the lowest (⁠$ρ0=0.75ρambient$⁠), a maximum compression ratio of ∼5.6 is reached at 6 × 10⁵ K. The value of this peak decreases with increasing $ρ0$ and reduces to 4.7 when $ρ0=2.0ρambient$⁠. The higher-temperature compression peak reaches a 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 non-relativistic ideal limit (⁠$ρ/ρ0=4.0$⁠), regardless of the initial density.

Figure 8 compares the principal Hugoniot curve of several elements including He,⁵³ C,⁴⁶ N,⁵⁷ O,⁵⁶ Ne,⁵⁵ and Si.⁴⁷ Interestingly, sodium together with helium shows a 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.⁷⁰ 

Figure 1 compares the principal Hugoniot curve of sodium from our first-principles calculations and those predicted by widely used EOS models, SESAME,⁹⁰ LEOS,^91,92 and average-atom Purgatorio (Lynx).^93,94 The compression maxima along the Hugoniot curves are closely related to ionization and interaction effects. SESAME and LEOS do not explicitly include information about the 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 10⁵ 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^32–34 and theoretical^41,43 data at low temperatures and pressures. In comparison with the experimental data, the local pseudopotential theory⁴¹ underestimates the pressure by up to 30 GPa. The EAM model⁴³ 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 highest-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⁶⁴ 

where $ξ=(P−P0)/ρ0$ and $η=1−ρ0/ρ$⁠. 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<102$ GPa) shock experiments on sodium,^33,34 and a “universal Hugoniot”⁹⁵ obtained from high-pressure (up to 0.2 Gbar) experimental Hugoniot data on Al, Fe, Cu, and Mo. We find that the u_s − u_p relation is well consistent with the linear relations predicted by experiments at $us<15$ km/s^33,34 but increasingly deviates from them when the shock velocities exceeds ∼500 km/s. This high velocity corresponds to a shock pressure of ∼1.9 Gbar along the principal Hugoniot curve, and temperature of $∼107$ K, at which the atoms are nearly fully ionized (see the discussion in Sec. III D). The deviation from the “universal Hugoniot” of fluid metals⁹⁵ 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-principles 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 (⁠$<104$ K) and high (⁠$>107$ K) temperature and pressure conditions (Figs. 6 and 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 10⁷ 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

where σ is the Stefan-Boltzmann constant and c is the speed of light in vacuum. We then re-construct the Hugoniot curves and show them in 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⁷ 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>107$ 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).

The structure of the Hugoniot curves observed in Sec. III C 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_Na–e(r), given by the formula

where the sum includes all electron-ion pairs and θ represents the Heaviside function.

The N_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_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×106$ 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.

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 × 10⁶ to 1 × 10⁶ 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 temperatures of $2×106$ 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×106$ 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 $ρambient$⁠. The value decreases at higher initial densities due to stronger particle interaction, and vice versa. Including a radiation contribution, shock Hugoniots are modified significantly above 10⁷ K and 1 Gbar and a much higher compression ratio can be reached.

See supplementary material for the tables of EOS data.

This research is supported by the U.S. Department of Energy, Grant Nos. 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, the National Center for Atmospheric Research, and the NERSC. This research is part of the Blue Waters sustained-petascale computing project (No. NSF ACI 1640776), which is supported by the National Science Foundation (Award Nos. 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.