Equation of state of warm-dense boron nitride combining computation, modeling, and experiment

Equation of state of warm-dense boron nitride combining computation, modeling, and experiment

Shuai Zhang zhang49@llnl.gov    Amy Lazicki jenei2@llnl.gov Lawrence Livermore National Laboratory, Livermore, California 94550, USA    Burkhard Militzer militzer@berkeley.edu Department of Earth and Planetary Science, University of California, Berkeley, California 94720, USA Department of Astronomy, University of California, Berkeley, California 94720, USA    Lin H. Yang    Kyle Caspersen    Jim A. Gaffney    Markus W. Däne    John E. Pask Lawrence Livermore National Laboratory, Livermore, California 94550, USA    Walter R. Johnson Department of Physics, 225 Nieuwland Science Hall, University of Notre Dame, Notre Dame, Indiana 46556, USA.    Abhiraj Sharma    Phanish Suryanarayana College of Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA    Duane D. Johnson Division of Materials Science & Engineering, Ames Laboratory, Ames, Iowa 50011, USA Department of Materials Science & Engineering, Iowa State University, Ames, Iowa 50011, USA    Andrey V. Smirnov Division of Materials Science & Engineering, Ames Laboratory, Ames, Iowa 50011, USA    Philip A. Sterne    David Erskine    Richard A. London    Federica Coppari    Damian Swift    Joseph Nilsen    Art J. Nelson    Heather D. Whitley whitley3@llnl.gov Lawrence Livermore National Laboratory, Livermore, California 94550, USA
July 20, 2019

The equation of state (EOS) of materials at warm dense conditions poses significant challenges to both theory and experiment. We report a combined computational, modeling, and experimental investigation leveraging new theoretical and experimental capabilities to investigate warm-dense boron nitride (BN). The simulation methodologies include path integral Monte Carlo (PIMC), several density functional theory (DFT) molecular dynamics methods [plane-wave pseudopotential, Fermi operator expansion (FOE), and spectral quadrature (SQ)], activity expansion (ACTEX), and all-electron Green’s function Korringa-Kohn-Rostoker (MECCA), and compute the pressure and internal energy of BN over a broad range of densities and temperatures. Our experiments were conducted at the Omega laser facility and the Hugoniot response of BN to unprecedented pressures (1200–2650 GPa). The EOSs computed using different methods cross validate one another in the warm-dense matter regime, and the experimental Hugoniot data are in good agreement with our theoretical predictions. By comparing the EOS results from different methods, we assess that the largest discrepancies between theoretical predictions are 4% in pressure and 3% in energy and occur at  K, slightly below the peak compression that corresponds to the -shell ionization regime. At these conditions, we find remarkable consistency between the EOS from DFT calculations performed on different platforms and using different exchange-correlation functionals and those from PIMC using free-particle nodes. This provides strong evidence for the accuracy of both PIMC and DFT in the high-pressure, high-temperature regime. Moreover, the recently developed SQ and FOE methods produce EOS data that have significantly smaller statistical error bars than PIMC, and so represent significant advances for efficient computation at high temperatures. The shock Hugoniot predicted by PIMC, ACTEX, and MECCA shows a maximum compression ratio of 4.550.05 for an initial density of 2.26 g/cm, higher than the Thomas-Fermi predictions by about 5%. In addition, we construct new tabular EOS models that are consistent with the first-principles simulations and the experimental data. Our findings clarify the ionic and electronic structure of BN over a broad range of temperatures and densities and quantify their roles in the EOS and properties of this material. The tabular models may be utilized for future simulations of laser-driven experiments that include BN as a candidate ablator material. (LLNL-JRNL-767019-DRAFT)

I Introduction

The equation of state (EOS) of materials from the condensed matter to warm dense matter and the plasma regime plays an indispensable role in radiation hydrodynamic simulations Colvin and Larsen (2013), which are required for the design and analysis of inertial confinement fusion (ICF) and high energy density (HED) experiments. In laser-driven capsule experiments, ablator materials are important to implosion dynamics and performance. Currently, the most widely used ablator materials are plastics, such as polystyrene derivatives and glow-discharge polymer, high density carbon (HDC), and beryllium. Materials with higher density and tensile strength, such as boron (B) and its compounds, offer the potential for improvements in performance and additional nuclear diagnostics in exploding pusher platforms. Ellison et al. (2018); Zhang et al. (2018a)

At ambient conditions, BN exists in two stable, nearly degenerate phases: hexagonal BN (h-BN) and cubic BN (c-BN), similar to the graphite and diamond phases of its isoelectronic material, carbon (C). Because of this similarity, BN is widely investigated for the synthesis of superhard materials and fabrication of thin films or heterostructures for various applications. Wang et al. (2017a) Nanostructured c-BN, whose hardness is almost twice that of bulk c-BN and close to that of diamond, has been synthesized at high-pressure and temperature conditions Solozhenko et al. (2012). Other applications for low-dimensional BN include nanoelectronic devices Wang et al. (2017a) and expanded h-BN for hydrogen storage Fu et al. (2017). It has also been demonstrated that the density and mechanical properties of BN can be tuned by constructing a mixture of its cubic and hexagonal phases. Frane et al. (2016)

There have been extensive theoretical and experimental studies on the structure Long et al. (2015); Éric Germaneau et al. (2013), stability BRITUN and Kurdyumov (2000); Sengupta et al. (2018); Taniguchi et al. (1997), EOS Aleksandrov et al. (1989); Wang et al. (2017b); Esler et al. (2010); Kawai et al. (2009); Hu et al. (2018); Marsh (1980), melting and phase diagram Lee et al. (2016); Turkevich (2002); Riedel (1994); de Koker (2012), and mechanical Le Godec et al. (2012); Zhang et al. (2011); Hamdi and Meskini (2010), optical Cunningham et al. (2018); Attaccalite et al. (2011), thermodynamic Wang et al. (2017b, 2009); Yang, W. et al. (2009); Hamdi and Meskini (2010), and transport Chakraborty et al. (2018); Jiang et al. (2018) properties of BN and its polymorphs. The phase transformation of rhombohedral BN (r-BN) was found to be dependent on the pressure transmitting medium Taniguchi et al. (1997), and the transition of h-BN into a wurtzite phase (w-BN) under plastic shear may be dramatically different from that under hydrostatic pressures Ji et al. (2012); Levitas et al. (2004). A large number of calculations using density functional theory (DFT) Kohn and Sham (1965); Hohenberg and Kohn (1964), and quantum Monte Carlo (QMC) simulations Esler et al. (2010); Ma et al. (2015); Atambo et al. (2015) have been performed on c-BN. Assisted by vibrational corrections, QMC results Esler et al. (2010) successfully reproduce the volume changes and Raman frequency shifts measured by static high-pressure experiments.

Experimentally, the diamond anvil cell or multi-anvil apparatus have been used to obtain the EOS of h-BN up to 12 GPa and 1000 K Lynch and Drickamer (1966); Solozhenko et al. (1995); Fuchizaki et al. (2008), c-BN to 160 GPa and 3300 K Datchi et al. (2007); Knittle et al. (1989); Goncharov et al. (2007), and of w-BN to 66 GPa Solozhenko et al. (1998). Shock compression measurements for BN up to 300 GPa have been reported for various initial densities (1.81–3.48 g/cmKawai et al. (2009); Hu et al. (2018); Marsh (1980); Coleburn and Forbes (1968), porosity Marsh (1980), and temperatures (293–713 K) Coleburn and Forbes (1968). Because of the limited data available at extremely high pressure and temperature conditions, existing tabular EOS models have traditionally relied on simplified electronic structure theory, such as the Thomas-Fermi (TF) theory. The goal of this work is to investigate the EOS of BN in the high-energy-density regime and provide new tabular models that are validated by first-principles simulations and experimental data.

In a recent study Zhang et al. (2018a), Zhang et al. computed the EOS of B based on first-principles quantum simulations over a wide range of temperatures and densities. The Hugoniot computed from those simulations shows excellent agreement with our experimental measurement on a planar laser shock platform. We have utilized the data to construct an EOS table (X52) for B. The work has also allowed us to study the performance of the polar direct-drive exploding pusher platform Ellison et al. (2018) and its sensitivity to the EOS.

In this work, we combine extensive theoretical calculations to build tabular models for the EOS of BN, which we then validate in the warm dense matter regime via comparison to experimental measurements of the BN Hugoniot. We also provide theoretical estimates of the uncertainty in the pressure and internal energy by comparing values from different simulation methods. Our theoretical methods include many-body path integral Monte Carlo (PIMC), several electronic structure theories based on pseudopotential DFT-molecular dynamics (DFT-MD), an activity expansion method, and an all-electron, Green’s function Korringa-Kohn-Rostoker (KKR) method. Our experiments consist of three measurements of the Hugoniot response of c-BN conducted at the Omega laser facility.

The paper is organized as follows: Sec. II introduces our simulation methods; Sec. III describes details of our shock experiments; Sec. IV introduces our EOS models; Sec. V compares and discusses our EOS and Hugoniot results from different theoretical methods and experiments and those between BN and C; finally we conclude in Sec. VI.

Ii First-principles simulation methods

Figure 1: Temperature-density diagram showing the parameter regions where the methods in this article are used for calculating the EOS of BN.

In this section, we introduce the theoretical methods that are used in this work to compute the internal energies and pressures of BN across a wide range of temperatures and densities in order to provide simulation data for construction of new tabular EOS models for BN. The theoretical methods applied here include PIMC, the activity expansion method as implemented in the ACTEX code, and several methods that are based on DFT. The DFT methods include both methods that sample the ionic positions via molecular dynamics and average-atom methods where the ionic positions are static. Figure 1 summarizes the temperature and density conditions at which each of the methods has been employed for calculations of BN in this study. In the following, we briefly describe the fundamental assumptions associated with each technique and comment on its accuracy. Additional details can be found in the cited references.

ii.1 Path Integral Monte Carlo

PIMC is a quantum many-body method for materials simulations that is based on sampling the finite temperature density matrix derived from the full many-body Hamiltonian, . In PIMC, particles are treated as quantum paths that are cyclic in imaginary time [0,=], where is the Boltzmann constant. Thermodynamic properties, such as the internal energy, are obtained by


in coordinate representation. is the partition function. is the density matrix. Trotter’s formula Trotter (1959) can be used to break up into slices, each corresponding to an imaginary time step . The method becomes exact in the limit of . Higher temperatures require fewer points, and convergence with respect to the imaginary time step must systematically be tested for each system studied. In practice, one starts with a solution of the two-body problem and only employs the PIMC method to sample higher-order correlations. This pair density matrix approach is described in Refs. Militzer, 2016; Pollock and Ceperley, 1984.

The application of PIMC to electronic structure calculations requires certain approximations due to the fermion sign problem. Fermionic symmetry requires that a negative sign arises from the anti-symmetrical wavefunction. This leads to the nearly complete cancellation of positive and negative contributions to the fermionic density matrix, which makes a direct numerical evaluation impractical for more than a few particles. The standard way to avoid this issue in PIMC simulations is to restrict the paths to the positive region of the trial density matrix, , by implementing the fixed-node approximation Ceperley (1991). The condition in 3-dimensional space defines the nodal surface, where is the number of particles. In high temperature simulations, is chosen to be a Slater determinant of free-particle density matrices


where denotes a plane wave with energy . The corresponding nodal surface is called free-particle nodes. The assumption of free-particle nodes is appropriate at high temperature. The PIMC method with free-particle nodes has been successfully developed and applied to 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 et al. (1999); Militzer and Graham (2006), helium Militzer (2006, 2005), and calculations of the EOS for a range of first-row elements Driver and Militzer (2012); Driver et al. (2015); Driver and Militzer (2016, 2015); Zhang et al. (2018a) and compounds Driver and Militzer (2012, 2017); Zhang et al. (2017a, 2018b). Recent developments Militzer and Driver (2015); Zhang et al. (2017b); Driver et al. (2018) have extended the applicability of PIMC to second-row elements at lower temperatures by appending localized orbitals to , opening a possible route toward accurate quantum many-body simulations of heavier elements.

In this study, we apply PIMC for the simulations of BN with free-particle nodes using the CUPID code Militzer (2000). All electrons and nuclei are treated explicitly as quantum paths. The Coulomb interactions are described via pair density matrices Natoli and Ceperley (1995); Militzer (2016), which are evaluated in steps of  Hartree (Ha). The nodal restriction is enforced in much smaller steps of . The calculations are performed over a wide range of densities 0.23–45.16 g/cm, or 0.1- to 20-times the ambient density  g/cm based on that of h-BN Hassel (1927), and temperatures –510 K. Each simulation cell consists of 24 atoms, which is comparable to our previous simulations for pure B Zhang et al. (2018a), nitrogen (N) Driver and Militzer (2016), and hydrocarbons Zhang et al. (2017a, 2018b). The cell size effects on the EOS are negligible at such high temperature conditions 111By comparing the EOS and the radial distribution function obtained using 24-atom cells to those using 96-atom cells in our DFT-MD calculations, we find negligible differences at temperatures above K. A comparison in is shown in Fig. 8..

ii.2 DFT-MD with plane-wave basis and projector augmented wave potentials

DFT-MD is a widely used method for accurately simulating condensed matter systems at finite temperatures. In DFT-MD, the ions are classical particles, which move according to Newton’s classical equations of motion. The forces are computed by solving the Kohn-Sham DFT equations for the electrons at each time step. The applicability and accuracy of DFT-MD for EOS calculations has been previously demonstrated for condensed phase materials in multiple studies (see Ref. Soderlind and Young, 2018 as an example). One difficulty lies in using this method for high temperatures, which is originated from significant thermal excitation of electrons and intractable computational cost.

Our DFT-MD simulations for BN are performed in two different ways. One way is by using the plane-wave basis and projector augmented wave (PAW) pseudopotentials Blöchl et al. (1994) (pwPAW), as implemented in the Vienna Ab initio Simulation Package (VASPKresse and Furthmüller (1996) and used in our previous studies (e.g, Refs. Zhang et al., 2016, 2017b, 2017a, 2018b, 2018a). Similar to our recent work on pure B Zhang et al. (2018a), we choose the hardest PAW potentials available in VASP, which freeze the 1s electrons in the core and have a core radius of 1.1 Bohr for both B and N. We choose the Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996) functional for describing electronic exchange and correlation interactions, a large cutoff energy of 2000 eV for the plane-wave basis, and the point to sample the Brillouin zone. The simulations are carried out using a Nosé thermostat Nosé (1984) to generate MD trajectories in the canonical ensemble. The MD time step is chosen to ensure total energy conservation and takes on values of 0.05-0.55 fs in these calculations, with smaller values corresponding to higher temperatures. We typically run for 5000 steps at each density-temperature () condition, which is found to be sufficient for convergence of the computed energies and pressures.

To ensure consistency with the all-electron PIMC energies, our pwPAW energies from VASP reported in this study are shifted by -79.017 Ha/BN. This is determined with all-electron calculations for isolated B and N atoms with OPIUM opi () using the PBE functional.

Our pwPAW calculations are performed at temperatures between 6.710 K and 5.05 K (0.6–43.5 eV). Due to limitations in applying the plane-wave expansion for orbitals at low densities and limitations in the applicability of the pseudopotentials that freeze the 1s electrons in the core at high densities, we consider a smaller range of densities ( up to 10) than that was examined via PIMC simulations. These conditions are relevant to shock-compression experiments and span the range in which Kohn-Sham DFT-MD simulations are feasible by conventional wavefunction based approaches. We performed calculations with both 24-atom and 96-atom cells to minimize the finite-size errors.

ii.3 DFT-MD with optimized norm-conserving Vanderbilt pseudopotentials and Fermi-operator expansion

As a check on the pwPAW calculations for the majority of the DFT-MD simulations and to enable extension of our DFT-MD calculations to higher density, we perform a separate set of DFT-MD simulations by utilizing optimized norm-conserving Vanderbilt (ONCV) Hamann (2013, 2017) pseudopotentials—a plane-wave method (pwONCV) at low temperatures and a Fermi operator expansion method (FOE) at high temperatures—in order to verify our pwPAW calculations and expand the range of applicability of Kohn-Sham DFT to higher temperatures. Detailed information about the ONCV pseudopotentials is described in Appendix A.

The pwONCV calculations at low temperature ( K) are similar to those using pwPAW. We applied a preconditioned conjugate gradient method Yang et al. (2007) to fully relax the electronic wavefunctions at each time step. An efficient fast Fourier transform (FFT) algorithm was used for the conversion of the wave functions between real and reciprocal spaces. Each simulation is performed either with frozen 1s core pseudopotentials (for ) or with all-electron pseudopotentials (for ), ensemble with over 5000 steps, time-step of 0.2 fs, and on 128-atom supercells.

At temperatures greater than  K, -shell ionization becomes significant Zhang et al. (2018a). We use all-electron ONCV potentials and FOE Goedecker (1999); Bowler and Miyazaki (2012), which takes advantage of the smooth Fermi-Dirac function at high temperature by approximating the function with polynomial expansion, to conduct Kohn-Sham DFT calculations. In the subspace-projected Hamiltonian approach, we adopted the Chebyshev filtered subspace iteration approach Saad (1992). As the ground-state electron density depends solely on the occupied eigenspace, the technique exploits the fast growth property of Chebyshev polynomial to magnify the relevant spectrum, thereby providing an efficient approach for the solution of the Kohn-Sham eigenvalue problem. The matrix-vector multiplications in the Chebyshev filtering procedure are performed on the FFT grids in Fourier space and only considered if the vector has a non-zero value in the matrix.

Three steps are involved in this method: (i) a Chebyshev filter to construct a subspace which is an approximation to the temperature-smearing occupied eigenspace in a given self-consistent iteration; (ii) FFT mesh to span the Chebyshev filtered subspace from real-space to Fourier space; (iii) FOE in terms of the subspace-projected Hamiltonian represented in the plane-wave basis to compute relevant quantities like the density matrix, electron density and band energy. The accuracy of the Chebychev polynomial expansion Goedecker and Colombo (1994a); Goedecker and Teter (1995a) depends on the electron temperature , and the width of the eigenspectrum . In particular, the degree of polynomial required to achieve the desired accuracy in the approximation Goedecker and Colombo (1994a) of the Fermi-Dirac distribution is . A more accurate estimate that takes into account the location of the Fermi level can be found in Ref. Suryanarayana, 2013. Chebychev polynomial orders of 40–60 and localization radii ranging from 1.056 to 2.88 Bohr were used in the FOE method.

To achieve the same level of accuracy as the plane-wave approach, our high- FOE simulations use PBE exchange-correlation functional and the same FFT meshes as the pwONCV method (real-space grid spacing ranges from 0.066 to 0.18 Bohr). The simulations were carried out using 32-atom supercells. Each simulation involves 3000–6000 steps (0.05–0.1 fs/step) to ensure sufficient statistics.

ii.4 DFT-MD using spectral quadrature

The spectral quadrature (SQ) method Suryanarayana (2013) is a density matrix based method for the solution of the Kohn-Sham equations that is particularly well suited for calculations at high temperature. In the SQ method, all quantities of interest, such as energies, forces, and pressures, are expressed as bilinear forms or sums of bilinear forms which are then approximated by quadrature rules that remain spatially localized by exploiting the locality of electronic interactions in real space Prodan and Kohn (2005), i.e., the exponential decay of the density matrix at finite temperature Goedecker (1998); Ismail-Beigi and Arias (1999); Benzi et al. (2013); Suryanarayana (2017). In the absence of truncation, the method becomes mathematically equivalent to the recursion method Haydock et al. (1972, 1975) with the choice of Gauss quadrature, while for Clenshaw-Curtis quadrature, the FOE Goedecker and Colombo (1994b); Goedecker and Teter (1995b) in Chebyshev polynomials is recovered. Being formulated in terms of the finite-temperature density matrix, the method is applicable to metallic and insulating systems alike, with increasing efficiency at higher temperature as the Fermi operator becomes smoother and density matrix becomes more localized Pratapa et al. (2016); Suryanarayana et al. (2018). scaling is obtained by exploiting the locality of the density matrix at finite temperature, while the exact diagonalization limit is obtained to desired accuracy with increasing quadrature order and localization radius. Convergence to standard planewave results, for metallic and insulating systems alike, is readily obtained Pratapa et al. (2016); Suryanarayana et al. (2018).

While mathematically equivalent to classical FOE methods for a particular choice of quadrature, the more general SQ formulation affords a number of advantages in practice Pratapa et al. (2016); Suryanarayana et al. (2018). These include: (1) The method is expected to be more robust since it explicitly accounts for the effect of truncation on the Chebyshev expansion. (2) The method computes only the elements of density matrix needed to evaluate quantities of interest—e.g., only diagonal elements to obtain densities and energies—rather than computing the full density matrix (to specified threshold) as in FOE methods. (3) The method computes the Fermi energy without storage or recomputation of Chebyshev matrices as required in FOE methods. (4) The method admits a decomposition of the global Hamiltonian into local sub-Hamiltonians in real space, reducing key computations to local sub-Hamiltonian matrix-vector multiplies rather than global full-Hamiltonian matrix-matrix multiplies as in FOE methods. Since the associated local multiplies are small (according to the decay of the density matrix) and independent of one another, the method is particularly well suited to massively parallel implementation; whereas the global sparse matrix-matrix multiplies required in FOE methods pose significant challenges for parallel implementation Bowler and Miyazaki (2012).

In the present work, we employ the massively parallel SQDFT code Suryanarayana et al. (2018) for high-temperature Kohn-Sham calculations. SQDFT implements the SQ method in real space using a high-order finite difference discretization wherein sub-Hamiltonians are computed and applied for each finite-difference grid point. For efficient MD simulations, Gauss quadrature is employed for the calculation of density and energy in each SCF iteration whereas Clenshaw-Curtis quadrature is employed for the calculation of atomic forces and pressure Pratapa et al. (2016). While applicable at any temperature in principle, the present implementation is most advantageous at temperatures in excess of  K, where the Fermi operator becomes sufficiently smooth and density matrix sufficiently localized to reduce wall times below those attainable by standard scaling methods for the system sizes considered here; though avenues exist to reduce this temperature substantially Xu et al. (2018).

Simulations were carried out for a series of 32-atom BN unit cells at densities from 6.77–13.55 g/cm and temperatures from 1010479–1347305 K. All-electron ONCV Hamann (2013) pseudopotentials were employed for B and N with cutoff radii of 0.60 and 0.65 Bohr, respectively. Exchange and correlation were modeled in the local density approximation (LDA) as parametrized by Perdew and Zunger Perdew and Zunger (1981). simulations were carried out using a Nosé-Hoover thermostat Nosé (1984); Hoover (1985) with 500 steps for equilibration followed by 3000–5000 steps for production (with time steps of 0.035–0.04 fs). A finite difference grid spacing of 0.1 Bohr (commensurate with unit cell dimensions), Gauss and Clenshaw-Curtis quadrature orders of 50 and 76, respectively, and localization radius of 1.3 Bohr were employed in the SQ calculations to obtain energies to 0.02% and pressures to 0.2% (discretization error) or less.

ii.5 All-electron, Green’s function Korringa-Kohn-Rostoker

In addition, we applied an all-electron, Green’s function KKR electronic-structure method (based on Kohn-Sham DFT) implemented within a scalar-relativistic approximation, i.e., spin-orbit is ignored beyond the core electrons. We use the Multiple-scattering Electronic-structure Calculation for Complex Applications (MECCA) code, a -space KKR code. Johnson et al. (2015) More technical details on high energy density applications using MECCA and the advantages using a Green function method can be found in reference Wilson et al. (2011). MECCA is applicable to the whole pressure and temperature range of interest in this paper, beyond that available from pseudopotential methods. However, as presently implemented, MECCA is a static DFT code that does not sample the ionic degrees explicitly, i.e., vibrational energies and corresponding entropy contributions cannot be obtained. As such, one must add these either from another calculation or some analytic model. Here, we apply the ideal-gas correction to the MECCA results to provide the most consistent comparisons with the other methods. This approach was used recently to address, for example, the principal Hugoniot curves for Be in a review of EOS models for ICF materials.Gaffney et al. (2018)

For current results, we used the atomic sphere approximation with periodic boundary conditions to incorporate interstitial electron contributions to Coulomb energy from all atomic Voronoi polyhedra. The KKR spherical-harmonic local basis included , i.e., , , and symmetries within the multiple-scattering contributions, and ’s up to 200 are included automatically until the free-electron Bessel functions contribute zero to the single-site wavefunction normalizations. The Green’s functions are integrated via complex-energy contours enclosing a subset of Matsubara poles at finite temperatures, as well as taking advantage of analytic continuation to decrease dramatically solution times. Johnson et al. (1985); Wilson et al. (2011) Various DFT exchange-correlation functionals are included through use of the libXC library. Marques et al. (2012) In this work we used the LDA functional of Vosko, Wilk, and Nusair. Vosko et al. (1980) Brillouin zone integrations for self-consistent charge iterations were performed with a 161616 Monkhorst-Pack Monkhorst and Pack (1976) -point mesh along the complex-energy contour for energies with an imaginary part smaller than 0.25 Rydberg, and a 101010 -point mesh otherwise. A denser mesh was used for the physical density of states calculated along the real-energy axes when needed.

Even though BN occurs in many phases near ambient conditions, for simplicity we chose to use a dense packed but cubic structure, the B2 phase (CsCl prototype) for all MECCA calculation to cover the broad range of pressures and temperatures.

ii.6 Activity expansion

Activity expansion calculations of the EOS are performed using the ACTEX code, which is based on an expansion of the plasma grand partition function in powers of the constituent particle activities (fugacities) Rogers and DeWitt (1973); Rogers (1974). The present calculations are similar to those used in previous work Zhang et al. (2018a) and include interaction terms beyond the Debye-Hückel, electron-ion bound states and ion-core plasma polarization terms, along with relativistic and quantum corrections Rogers (1979, 1981). EOS data generated with the ACTEX code, as well as OPAL opacity tables which use the state populations computed from ACTEX, have been extensively checked by comparison with astronomical observations Rogers and Iglesias (1994) and with laser-driven experiments Rogers and Young (1997).

As with previous studies Zhang et al. (2018a), we cut off ACTEX calculations at temperatures below the point where many-body terms become comparable to the leading-order Saha term ( K). This ensures that the activity expansion method is valid while allowing investigation of the predicted peak compression on the Hugoniot.

Iii Shock Hugoniot experiment

Figure 2: (a) experimental configuration (not drawn to scale), (b) image of a typical c-BN crystal glued to the quartz plate, viewed from the perspective of the VISAR diagnostic and (c) image of the VISAR data from shot 75265, with the analyzed velocities shown as red and blue traces (corresponding the two interferometer legs). The dashed traces are the apparent velocities and the solid traces are corrected for the index of refraction in quartz and cBN.

Experiments to constrain the EOS of BN were performed at the Omega laser facility at the Laboratory for Laser Energetics in Rochester, NY. Samples were c-BN crystals of greater than 99% purity (by weight) and density of 3.45(0.03) g/cm, obtained from Saint-Gobain Ceramic Materials. Pale amber-colored and -oriented (identified by their morphology) optically transparent single crystals were characterized using x-ray photoelectron spectroscopy (XPS) and Raman spectroscopy as in Feng et al. (2013). XPS analysis was performed with a PHI Quantum 2000 system, using focused (11 mm) monochromatic Al x-rays (1486.3 eV). XPS revealed a large amount of C, O and Si contamination, but a 60 second 3 kV Ar ion beam sputter (estimated to remove about 2-5 nm from the surface), dropped the concentration of contaminants by nearly 50%, indicating that these form primarily a surface contamination (a m contaminated surface layer will have no effect on our measurement). After etching, XPS identified a B:N ratio of 1.08:1. Room temperature Raman spectroscopy at 514.5 nm showed the TO and LO phonons of c-BN at 1057.7 and 1309.1 cm, with no sign of the defect bands observed for amber crystals in Ref. Feng et al., 2013, indicating a high bulk purity. An extremely weak peak at 1122.3 cm suggests a negligible contamination of BC.

Crystals with parallel facets separated by 150 m and lateral dimensions of 150-250 m were affixed to 90 m-thick -cut -quartz (density of 2.65 g/cm) windows with micron-scale layers of epoxy. A 3-m thick layer of Au was deposited on the other side of the quartz window, to absorb ablation plasma x-rays and reduce x-ray preheat of the BN samples to negligible levels, and a 25 m-thick layer of plastic was deposited onto the Au to form the laser ablator (Fig. 2(a)).

Samples were ablated directly using 12 beams at of the Omega laser with a 1-ns top-hat pulse shape and distributed phase plates forming a 800 m spot size. Laser energies were tuned to drive the target at intensities ranging from to  TW/cm.

A reflecting shock wave could be tracked continuously as it propagated through the quartz and c-BN samples, using a line-imaging velocimeter (VISAR: Velocity Interferometer System for Any Reflector) Celliers et al. (2004). The in-situ apparent velocities are corrected for the index of refraction of the quartz (1.54687) Ghosh (1999) and c-BN (2.126) Gielisse et al. (1967) at 532 nm, which is the wavelength of the VISAR probe laser.

The shock velocities in the quartz and c-BN at the interface between the two are used in the impedance-matching technique, to determine the EOS data point for c-BN. Because of a finite glue bond thickness between the two materials, the shock velocity in the c-BN must be extrapolated to the quartz surface. The quartz Hugoniot standard is taken from Brygoo et al. (2015) and the reshock model from Knudson and Desjarlais (2013). The shock impedance in cBN at these conditions is higher than quartz, but sufficiently close that the accuracy of the off-Hugoniot quartz model has a small effect on the result (differs by 1% from the result obtained by simply assuming a reflected Hugoniot for the reshock state).

Quartz BN
(km/s) (km/s) (km/s) (GPa) (g/cm)
75265 31.27(0.47) 31.95(0.29) 18.97(0.47) 2091(53) 8.49(0.34)
75263 34.99(0.34) 35.04(0.31) 21.87(0.37) 2643(48) 9.18(0.30)
75264 24.51(0.61) 25.29(0.35) 13.92(0.58) 1214(52) 7.67(0.44)
Table 1: Measured quartz and c-BN shock velocities () and analyzed c-BN particle velocity (), pressure () and density ().

The results of these measurements are recorded in Table 1. Factors contributing to the uncertainty in the Omega measurements include: uncertainty in the quartz and c-BN wave velocities, uncertainty in the extrapolation of the c-BN velocity across the epoxy layer, uncertainty in the initial density of c-BN, and systematic uncertainty in the quartz standard EOS. Uncertainty in the c-BN index of refraction is not quantified so is not included in the error bar.

Iv Construction of EOS models for BN

Before describing the results of the first principles simulations and experiments in detail, we describe the new EOS models and make comparisons to a subset of the calculations. We construct new EOS tables (X2151 and X2152) for BN under the QEOS framework More et al. (1988); Young and Corey (1995). QEOS is a self-contained quasi-single-phase set of thermodynamic models that are widely applicable and guarantee the correct physical limits at both high/low temperature and high/low density. The standard QEOS model based on TF theory also guarantees thermodynamic consistency. In our QEOS framework, we decompose the EOS into separate contributions corresponding to the cold curve, the ion thermal term that describes contributions to the EOS from the ionic degrees of freedom, and the electron thermal term that describes the contributions to the EOS from thermal distribution of the electrons. The cold curve is generally taken from experimental data static DFT calculations, while the electron thermal term is generated using fast electronic structure methods, namely, TF theory and DFT calculations for the average atom-in-jellium model (Purgatorio) described in Appendix B. The ion thermal term is often derived using a form proposed by Cowan More et al. (1988); Young and Corey (1995) and can be modified to fit both experimental data and data from many-body calculations. In condensed phases (at high densities and low temperatures), the EOS, and hence the shock response of materials, is dominated by the cold curve, whereas the ion thermal term dominates the EOS through much of the high-velocity shock regime that is currently accessible in planar experiments at Omega and the National Ignition Facility. The behavior of the EOS and the Hugoniot near peak compression, on the other hand, is mostly dominated by the electron thermal term. The Hugoniot response that a model predicts near peak compression is therefore determined mostly by the underlying electron thermal model, and thus notable differences are seen between TF-based QEOS models and Purgatorio-based QEOS models.

The QEOS framework was chosen due to the lack of data necessary to constrain a more complicated multi-phase EOS representation and because the focus of the current study is in the liquid/plasma region relevant to high velocity, laser-driven shocks. Both X2151 and X2152 tables have reasonably similar parameterization except for the electron-thermal model. At the time when the X2151 table was constructed there was only a Purgatorio 222We performed Purgatorio calculations for boron and for nitrogen in order to generate the electron thermal term of the EOS, which is used for constructing EOS tables based on the QEOS model. Our Purgatorio calculations use the Coulomb potential and Hedin-Lundqvist Hedin and Lundqvist (1971) form of exchange-correlation functional under LDA. electron-thermal model for B, therefore the full electron-thermal model for BN is a mixture of a Purgatorio electron-thermal model for B and a TF electron-thermal model for N. Once a N Purgation electron-thermal model became available, the X2152 table was constructed, where the hybrid TF-Purgatorio electron-thermal model from X2151 was exchanged with a fully Purgatorio electron-thermal model (some adjustments to other EOS parameters were needed to improve the fit for X2152). Therefore, examining the L2150 (legacy TF EOS), X2151, and X2152 gives a demonstration of how the Hugoniot varies from a fully mean-field TF description of ionization, to a hybrid treatment, to a fully quantum atom-in-jellium description.

2.258 g/cm reference density
295 K reference temperature
37 GPa bulk modulus
369 GPa bulk modulus
910 erg/cm cohesive energy
2200 K melt temperature @ 1 bar
1675 K Debye temperature @
1/3 Cowan exponent
Table 2: Key parameters used in the X2152 EOS table.

In both X2151 and X2152, the equilibrium conditions were chosen to be in the hexagonal phase, with a density of 2.258 g/cm, at 295 K and 1 atm. The cold curves are identical in the two models and were fit to calculations from this study and Hugoniot measurements from the Marsh compendium Marsh (1980). Since the ground state phase was taken to be hexagonal the transformation to the cubic phase was represented by employing break-points Young and Corey (1995) to transition from the hexagonal cold-curve to the cubic cold-curve at 10 GPa (the wurtzite phase is essentially combined with the cubic phase in this QEOS form). This transformation pressure is slightly higher than what is reported (1–6 GPa Solozhenko et al. (1999)) but was chosen so that the density where the transformation begins is notably denser than the reference density; this was a practical choice to enhance the stability of the EOS when employed during hydrodynamic simulations. The first-principles isochores calculated for this work were used to constrain the ion-thermal models; specifically, the density dependent Grüneisen model, and the Cowan liquid model. The largest difference between X2151 and X2152 (outside of the electron-thermal model) is that the best ion-thermal fit for X2151 (hybrid electron-thermal) was found using a Cowan exponent of 0.5, conversely the best fit for X2152 (purely Purgatorio) was determined using the canonical value of 1/3. All other EOS parameters (melt temperature, Debye temperature, etc.) were taken directly from known literature. The thermodynamic parameters in the ion thermal model are determined by fitting the pressure data from PIMC, DFT-MD, and ACTEX, taking into account the range of applicability of each method. The key parameters used in X2152 are shown in Table 2. In order to avoid problems with energy offsets (energy zeros) in various techniques, only the pressure data are used for constructing the LEOS tables. The fidelity of this procedure is discussed here.

We note that the EOS obtained using different electronic structure theories can vary depending on the underlying physics. For example, orbital-free (OF) MD, which significantly reduces computational cost of standard DFT-MD by constructing the energy functional in a form that is independent of electronic wavefunctions, predicts CH to be less compressible at the compression maximum than predicted by PIMC and Purgatorio Zhang et al. (2017a, 2018b). Zhang et al. Zhang et al. (2018b) found that this is because the internal energies calculated by OFMD are lower than PIMC, although the pressures are similar, at the same temperatures. Comparing a recent work Danel et al. (2018) on carbon EOS using OFWMD (with W standing for Weizsäcker) to the most recent, Purgatorio-based LEOS 9061 table 333LEOS 9061 is a multi-phase, Purgatorio-based table for carbon developed by fitting the ionic thermodynamic parameters to the first-principles DFT and PIMC data reported in Ref. Benedict et al., 2014, the peak compression predicted by OFWMD is also smaller (4.5 by OFWMD versus 4.6 by LEOS 9061). In addition, OFMD calculations for silicon Hu et al. (2016) shows a single compression maximum along the Hugoniot, whereas PIMC predicts two peaks corresponding to and shell ionization respectively.

Figure 3: (a) Pressure- and (b) temperature-compression Hugoniot of BN predicted by different LEOS models in comparison with PIMC and DFT-MD (pwPAW). The initial density of all Hugoniot curves are set to be 2.15 g/cm. Note that the deviations at above GPa and 2 K are due to the electron relativistic effect, which is included in the Purgatorio tables (thus fully in X2152 and partially in X2151) but not in L2150 or PIMC.

We examine the internal energy differences by comparing the Hugoniot curves for BN based on three LEOS tables (LEOS 2150, X2151, and X2152), for which the electron thermal free energy are constructed differently, as we have explained previously in this section. The results are shown in Fig. 3. Consistent with previous studies, we find that the TF-based model (L2150) predicts a lower peak compression with a broader shape along the vertical axis than the fully Purgatorio-based model (X2152). As expected, the model which combines TF and Purgatorio models lies between the two. Both the shape and the magnitude of the peak compression are intimately related to the -shell ionization of B and N. The TF model is broad due to the neglect of the shell effects, and we observe that the peak compression becomes sharper as one accounts for the -shell ionization of B (X2151), and sharper still when we also account for the shell structure of N (X2152).

Figure 4: Comparison of the pressure and the energy terms of the Hugoniot function along the 2 K isotherm, which is near the compression maximum. Shaded areas denote the error bar of the PIMC data.

The differences in the maximum compression predicted by the different models can be explained by decomposing the Hugoniot function [left-hand side of the Hugoniot equation , where and denote the energy, pressure, and volume of the sample in the shocked and the initial states, respectively] into the energy term and the pressure term and comparing the two as functions of density along isotherms. Figure 4 shows such comparisons at 2 K, which is near the compression maximum along the shock Hugoniot (Fig. 3). The density at which the energy and the pressure curves cross is the Hugoniot density at this temperature. We find that the pressure curves of X2151 and X2152 are on top of each other, but their energies are different. The energies of X2151 are lower, leading to a smaller compression ratio than X2152. In comparison, X2152 data are similar to PIMC in both energy and pressure. This indicates that when constructing an EOS model by merely fitting pressure, it is important to make the electronic contribution fully Purgatorio-based. This is not surprising because Purgatorio is essentially a DFT method. The EOS consistency here demonstrates that the agreement in EOS between PIMC and DFT is not accidental, but represents a consistent description of the electronic interaction in both methods. In addition, Fig. 4 shows the non-smoothness and error bar of the PIMC data at 2 K, which leads to an uncertainty in the compression ratio of 0.05 (or 1%). This represents the level of uncertainty in our reported compression maximum along the Hugoniot by PIMC. At both higher and lower temperatures, the uncertainties are smaller because of the smaller error of the EOS data and higher smoothness of the data along isotherms.

V Results and Discussion

v.1 Isochore Comparisons

Figure 5: Comparison of the pressure-temperature profiles of BN along several isochores from PIMC, DFT-MD (PAW, frozen 1s), DFT-MD (ONCV, frozen 1s), FOE (all-electron), SQ (all-electron), ACTEX, MECCA, and X2152. Subplot (b) is a zoom-in version of (a).

In order to evaluate the performance of recent extensions of DFT methods to high temperature, we compare the computed EOS data from PIMC, pwPAW, pwONCV, FOE, SQ, ACTEX, and MECCA. We choose the X2152 model along several isochores between 0.23 and 45.16 g/cm in Fig. 5 for the basis of performing the comparison. Figure 5(b) highlights the comparison in the temperature range of 10–10 K. This is the regime where electrons are significantly ionized, providing an important testbed for different methods.

Figure 6: Percent difference in internal energy of BN between ACTEX and X2152 along several isochores. The compression ratio (with respect to =2.258 g/cm) are labeled at the top of the plotting area. The reference points for ACTEX and X2152 are both at and ambient temperature.

We find that, at temperatures greater than 2 K, PIMC, ACTEX, and MECCA results show excellent agreement with each other, while the ACTEX predictions are slightly higher than the other two methods only at higher densities. At densities above 4.52 g/cm and temperatures below 1.35 K, deviations of ACTEX from the other methods are evident, which indicates a cut-off temperature () below which the ACTEX method breaks down. This is where the two-body term at order 2 in the activity becomes comparable to the Saha term, which we use as a simple measure of the point where higher order terms start to contribute. Since those terms are not included in ACTEX, we can consider this to be the limit of the current theory. Moreover, we have plotted the percent differences between ACTEX and X2152 data (see Fig. 6 for the comparison in energy; pressure plots look similar), and found the cutoff is dependent on the density: gradually increases from 10 K to 4 K as density increases from 0.1- to 20-times . Above , the agreement between ACTEX and X2152 data is excellent, with differences below 2% in general.

Figure 7: Percent contributions of the ion thermal (left) and electron thermal (right) terms to the total pressure of BN. The remaining contributions are from the cold curve. The temperature-density conditions corresponding to several isochores along which we performed EOS calculations are shown with ’+’ symbols.

Our pressure-temperature profiles by MECCA are overall consistent with those by PIMC, pwPAW, pwONCV, FOE, SQ, and ACTEX. The agreement is best at densities higher than 4.5 g/cm and temperatures higher than  K, where the contributions to the EOS from the ions (the ion thermal contributions) are less significant than those from the thermal electrons (see Fig. 7).

Figure 8: Comparison of the nuclear pair correlation function obtained from DFT-MD (pwPAW) for BN using 24-atom (red) and 96-atom (dark) cells at two different densities and three temperatures. The reference density is 2.26 g/cm. The peaks at K indicates a polymeric structure of the liquid. Differences between small and large cells are evident at 4000 K, indicating a significant finite size effect. This effect is stronger at higher densities and becomes negligible at temperatures higher than 5 K.

At intermediate-low densities (0.23-2.3 g/cm), we observe a discrepancy between MECCA and the DFT-MD/X2152 data, and it grows larger as temperature decreases further below K. This is because the MECCA simulations are performed using static configurations with 2 atoms in the B2 (cesium chloride) structure, which do not include ion motion, and we have thus approximated the ion thermal effect by adding ideal gas corrections to the pressures and energies. However, at the low-temperature conditions, the nuclei show significant correlations by forming polymers, such as N-N pairs or B-N structures that are characterized by the strong fluctuations in the radial pair distribution function at  K and shown in Fig. 8(a)-(c). Therefore, by disregarding the vibrational and rotational contributions, the ideal gas model underestimates the EOS at these conditions. As temperature exceeds  K, the features in the pair distribution function quickly smooth out because the polymeric structures are de-stabilized by thermal effects, which makes the ideal gas approximation for the ions work better and explains the improved agreement between the EOS from DFT-MD and MECCA. Moreover, we note that the agreement between the EOS from X2152 and MECCA can be improved by replacing the ideal-gas correction with the ion thermal model from X2152. The differences at reduce more by applying a constant shift to the MECCA pressures to anchor the pressure-zero point at and 300 K. These findings explain the good consistency between the shock Hugoniot predicted by X2152 and MECCA EOS data, which we address in Sec. V.3.

At densities higher than 2.26 g/cm, the radial distribution function also show significant pair correlations at temperatures below K (Fig. 8(d)-(f)). However, the agreement between the EOS from MECCA and those from DFT-MD are far better than at lower densities. This is the regime where the cold curve contribution dominates the EOS, as Fig. 7 implies. The excellent agreement between MECCA and DFT-MD EOS indicates the effects of the simulation cell and the non-ideal ion thermal contribution are less significant in the more strongly compressed () regime.

At 2.26 g/cm and K, We also observe differences between X2152 and DFT-MD. This can be explained by the differences in the cold curve between X2152 and DFT-MD. The energy minimum in X2152 is set to g/cm corresponding to h-BN, while DFT-MD tends to stabilize c-BN because of the cubic simulation cell being implemented for the liquid simulations. In fact, we found that altering the cold-curve in X2152 such that the is more in line with the ambient density of c-BN allows for better agreement with these low temperature points.

Figure 9: EOS differences of PIMC (red), FOE (black), and MECCA (blue) relative to SQ along two isotherms (1.01 and 1.35 K). Because of the different references chosen in the EOS datasets, all energies have been shifted by the corresponding value at 15.80 g/cm and 1.35 K. The energy differences are normalized by the corresponding ideal gas values ( per BN). The statistical error bars correspond to the 1 uncertainty of the FOE and PIMC data.
(g/cm) (K) (GPa) (Ha/BN) (GPa) (Ha/BN) (GPa) (Ha/BN) (GPa) (Ha/BN)
6.77 1010479 218079 -18.3750.016 2186015 -19.6880.019 21446628 -19.0391.317 21510 -17.589
9.03 1010479 2929712 -20.0450.018 2935518 -21.0060.069 28664775 -20.7911.217 28721 -19.528
10.16 1010479 3313617 -20.6640.021 3321229 -21.5220.040 32070860 -21.9411.201 32411 -20.212
11.29 1010479 3702715 -21.1490.017 3697920 -21.8660.069 36758956 -20.9781.201 36149 -20.783
13.55 1010479 4494623 -21.9210.022 4504040 -22.5110.040 467181087 -19.8571.138 43755 -21.692
15.80 1010479 5317639 -22.3790.032 5331758 -23.4720.046 545621281 -20.6911.152 51570 -22.342
6.77 1347305 3109712 7.5530.020 3076920 5.9130.079 30240577 6.2261.210 30855 8.040
9.03 1347305 4136915 4.5800.019 4129122 3.6340.150 41816759 5.7131.190 41022 5.073
10.16 1347305 4662118 3.5280.022 4665427 2.6130.066 47342858 5.1191.201 46111 3.884
11.29 1347305 5183826 2.5650.029 5190441 2.0570.160 52711964 4.0611.212 51226 2.863
13.55 1347305 6253722 1.1370.021 6263342 0.4150.090 613651153 0.3781.206 61566 1.215
15.80 1347305 7336030 0.0000.024 7358259 0.0000.101 729051299 0.0001.166 72125 0.000
Table 3: Comparison of computed internal energies and pressures from SQ, FOE, PIMC, and MECCA. The energies have been shifted by setting the reference to their respective values at 15.80 g/cm and 1.35 K, at which the pressures are close to each other. The errors in the SQ, FOE, and PIMC data are the statistical 1 error bar determined by blocking analysis Allen and Tildesley (1987).

We compare the EOS data from SQ with those from PIMC, FOE, and MECCA along two different isotherms: 1.01 and 1.35 K. Their values are listed in Tab. 3 and the differences shown in Fig. 9(a) and (b) for pressures and energies, respectively. Our FOE and SQ pressures are in excellent agreement with each other (differences are less than 1%). This can be explained by the use of all-electron ONCV potentials and the DFT-MD nature of both methods. The FOE energies are slightly lower than the SQ values by 1-2% of the corresponding ideal gas values. The small differences can be attributed mainly to different discretization errors in the two approaches, whereas differences associated with trajectory lengths, pseudopotentials, and exchange-correlation functionals were determined to be an order of magnitude smaller.

Our PIMC data at these temperatures scatter around the DFT values, because of the longer paths and larger error bars at such conditions. The differences between PIMC and SQ are 4% in pressure and Ha/atom (or when normalized by the ideal gas value) in energy, which is typical of what we found about differences between PIMC and DFT-MD in previous work on B Zhang et al. (2018a) and hydrocarbon systems Zhang et al. (2018b, a). MECCA data also agree with SQ and FOE at these conditions, with differences in pressure and Ha/atom (or when normalized by corresponding ideal gas values) in energy. The cross validation of the different DFT methods and their consistency with PIMC predictions strongly suggest both the PIMC and the DFT-MD approaches, albeit carrying approximations in each, are reliable for studying the EOS of warm dense matter.

Figure 9 and Table 3 also show the standard error bars of our EOS data, determined by statistical averaging of the MD (for FOE and SQ) or PIMC data blocks. At the temperatures of 1.01–1.35 K, PIMC errors are 2–3% in pressure and 0.6 Ha/atom in energy; FOE errors are 0.05–0.8% in pressure and 0.01–0.08 Ha/atom in energy. In comparison, the statistical error bars of the SQ data are significantly smaller (see Tab. 3). These results, for the first time, establish SQ as an accurate method capable of calculating the EOS of partially ionized, warm-dense plasmas with high precision and accuracy comparable to PIMC.

v.2 Comparison between theory and experiment

Figure 10: Comparison of the Hugoniot of BN from experiment to predictions from PIMC and DFT-MD (pwPAW) simulations and the X2152 model in (a) pressure-density and (b) pressure-compression ratio representations. The initial densities of corresponding Hugoniots are shown in the legend. In (a), equal-temperature conditions along the two Hugoniot curves are connected with lines (as guides to the eyes) to approximate the location of isotherms. The corresponding temperatures are labeled in colored texts. Note that the deviation between PIMC and X2152 curves at above GPa is due to the electron relativistic effect, which is considered in X2152 but not in PIMC.

In this section, we compare our experimental measurements of the pressure-density relation of BN with our theoretical predictions. The experimental data are along the Hugoniot curve, which varies depending on the properties of the sample material. Figure 10 compiles the experimental and theoretical Hugoniot curves corresponding to two different initial densities (): Omega data with of 3.45 g/cm and the Rusbank data Marsh (1980) with of 2.15 g/cm. The corresponding theoretical predictions by X2152 are shown with dark curves. We also show the PIMC and the DFT-MD predictions for 3.45 g/cm and 2.15 g/cm.

The comparison in Fig. 10 shows very good consistency between the measurements and the theoretical predictions. Assisted by the theoretical predictions, we are able to estimate Hugoniot temperatures for the experimental data. We label the Hugoniot temperatures for selected DFT-MD data points with blue-colored text in Fig. 10. We find the Omega data points are in the temperature range of 10–10 K. Our results also show that the PIMC and DFT-MD predicted Hugoniot are in remarkable agreement with X2152 for both initial densities, which spans the Hugoniot curves over a wide range in the phase space. This further shows the validity of the fitting and construction procedure and the quality of our X2152 table. Our calculations and the X2152 model predicts BN to have a maximum compression ratio of 4.59 at 9.8 GPa for  g/cm and 4.47 at 1.8 GPa for  g/cm. We also note that the pressure-density Hugoniots predicted by our different tabular models are very similar (see Fig. 3) at the pressure regime (–3 GPa) explored in our current experiments. We expect future, accurate experiments at higher pressures (e.g., near the compression maximum) to further check our predictions.

v.3 Comparison of different EOS methods

Figure 11: Comparison of the pressure-compression Hugoniot of BN from different theories and LEOS models. The initial density of every Hugoniot curve is 2.26 g/cm. Two sets of DFT-MD (pwPAW) Hugoniots constructed with a difference of the cohesive energy (7.1 eV/atom Ooi et al. (2006)) in the initial energy are also shown for comparison. Note that all MECCA pressures in the EOS have been shifted relative to the value at the initial density and 300 K. Also note that the deviation at above GPa and 2 K is due to the electron relativistic effect, which is considered in X2152 and ACTEX but not in MECCA.

Finally, we make a comprehensive comparison of the shock Hugoniot curves for BN predicted by our different EOS methods. The pressure-compression and temperature-compression Hugoniot curves from ACTEX, TF, MECCA, and X2152 are shown in Fig. 11. We note that ACTEX and X2152 each intrinsically accounts for electron relativistic effects, thus the Hugoniot deviates from the nonrelativistic ideal electron gas limit of 4 at very high temperatures ( K). In comparison, the relativistic correction has not been applied to the TF or MECCA calculations.

At pressures of GPa and temperatures K, ACTEX, X2152, and MECCA yield very similar Hugoniot profiles and a maximum compression of 4.55 for of 2.26 g/cm, while the peak is more broadened according to the TF model and the maximum compression ratio is lower by 0.2. The peak is associated with the shell ionization of B and N, which is smoothed out in the TF model because electronic shell effects are missing in this approach but captured by the other methods. The slightly larger compression predicted by ACTEX than X2152 is consistent with the larger values of the ACTEX EOS data than X2152 (Figs. 5 and  6). The slightly lower compression predicted by MECCA than X2152 can be explained by the non-perfect reconciliation in pressure and energy terms in the Hugoniot function (MECCA pressures are slightly lower while energies are similar in comparison to SQ and PIMC, as shown in Fig. 9).

In the low-temperature condensed matter regime, we find that, with a constant pressure shift in the EOS, our MECCA predictions for the Hugoniot are in good consistency with those of X2152. This indicates the efficacy of using the ideal gas model to approximate the ion thermal effect when constructing EOS using small-size, fixed-lattice models (as in MECCA). Our TF results predict BN to be stiffer in this regime because the initial energy in TF is estimated using an average-atom method (described in Appendix B), which may be higher than the actual value because of the excess energy release due to bonding. We also show differences between X2152 and our DFT-MD (pwPAW) predictions, in particular in Hugoniot temperatures (Fig. 11(b)). This is because of the EOS differences between h-BN and c-BN that we have elaborated previously in Sec. V.1.

v.4 EOS and Hugoniot of isoelectronic materials

Figure 12: (a) Pressure- and (b) temperature-density Hugoniot of BN in comparison with C. The electron thermal contribution to both tables are based on Purgatorio. The initial density of both materials are set to be 2.26 g/cm.

Our EOS models and results for BN enable us to investigate the difference with C—an isoelectronic material of BN. Figure 12 compares the Hugoniot of BN and of C based on X2152 and LEOS 9061, setting their initial densities to be the same (2.26 g/cm). LEOS 9061 is the a multi-phase EOS table constructed for C by using a Purgatorio table for the electron thermal term and fitting DFT and PIMC data Benedict et al. (2014) to obtain the ion thermal term, similar to our present work on BN.

The Hugoniot comparison shows that, at temperature regimes of both  K and K, the compression ratio of BN is higher than C. The compression peak is thus slightly narrower for C. This is because the level of C is in between those of B and N. The differences between BN and C in the low-pressure condensed-matter region ( K) reflect differences in the cold-curve and ion thermal contributions to the EOS. These differences are physically consistent with the influence of different types of interactions between atoms in the two materials. BN has slightly higher ionic character than C due to the differences between the electronegativity of B and N, associated with dipolar interactions between the non-identical atoms.

v.5 Zero-point motion effects

We have also examined the effect of Zero-point motion (ZPM) on the EOS and Hugoniot of BN. In order to do this, we implement the Debye model Wal () to estimate the magnitude of the EOS contributions due to ZPM. This correction reasonably account for the nuclear quantum effects that have been neglected in the our Born-Oppenheimer MD simulations. According to the Debye model, the harmonic vibration energy can be approximated by , where is the volume-dependent Debye temperature and is related to the ambient-density via with being the Grüneisen parameter, and the corresponding pressure . We take the values K and for c-BN from previous measurements and calculations Wang et al. (2017b); de Koker (2012), apply the corrections to our EOS data from DFT-MD (pwPAW) and evaluate the changes in the Hugoniot curve. The results are summarized in Fig. 13.

Figure 13: Zero-point motion effects on the pressure of BN as a function of density along several isotherms. The inset shows the percent increase in pressure for the EOS (black) and along the Hugoniot (red) and percent decrease in compression ratio along the Hugoniot (blue).

Our results show that ZPM causes a pressure increase by over 10% at 6.7 K and ambient density. This percentage difference decreases gradually to at 20 g/cm. The differences dramatically decrease as temperature becomes higher, more so at lower densities. The effect of ZPM on Hugoniot, however, is small. For example, the compression ratio decreases by up to 0.01 (0.4%) for the temperature range 6.7–5.1 K considered in our DFT-MD (pwPAW) simulations. This is similar to what we have seen in carbon-hydrogen systems Zhang et al. (2018b). These findings indicate that the ZPM should be carefully addressed when studying the the low- materials in the condensed matter regime, but is negligible for studying the shock Hugoniot of them in the high-energy-density plasma state.

Vi Conclusions

In this work, we present a comprehensive study of the EOS of BN over a wide range of pressures and temperatures by implementing several computational methods, including PIMC, DFT-MD using standard plane-wave basis and PAW or ONCV potentials, ACTEX, FOE, SQ, MECCA, and TF. We use the PIMC, DFT-MD, and ACTEX data to construct two new EOS tables (X2152 and X2151) for BN using the QEOS model.

Our EOS data by PIMC, FOE, SQ, and MECCA show good consistency at K where 1s electrons are ionized. Our findings establish SQ as an accurate method capable of calculating the EOS with high precision and accuracy comparable to PIMC. Our detailed EOS comparison provides strong evidences that cross validate both the PIMC and the DFT-MD approaches for EOS studies of the partially ionized, warm-dense plasmas.

At 2.5–3.2 K and 1.0–1.3 GPa, our PIMC, ACTEX, and MECCA calculations uniformly predict a maximum compression of 4.55 along the shock Hugoniot for h-BN (=2.26 g/cm), which originates from shell ionization. This compression is underestimated by TF models by 0.2. The maximum compression decreases to 4.47 for c-BN (=3.45 g/cm) and increases to 4.59 for =2.15 g/cm.

We also report Hugoniot data up to GPa from experiments at the Omega laser facility. The measured data show good agreement with our theoretical predictions based on DFT-MD.

By comparing QEOS models with the electron thermal term constructed in different ways (Purgatorio, TF, or hybrid), we find that the shock Hugoniot can be well reproduced by fitting the QEOS models to the pressures in the EOS calculated from first principles. Consistent with our previous studies, we find that the Purgatorio-based EOS models provide the best agreement with both internal energies and pressures from first principles calculations. Because the largest differences in the Hugoniot response of the models occurs near peak compression, performing experiments for materials near peak compression Swift et al. (2012); Kritcher et al. (2016); Nilsen et al. (2016); Swift et al. (2018); Döppner et al. (2018) would provide a rigorous experimental test of our understanding of electronic structure in high energy density plasmas. It would also be worthwhile to pursue experiments that provide measurements of the temperature and the pressure in either Hugoniot or off-Hugoniot experiments, which would provide data to validate the first principle calculations.

We find the shock Hugoniot profiles of isoelectronic materials BN and C are very similar, with the compression peak of C being slightly sharper. This is explained by the differences between the 1s level of C and those of B and N. Based on the similarities of these materials in the laser-induced shock regime, BN ablators would be expected to behave similarly to HDC ablators. While the impact of the condensed phase microstructure of the materials may also be an important consideration in the compressive, ICF regime where much of the ablator is still present during the implosion phase, the microstructure should be less consequential to the behavior of exploding pushers where most of the ablator has been vaporized.

Vii Appendix

vii.1 Optimized norm-conserving Vanderbilt pseudopotentials

We employed ONCV pseudopotentials Hamann (2013) for a subset of DFT-MD calculations, in addition to the FOE and SQ calculations. Fully nonlocal two-projector norm-conserving pseudopotentials were generated. The resulting potentials have an accuracy in electronic structure properties comparable to VASP PAW and all-electron calculations. Due to the wide range of density and temperature grids used in the EOS table generation, we have constructed two versions of ONCV pseudopotentials for B and N to reduce projector overlap and core-state ionization under these extreme conditions. The first set of ONCV pseudopotentials have and valence states for B and and valence states for N, respectively. The second set of ONCV pseudoptentials are all-electron pseudopotentials that include valence. The parameters associated with the corresponding psuedopotentials are listed in Table 4. To cross check the accuracy of the ONCV pseudopotentials we compared calculated pressures with regularized Coulomb potentials ( Bohr and kinetic-energy cutoff of 6000 Ha) for solid c-BN phase at each density-temperature point in the DFT-MD simulations. The overall agreement between ONCV pseudopotentials and regularized Coulomb potentials is within 1% except a few points slightly greater. As an example, Figure 14 shows the percent difference of pressure between all-electron ONCV pseudoptentials and Coulomb potentials for c-BN within the density-temperature grid employed in the DFT-MD simulations. The pressure difference ranges from to , with the larger differences in the low-temperature, low-density regions.

Species Valence Note
(Bohr) (Ha)
B 1.125 35 pwONCV
B 0.6 160 FOE
B 0.6 SQ
N 1.2 35 pwONCV
N 0.65 160 FOE
N 0.65 SQ
Table 4: Parameters used to generate ONCV psuedopotentials for B and N. Bulk properties calculated from these pseudopotentials were benchmark against VASP PAWs and regularized Coulomb potentials. and denote the local potential core radius and the kinetic energy cutoff, respectively. The potentials for SQ are similar to those in FOE, but used higher continuity at to remove cusps and improve convergence.
Figure 14: Percent pressure difference between calculations using ONCV all-electron pseudopotentials and regularized Coulomb potentials for BN in the cubic phase. For most of the phase points examined in this study, the difference is within 1% except a few cases where the difference is slightly greater.

vii.2 Mean-field Thomas-Fermi and average-atom in jellium (Purgatorio)

Our EOS models are developed on a broad grid in phase space, spanning many decades in both temperature and pressure. As such, we require efficient methods for computing the electron thermal contribution to the EOS. In this work, we apply two methods for this purpose, both of which are based on density functional theory. Our TF calculations are based on the generalized theory of Feynman et al. (1949). In contrast to the TF approach, which assumes a uniform Fermi distribution of states and thus does not explicitly include discretized states, Purgatorio solves the electronic structure problem for an atom-in-jellium within LDA self-consistently, and thus allows for the inclusion of discretized states.Wilson et al. (2006); Sterne et al. (2007)

For computing the EOS of mixtures, such as BN, from either Purgatorio or TF, we apply a constant electron pressure mixing rule, following the prescription outlined in Ref. Johnson and Nilsen, 2014. Briefly, if and represent concentrations of the two ions, then the Wigner-Seitz (WS) volume per ion of the plasma is required to be the weighted sum of the WS volumes of its two constituent ions:


In the above, , and are the densities of the plasma and its ionic components, and are atomic weights of the constituent ions and is the Avagadro constant. This equation is supplemented by the requirement that the free electron density of the plasma be unique:


Moreover, since the pressure in the TF theory depends only on and , it follows that the electron density in the plasma is also unique . In the TF method, the free electron density associated with ion is determined by solving the TF equations for the ion at specified values of temperature and density . At a given value of , Eqs. 3-4 provide two equations that can be solved to give values of the unknown densities and . Inasmuch as is a monotonic function of , it follows that the chemical potential is also unique .

To create an EOS table for two-ion plasmas, we first choose a grid uniformly spaced on a logarithmic scale. For each temperature on the grid, we solve the TF equations for the two ions on density sub-grids ranging from 1/2 to 5 times the respective cold-matter densities. The properties of ion 2: , , and , considered as functions of electron density are interpolated onto the electron density grid of ion 1. In this way, Eq. 4 is automatically satisfied at each point on the grid. We can verify that this procedure leads to and for the interpolated values. Furthermore, we can now determine the density of the two-ion plasma at each point on the grid using Eq. 3. In this way, an EOS table is created for as a function of and . The approach is similar for a Purgatorio-based EOS table for a multi-component material: we perform Purgatorio calculations for the individual elements on a (, ) grid and mix the tables according to the pressure equality denoted in Eq. 4.

This work was in part performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. Computational support was provided by LLNL high-performance computing facility (Quartz and Jade) and the Blue Waters sustained-petascale computing project (NSF ACI 1640776). Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. B.M. is supported by the U. S. Department of Energy (grant DE-SC0016248) and the University of California. S.Z. is partially supported by the PLS-Postdoctoral Grant of LLNL. D.D.J. and A.V.S. were partially funded for KKR results by the U. S. Department of Energy, Office of Science, Fusion Energy Sciences through Ames Laboratory, which is operated by Iowa State University for the U.S. DOE under contract DE-AC02-07CH11358. J.P. thanks D.R. Hamann for helpful discussions and code facilitating the construction of robust all-electron ONCV potentials. We appreciate Zsolt Jenei for Raman spectroscopy and John Klepeis, Tadashi Ogitsu and John Castor for useful discussion. This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.


  • Colvin and Larsen (2013) J. Colvin and J. Larsen, Extreme Physics: Properties and Behavior of Matter at Extreme Conditions (Cambridge University Press, 2013).
  • Ellison et al. (2018) C. L. Ellison, H. D. Whitley, C. R. Brown, S. R. Copeland, W. J. Garbett, H. P. Le, M. B. Schneider, Z. B. Walters, H. Chen, J. I. Castor, R. S. Craxton, M. G. Johnson, E. M. Garcia, F. R. Graziani, G. E. Kemp, C. M. Krauland, P. W. McKenty, B. Lahmann, J. E. Pino, M. S. Rubery, H. A. Scott, R. Shepherd,  and H. Sio, Phys. Plasmas 25, 072710 (2018).
  • Zhang et al. (2018a) S. Zhang, B. Militzer, M. C. Gregor, K. Caspersen, L. H. Yang, J. Gaffney, T. Ogitsu, D. Swift, A. Lazicki, D. Erskine, R. A. London, P. M. Celliers, J. Nilsen, P. A. Sterne,  and H. D. Whitley, Phys. Rev. E 98, 023205 (2018a).
  • Wang et al. (2017a) J. Wang, F. Ma, W. Liang,  and M. Sun, Materials Today Physics 2, 6 (2017a).
  • Solozhenko et al. (2012) V. L. Solozhenko, O. O. Kurakevych,  and Y. Le Godec, Advanced Materials 24, 1540 (2012).
  • Fu et al. (2017) P. Fu, J. Wang, R. Jia, S. Bibi, R. I. Eglitis,  and H.-X. Zhang, Computational Materials Science 139, 335 (2017).
  • Frane et al. (2016) W. L. D. Frane, O. Cervantes, G. F. Ellsworth,  and J. D. Kuntz, Diamond and Related Materials 62, 30 (2016).
  • Long et al. (2015) J. Long, C. Shu, L. Yang,  and M. Yang, Journal of Alloys and Compounds 644, 638 (2015).
  • Éric Germaneau et al. (2013) Éric Germaneau, G. Su,  and Q.-R. Zheng, Journal of Physics: Condensed Matter 25, 125504 (2013).
  • BRITUN and Kurdyumov (2000) V. F. BRITUN and A. V. Kurdyumov, High Pressure Research 17, 101 (2000).
  • Sengupta et al. (2018) N. Sengupta, J. E. Bates,  and A. Ruzsinszky, Phys. Rev. B 97, 235136 (2018).
  • Taniguchi et al. (1997) T. Taniguchi, T. Sato, W. Utsumi, T. Kikegawa,  and O. Shimomura, Applied Physics Letters 70, 2392 (1997).
  • Aleksandrov et al. (1989) I. V. Aleksandrov, A. F. Goncharov, S. M. Stishov,  and E. V. Yakovenko, Soviet Journal of Experimental and Theoretical Physics Letters 50, 127 (1989).
  • Wang et al. (2017b) Q. Wang, L. Chen, L. Xiong,  and H. Gong, Journal of Physics and Chemistry of Solids 104, 276 (2017b).
  • Esler et al. (2010) K. P. Esler, R. E. Cohen, B. Militzer, J. Kim, R. J. Needs,  and M. D. Towler, Phys. Rev. Lett. 104, 185702 (2010).
  • Kawai et al. (2009) N. Kawai, M. Yokoo, K.-i. Kondo, T. Taniguchi,  and F. Saito, Journal of Applied Physics 106, 033508 (2009).
  • Hu et al. (2018) X. Hu, G. Yang, B. Zhao, P. Li, J. Yang, C. Leng, H. Liu, H. Huang,  and Y. Fei, Journal of Applied Physics 123, 175903 (2018).
  • Marsh (1980) S. P. Marsh, LASL Shock Hugoniot Data (University of California Press, Berkeley, 1980).
  • Lee et al. (2016) H.-F. Lee, K. Esfarjani, Z. Dong, G. Xiong, A. A. Pelegri,  and S. D. Tse, ACS Nano 10, 10563 (2016), pMID: 27797465.
  • Turkevich (2002) V. Z. Turkevich, Journal of Physics: Condensed Matter 14, 10963 (2002).
  • Riedel (1994) R. Riedel, Advanced Materials 6, 549 (1994).
  • de Koker (2012) N. de Koker, Journal of Physics: Condensed Matter 24, 055401 (2012).
  • Le Godec et al. (2012) Y. Le Godec, O. O. Kurakevych, P. Munsch, G. Garbarino, M. Mezouar,  and V. L. Solozhenko, Journal of Superhard Materials 34, 336 (2012).
  • Zhang et al. (2011) J. S. Zhang, J. D. Bass, T. Taniguchi, A. F. Goncharov, Y.-Y. Chang,  and S. D. Jacobsen, Journal of Applied Physics 109, 063521 (2011).
  • Hamdi and Meskini (2010) I. Hamdi and N. Meskini, Physica B: Condensed Matter 405, 2785 (2010).
  • Cunningham et al. (2018) B. Cunningham, M. Grüning, P. Azarhoosh, D. Pashov,  and M. van Schilfgaarde, Phys. Rev. Materials 2, 034603 (2018).
  • Attaccalite et al. (2011) C. Attaccalite, M. Grüning,  and A. Marini, Phys. Rev. B 84, 245110 (2011).
  • Wang et al. (2009) H. Wang, H. Xu, X. Wang,  and C. Jiang, Physics Letters A 373, 2082 (2009).
  • Yang, W. et al. (2009) Yang, W., Sun, J. X.,  and Yu, F., Eur. Phys. J. B 71, 211 (2009).
  • Chakraborty et al. (2018) P. Chakraborty, G. Xiong, L. Cao,  and Y. Wang, Carbon 139, 85 (2018).
  • Jiang et al. (2018) P. Jiang, X. Qian, R. Yang,  and L. Lindsay, Phys. Rev. Materials 2, 064005 (2018).
  • Ji et al. (2012) C. Ji, V. I. Levitas, H. Zhu, J. Chaudhuri, A. Marathe,  and Y. Ma, Proceedings of the National Academy of Sciences 109, 19108 (2012).
  • Levitas et al. (2004) V. I. Levitas, J. Hashemi,  and Y. Z. Ma, EPL (Europhysics Letters) 68, 550 (2004).
  • Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  • Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  • Ma et al. (2015) F. Ma, W. Purwanto, S. Zhang,  and H. Krakauer, Phys. Rev. Lett. 114, 226401 (2015).
  • Atambo et al. (2015) M. O. Atambo, N. W. Makau, G. O. Amolo,  and R. Maezono, Materials Research Express 2, 105902 (2015).
  • Lynch and Drickamer (1966) R. W. Lynch and H. G. Drickamer, The Journal of Chemical Physics 44, 181 (1966).
  • Solozhenko et al. (1995) V. Solozhenko, G. Will,  and F. Elf, Solid State Communications 96, 1 (1995).
  • Fuchizaki et al. (2008) K. Fuchizaki, T. Nakamichi, H. Saitoh,  and Y. Katayama, Solid State Communications 148, 390 (2008).
  • Datchi et al. (2007) F. Datchi, A. Dewaele, Y. Le Godec,  and P. Loubeyre, Phys. Rev. B 75, 214104 (2007).
  • Knittle et al. (1989) E. Knittle, R. M. Wentzcovitch, R. Jeanloz,  and M. L. Cohen, Nature (London) 337, 349 (1989).
  • Goncharov et al. (2007) A. F. Goncharov, J. C. Crowhurst, J. K. Dewhurst, S. Sharma, C. Sanloup, E. Gregoryanz, N. Guignot,  and M. Mezouar, Phys. Rev. B 75, 224114 (2007).
  • Solozhenko et al. (1998) V. L. Solozhenko, D. Häusermann, M. Mezouar,  and M. Kunz, Applied Physics Letters 72, 1691 (1998).
  • Coleburn and Forbes (1968) N. L. Coleburn and J. W. Forbes, The Journal of Chemical Physics 48, 555 (1968).
  • Trotter (1959) H. Trotter, Proc. Am. Math. Soc. 10, 545 (1959).
  • Militzer (2016) B. Militzer, Comput. Phys. Commun. 204, 88 (2016).
  • Pollock and Ceperley (1984) E. Pollock and D. Ceperley, Phys. Rev. B 30, 2555 (1984).
  • Ceperley (1991) D. M. Ceperley, J. Stat. Phys. 63, 1237 (1991).
  • Pierleoni et al. (1994) C. Pierleoni, D. M. Ceperley, B. Bernu,  and W. R. Magro, Phys. Rev. Lett. 73, 2145 (1994).
  • Magro et al. (1996) W. R. Magro, D. M. Ceperley, C. Pierleoni,  and B. Bernu, Phys. Rev. Lett. 76, 1240 (1996).
  • Militzer and Ceperley (2001) B. Militzer and D. M. Ceperley, Phys. Rev. E 63, 066404 (2001).
  • Militzer et al. (2001) B. Militzer, D. M. Ceperley, J. D. Kress, J. D. Johnson, L. A. Collins,  and S. Mazevet, Phys. Rev. Lett. 87, 275502 (2001).
  • Hu et al. (2010) S. X. Hu, B. Militzer, V. N. Goncharov,  and S. Skupsky, Phys. Rev. Lett. 104, 235003 (2010).
  • Militzer and Ceperley (2000) B. Militzer and D. M. Ceperley, Phys. Rev. Lett. 85, 1890 (2000).
  • Hu et al. (2011) S. X. Hu, B. Militzer, V. N. Goncharov,  and S. Skupsky, Phys. Rev. B 84, 224109 (2011).
  • Militzer et al. (1999) B. Militzer, W. Magro,  and D. Ceperley, Contrib. Plasm. Phys. 39, 151 (1999).
  • Militzer and Graham (2006) B. Militzer and R. L. Graham, J. Phys. Chem. Solids 67, 2136 (2006).
  • Militzer (2006) B. Militzer, Phys. Rev. Lett. 97, 175501 (2006).
  • Militzer (2005) B. Militzer, J. Low Temp. Phys. 139, 739 (2005).
  • Driver and Militzer (2012) K. P. Driver and B. Militzer, Phys. Rev. Lett. 108, 115502 (2012).
  • Driver et al. (2015) K. P. Driver, F. Soubiran, S. Zhang,  and B. Militzer, J. Chem. Phys. 143, 164507 (2015).
  • Driver and Militzer (2016) K. P. Driver and B. Militzer, Phys. Rev. B 93, 064101 (2016).
  • Driver and Militzer (2015) K. P. Driver and B. Militzer, Phys. Rev. B 91, 045103 (2015).
  • Driver and Militzer (2017) K. P. Driver and B. Militzer, Phys. Rev. E 95, 043205 (2017).
  • Zhang et al. (2017a) S. Zhang, K. P. Driver, F. Soubiran,  and B. Militzer, Phys. Rev. E 96, 013204 (2017a).
  • Zhang et al. (2018b) S. Zhang, B. Militzer, L. X. Benedict, F. Soubiran, P. A. Sterne,  and K. P. Driver, J. Chem. Phys. 148, 102318 (2018b).
  • Militzer and Driver (2015) B. Militzer and K. P. Driver, Phys. Rev. Lett. 115, 176403 (2015).
  • Zhang et al. (2017b) S. Zhang, K. P. Driver, F. Soubiran,  and B. Militzer, J. Chem. Phys. 146, 074505 (2017b).
  • Driver et al. (2018) K. P. Driver, F. Soubiran,  and B. Militzer, Phys. Rev. E 97, 063207 (2018).
  • Militzer (2000) B. Militzer, Ph.D. Thesis, University of Illinois at Urbana-Champaign (2000).
  • Natoli and Ceperley (1995) V. Natoli and D. M. Ceperley, J. Comp. Phys. 117, 171 (1995).
  • Hassel (1927) O. Hassel, Norsk Geologisk Tidsskrift 9, 266 (1927).
  • (74) By comparing the EOS and the radial distribution function obtained using 24-atom cells to those using 96-atom cells in our DFT-MD calculations, we find negligible differences at temperatures above K. A comparison in is shown in Fig. 8.
  • Soderlind and Young (2018) P. Soderlind and D. A. Young, Computation 6, 13 (2018).
  • Blöchl et al. (1994) P. E. Blöchl, O. Jepsen,  and O. K. Andersen, Phys. Rev. B 49, 16223 (1994).
  • Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  • Zhang et al. (2016) S. Zhang, K. P. Driver, F. Soubiran,  and B. Militzer, High Energ. Dens. Phys. 21, 16 (2016).
  • Perdew et al. (1996) J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  • Nosé (1984) S. Nosé, J. Chem. Phys. 81, 511 (1984).
  • (81) See http://opium.sourceforge.net for information about the OPIUM code.
  • Hamann (2013) D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  • Hamann (2017) D. R. Hamann, Phys. Rev. B 95, 239906 (2017).
  • Yang et al. (2007) L. H. Yang, R. Q. Hood, J. E. Pask,  and J. E. Klepeis, J. Comput.-Aided Mater. Des. 14, 337 (2007).
  • Goedecker (1999) S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999).
  • Bowler and Miyazaki (2012) D. R. Bowler and T. Miyazaki, Rep. Progr. Phys. 75, 036503 (2012).
  • Saad (1992) Y. Saad, Numerical Methods for Large Eigenvalue Problems (John Wiley, New York, 1992).
  • Goedecker and Colombo (1994a) S. Goedecker and L. Colombo, Phys. Rev. Lett. 73, 122 (1994a).
  • Goedecker and Teter (1995a) S. Goedecker and M. Teter, Phys. Rev. B 51, 9455 (1995a).
  • Suryanarayana (2013) P. Suryanarayana, Chem. Phys. Lett. 584, 182 (2013).
  • Prodan and Kohn (2005) E. Prodan and W. Kohn, Proc. Nat. Acad. Sci. US 102, 11635 (2005).
  • Goedecker (1998) S. Goedecker, Phys. Rev. B 58, 3501 (1998).
  • Ismail-Beigi and Arias (1999) S. Ismail-Beigi and T. A. Arias, Phys. Rev. Lett. 82, 2127 (1999).
  • Benzi et al. (2013) M. Benzi, P. Boito,  and N. Razouk, SIAM Review 55, 3 (2013).
  • Suryanarayana (2017) P. Suryanarayana, Chem. Phys. Lett. 679, 146 (2017).
  • Haydock et al. (1972) R. Haydock, V. Heine,  and M. J. Kelly, J. Phys. C: Solid State 5, 2845 (1972).
  • Haydock et al. (1975) R. Haydock, V. Heine,  and M. J. Kelly, J. Phys. C: Solid State 8, 2591 (1975).
  • Goedecker and Colombo (1994b) S. Goedecker and L. Colombo, Phys. Rev. Lett. 73, 122 (1994b).
  • Goedecker and Teter (1995b) S. Goedecker and M. Teter, Phys. Rev. B 51, 9455 (1995b).
  • Pratapa et al. (2016) P. P. Pratapa, P. Suryanarayana,  and J. E. Pask, Comput. Phys. Commun. 200, 96 (2016).
  • Suryanarayana et al. (2018) P. Suryanarayana, P. P. Pratapa, A. Sharma,  and J. E. Pask, Computer Physics Communications 224, 288 (2018).
  • Xu et al. (2018) Q. Xu, P. Suryanarayana,  and J. E. Pask, J. Chem. Phys. 149, 094104 (2018).
  • Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  • Hoover (1985) W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
  • Johnson et al. (2015) D. Johnson, A. Smirnov,  and S. Khan, MECCA: Multiple-scattering Electronic-structure Calculations for Complex Alloys (KKR-CPA Program, ver. 2.0) (Iowa State University and Ames Laboratory, Ames, 2015).
  • Wilson et al. (2011) B. Wilson, D. Johnson,  and A. Alam, High Energy Density Physics 7, 61 (2011).
  • Gaffney et al. (2018) J. Gaffney, S. Hu, P. Arnault, A. Becker, L. Benedict, T. Boehly, P. Celliers, D. Ceperley, O. ÄŒertík, J. Clérouin, G. Collins, L. Collins, J.-F. Danel, N. Desbiens, M. Dharma-wardana, Y. Ding, A. Fernandez-Pañella, M. Gregor, P. Grabowski, S. Hamel, S. Hansen, L. Harbour, X. He, D. Johnson, W. Kang, V. Karasiev, L. Kazandjian, M. Knudson, T. Ogitsu, C. Pierleoni, R. Piron, R. Redmer, G. Robert, D. Saumon, A. Shamp, T. Sjostrom, A. Smirnov, C. Starrett, P. Sterne, A. Wardlow, H. Whitley, B. Wilson, P. Zhang,  and E. Zurek, High Energy Density Physics 28, 7 (2018).
  • Johnson et al. (1985) D. D. Johnson, F. J. Pinski,  and G. M. Stocks, Journal of Applied Physics 57, 3018 (1985).
  • Marques et al. (2012) M. A. Marques, M. J. Oliveira,  and T. Burnus, Computer Physics Communications 183, 2272 (2012).
  • Vosko et al. (1980) S. H. Vosko, L. Wilk,  and M. Nusair, Canadian Journal of Physics 58, 1200 (1980)https://doi.org/10.1139/p80-159 .
  • Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  • Rogers and DeWitt (1973) F. J. Rogers and H. E. DeWitt, Phys. Rev. A 8, 1061 (1973).
  • Rogers (1974) F. J. Rogers, Phys. Rev. A 10, 2441 (1974).
  • Rogers (1979) F. J. Rogers, Phys. Rev. A 19, 375 (1979).
  • Rogers (1981) F. J. Rogers, Phys. Rev. A 24, 1531 (1981).
  • Rogers and Iglesias (1994) F. J. Rogers and C. A. Iglesias, Science 263, 50 (1994).
  • Rogers and Young (1997) F. J. Rogers and D. A. Young, Phys. Rev. E 56, 5876 (1997).
  • Feng et al. (2013) S. Feng, L. Hou, X. Liu, Y. Gao, X. Li, Q. Wang, Z. Chen, G. Jia,  and J. Zheng, Applied Surface Science 285, 817 (2013).
  • Celliers et al. (2004) P. M. Celliers, D. K. Bradley, G. W. Collins, D. G. Hicks, T. R. Boehly,  and W. J. Armstrong, Rev. Sci. Instrum. 75, 4916 (2004).
  • Ghosh (1999) G. Ghosh, Optics Communications 163, 95 (1999).
  • Gielisse et al. (1967) P. J. Gielisse, S. S. Mitra, J. N. Plendl, R. D. Griffis, L. C. Mansur, R. Marshall,  and E. A. Pascoe, Phys. Rev. 155, 1039 (1967).
  • Brygoo et al. (2015) S. Brygoo, M. Millot, P. Loubeyre, A. E. Lazicki, S. Hamel, T. Qi, P. M. Celliers, F. Coppari, J. H. Eggert, D. E. Fratanduono, D. G. Hicks, J. R. Rygg, R. F. Smith, D. C. Swift, G. W. Collins,  and R. Jeanloz, Journal of Applied Physics 118, 195901 (2015).
  • Knudson and Desjarlais (2013) M. D. Knudson and M. P. Desjarlais, Phys. Rev. B 88, 184107 (2013).
  • More et al. (1988) R. M. More, K. H. Warren, D. A. Young,  and G. B. Zimmerman, Phys. Fluids 31, 3059 (1988).
  • Young and Corey (1995) D. A. Young and E. M. Corey, J. Appl. Phys. 78, 3748 (1995).
  • (126) We performed Purgatorio calculations for boron and for nitrogen in order to generate the electron thermal term of the EOS, which is used for constructing EOS tables based on the QEOS model. Our Purgatorio calculations use the Coulomb potential and Hedin-Lundqvist Hedin and Lundqvist (1971) form of exchange-correlation functional under LDA.
  • Solozhenko et al. (1999) V. L. Solozhenko, V. Z. Turkevich,  and W. B. Holzapfel, The Journal of Physical Chemistry B 103, 2903 (1999)https://doi.org/10.1021/jp984682c .
  • Danel et al. (2018) J.-F. Danel, L. Kazandjian,  and R. Piron, Phys. Rev. E 98, 043204 (2018).
  • (129) LEOS 9061 is a multi-phase, Purgatorio-based table for carbon developed by fitting the ionic thermodynamic parameters to the first-principles DFT and PIMC data reported in Ref. \rev@citealpnumBenedict_2014.
  • Hu et al. (2016) S. X. Hu, B. Militzer, L. A. Collins, K. P. Driver,  and J. D. Kress, Phys. Rev. B 94, 094109 (2016).
  • Allen and Tildesley (1987) M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford Science Publications, 1987).
  • Ooi et al. (2006) N. Ooi, V. Rajan, J. Gottlieb, Y. Catherine,  and J. B. Adams, Modelling and Simulation in Materials Science and Engineering 14, 515 (2006).
  • Benedict et al. (2014) L. X. Benedict, K. P. Driver, S. Hamel, B. Militzer, T. Qi, A. A. Correa, A. Saul,  and E. Schwegler, Phys. Rev. B 89, 224109 (2014).
  • (134) D. Wallace, Statistical Physics of Crystals and Liquids (World Scientific, 2002) pp. 155-217.
  • Swift et al. (2012) D. Swift, J. Hawreliak, D. Braun, A. Kritcher, S. Glenzer, G. W. Collins, S. Rothman, D. Chapman,  and S. Rose, AIP Conf. Proc. 1426, 477 (2012).
  • Kritcher et al. (2016) A. L. Kritcher, T. Doeppner, D. Swift, J. Hawreliak, J. Nilsen, J. Hammer, B. Bachmann, G. Collins, O. Landen, C. Keane, S. Glenzer, S. Rothman, D. Chapman, D. Kraus,  and R. Falcone, J. Phys. Conf. Ser. 688, 012055 (2016).
  • Nilsen et al. (2016) J. Nilsen, B. Bachmann, G. Zimmerman, R. Hatarik, T. Döppner, D. Swift, J. Hawreliak, G. Collins, R. Falcone, S. Glenzer, et al., High Energ. Dens. Phys. 21, 20 (2016).
  • Swift et al. (2018) D. C. Swift, A. L. Kritcher, J. A. Hawreliak, A. Lazicki, A. MacPhee, B. Bachmann, T. Döppner, J. Nilsen, G. W. Collins, S. Glenzer, S. D. Rothman, D. Kraus,  and R. W. Falcone, Review of Scientific Instruments 89, 053505 (2018).
  • Döppner et al. (2018) T. Döppner, D. Swift, A. Kritcher, B. Bachmann, G. Collins, D. Chapman, J. Hawreliak, D. Kraus, J. Nilsen, S. Rothman, L. Benedict, E. Dewald, D. Fratanduono, J. Gaffney, S. Glenzer, S. Hamel, O. Landen, H. Lee, S. LePape, T. Ma, M. MacDonald, A. MacPhee, D. Milathianaki, M. Millot, P. Neumayer, P. Sterne, R. Tommasini,  and R. Falcone, Phys. Rev. Lett. 121, 025001 (2018).
  • Feynman et al. (1949) R. P. Feynman, N. Metropolis,  and E. Teller, Phys. Rev. 75, 1561 (1949).
  • Wilson et al. (2006) B. Wilson, V. Sonnad, P. Sterne,  and W. Isaacs, J. Quant. Spectrosc. Radiat. Transfer 99, 658 (2006).
  • Sterne et al. (2007) P. Sterne, S. Hansen, B. Wilson,  and W. Isaacs, High Energy Density Physics 3, 278 (2007).
  • Johnson and Nilsen (2014) W. Johnson and J. Nilsen, Phys. Rev. E 89, 023107 (2014).
  • Hedin and Lundqvist (1971) L. Hedin and B. I. Lundqvist, J. Phys. C: Solid State Phys. 4, 2064 (1971).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description