Analytic equations of state (EOS) are intended to reproduce theoretical and experimental data in a single phase portion of the thermodynamic space. We devise a complete and thermodynamically consistent model with four distinct features: (1) a reference isotherm that remains thermodynamically stable, (2) a flexible specific heat model based on a fourth-order rational polynomial, (3) a Grüneisen parameter that depends on specific volume and temperature, and (4) pressure and internal energy functions that can be inverted analytically in temperature. The model aims to improve the accuracy of existing equations of state while remaining computationally efficient. To demonstrate its features, we include calibrations for single-crystal pentaerythritol tetranitrate (PETN), liquid nitromethane (NM), and hexagonal close-packed beryllium (Be) metal. The parameter optimization uses the specific heat capacity, Grüneisen parameter, and static compression curves obtained from density functional theory for the crystalline solids and molecular dynamics simulations for liquid NM. We also present a velocity autocorrelation function that yields accurate phonon densities of states for the EOS calibration from the molecular dynamics trajectories. Each of the three calibrations is constrained to enforce the ambient state from experimental measurements and validated against experimental Hugoniot data from multiple sources. We also include one-dimensional hydrodynamic simulations of the isentropic compression experiments for beryllium conducted at the *Z* facility.

## I. INTRODUCTION

Modeling the response of condensed phase materials at high pressure is central to fields such as aerospace, astrophysics, geophysics, and shock and detonation physics. At the continuum level, the governing partial differential equations describe the dynamics of the system where the closure conditions are provided by constitutive models. It is common in high-pressure physics to ignore the contribution of shear deformations to the free energy.^{1} Once the yield condition for plastic flow is exceeded, solids begin to behave like fluids and only a thermodynamically consistent equation of state (EOS) is required.^{2}

Historically, the majority of high-pressure EOS are empirical in nature and rely on limited sets of experimental data to reproduce the relevant physics. Examples of such models include the Davis reactants,^{3} the linear shock-particle velocity, the Hayes,^{4} or the Tillotson^{5} equations of state. Many of these models are of the Mie–Grüneisen type with the full EOS surface based on a reference curve that can be a Hugoniot, an isentrope, or an isotherm. Mie–Grüneisen EOS forms are a simple linear extrapolation off the reference curve. Presuming one samples near the reference curve, an accurate EOS can be attained, but that might not be true in regions far from the reference.

One of the common deficiencies across these models is the use of a reference curve that leads to stability problems under tension.^{6} Solids and liquids can go into expanded states where the material will eventually fail under certain tensile loading conditions (spallation/cavitation). The EOS must remain qualitatively correct to allow for the dynamic evolution of failure models such as Rajendran–Dietenberger–Grove (RDG)^{7} or tensile plasticity (TEPLA).^{8} If the material enters a region where the EOS has unphysical properties—e.g., negative compressibility—the simulation will break down.

A second limitation stems from a lack of thermal data that forces the use of simplifying assumptions—such as a constant specific heat capacity, $ C V $, or a Grüneisen parameter, $\Gamma $, that is only a function of the specific volume, $V$. These approximations are often reasonable for atomic crystals, such as metals, above room temperature.^{9} If the flow regime is dominated by the compressive pressure, any uncertainties in the value of $ C V $ and $\Gamma $ are usually not critical.^{10} Conversely, materials like polymers and energetic molecular crystals are characterized by large molecules with many internal degrees of freedom. In the presence of temperature-dependent chemical reactions, the EOS must be able to capture accurately their thermal response.

Winey *et al.*^{11} recognized these limitations and developed a thermodynamically consistent equation of state for liquid nitromethane (NM) using analytic expressions for $ C V $, the coefficient of thermal pressure, and the isothermal bulk modulus. They report temperature measurements of NM under stepwise loading conditions using Raman spectroscopy. Despite the relatively large uncertainties ( $\xb1100$ K), they show how simplifying assumptions—such as constant $ C V $ or $\Gamma /V$—can lead to significant inaccuracies and stress the need for additional temperature measurements.

The lack of experimental data for constraining an EOS parameterization at very high temperatures and compressions can be mitigated through the use of atomic-scale calculations. Provided an accurate description of the binding energy and interatomic forces for a material is available, the contributions to the Helmholtz free energy from the static lattice energy and the phonon densities of states can be computed using standard methodologies.^{12,13} For crystalline materials, the calculation of the static lattice energy and phonon frequencies for EOS development using density functional theory (DFT)^{14,15} are both routine and accurate given the availability well-validated, high-performance codes. Example equations of state include the ones developed by Cawkwell *et al.*^{16} and Sergeev *et al.*^{17} for pentaerythritol tetranitrate (PETN) and Cawkwell *et al.*^{18} for cyclotetramethylene tetranitramine (HMX).

The calculation of the static lattice energies and vibrational frequencies for systems without long-range order, such as liquids, requires instead molecular dynamics (MD) simulations that typically employ a more approximate description of the interatomic forces than that provided by DFT. Nevertheless, fast semi-empirical, empirical, and machine-learning (ML)-based interatomic potentials and force fields for metals and organic materials are widely available for use in these types of simulations.^{19,20} Given the availability of high-quality calculated data, there is a need for analytic EOS models that can incorporate these advancements in a computationally efficient manner and remain qualitatively correct outside their calibration domain.

The goal is to derive a complete and thermodynamically consistent equation of state that can accurately reproduce the volumetric and thermal response of a broad class of solids and liquids in a single phase domain. The EOS is constructed from a robust cold curve^{6} and a thermal model associated with the lattice vibrations, and optionally, the contributions to the free energy from delocalized electrons in metals. Since the EOS is specifically designed to be utilized in hydrodynamic codes, it must be analytically invertible in pressure, specific internal energy, and temperature. Parametric conditions are also proposed to obtain thermodynamic stability and convexity in the phase space region of interest.

In Sec. II, we present the derivation of the equation of state using a Helmholtz free energy formalism. The model properties and relevant functional forms are proposed in Sec. III. The static DFT and MD calculations used to obtain data for the calibrations of single-crystal pentaerythritol tetranitrate (PETN), liquid nitromethane, and hexagonal close-packed (HCP) beryllium metal are described in Sec. IV. Additionally, we discuss the limitations of the Dickey–Paskin velocity autocorrelation function^{21} for obtaining the phonon DOS and propose an alternative formulation that yields phonon DOS with the appropriate normalization for EOS development. To demonstrate the qualitative features of the EOS, we show calibration examples for PETN, liquid NM, and Be in Sec. V. The three EOS calibrations have been validated against results from high pressure experiments, including single-shock and isentropic compression data, in Sec. VI. Section VII summarizes the main conclusions and the avenues for future work. The supplementary material includes a self-contained section dedicated to code developers interested in the implementation and verification of the equation of state as well as details of the *Z* compression experiment.

## II. HELMHOLTZ FREE ENERGY MODEL

^{22}

^{,}$\Gamma (V,T)$ is determined by $\Gamma (V, T 0 )$ and $ C V (V,T)$ in the general case due to the compatibility condition of Sheffield and Duvall,

^{23}

### A. Cold curve

^{6}to represent the $T=0$ K isotherm

^{6}Equation (7) leads to the following expression for the isothermal bulk modulus:

^{24}In contrast to commonly used isotherms, such as the third-order Birch–Murnaghan

^{25}or the Rose–Vinet,

^{26}the isothermal bulk modulus from Eq. (8) remains positive in all volume space avoiding thermodynamic and numerical stability issues in hydrodynamic simulations.

^{6}The free energy term associated with the potential of the static lattice and the zero point energy can be obtained in closed form from Eq. (7),

### B. Thermal model

^{2}). For insulator-type solids, $ C V \u221e $ can be set to recover the Dulong–Petit limit while a larger value is expected for metals due to the contributions from their delocalized electrons.

^{27,28}For the low temperature asymptotics, the leading term is linear in $T$ except when $a(V)$ vanishes and the Debye-like $ T 3 $ term dominates. Figure 1 depicts the normalized specific heat model, $ C V / C V \u221e $, in Eq. (10) using different $a$ values. The model remains a strictly monotonic function of $\tau $ for $0\u2264a<16/3$. The asymptotic behavior given in Eq. (11) is due to matching the cubic term in the numerator and denominator of the polynomial functions. Such choice also enables the analytic inversion of the internal energy and the pressure functions in terms of the temperature (Sec. III A).

## III. MODEL PROPERTIES

### A. Temperature inversion for internal energy and pressure

### B. Grüneisen parameter

^{1}It can be interpreted as the variation of pressure with thermal energy per unit volume of a material held a constant volume.

^{29}It measures the spacing of the isentropes in the $log\u2061P$– $log\u2061V$ plane and their negative slope in the $log\u2061T$– $log\u2061V$ plane.

^{9}Since the model in Sec. II uses the $ T 0 =0$ K isotherm as reference and $ C V (V, T 0 )=0$, the compatibility condition in Eq. (5) indicates that Grüneisen parameter is solely determined by $ C V $,

^{2}

### C. Temperature scale and linear coefficient

The functional forms in Eq. (25) and (26) with $n\u22650$ guarantee that $\Gamma (V,T)\u22650$ in all phase space. If $n<0$, a separate check should be conducted to verify that $\Gamma $ takes the appropriate values in the $V$– $T$ regions of interest. In Secs. III D and III E, we use thermodynamic stability and convexity at the infinite temperature limit to derive analytic constraints for the fitting parameters.

### D. Thermodynamic stability

^{2}summarizes the conditions for thermodynamic stability in terms of the specific heat capacities and bulk moduli,

#### 1. Atomic crystals with *a* = const.

#### 2. Molecular crystals with Γ_{0} = 0

The bounds in Eqs. (32) and (33) become useful during the EOS calibration but care must be taken to ensure thermodynamic stability in the $V$– $T$ range of interest for the general case. Section V shows a calibration example where $ K T <0$ along the principal isentrope despite satisfying the inequality in Eq. (37). This is generally the case when $ \Gamma 0 \u223c O (1)$ and indicates a deficiency of the model.

### E. Convexity

*fundamental derivative*,

^{30}

^{,}

^{9}This is the expected wave behavior of most materials except in the vicinity of phase transitions.

^{31}Similar to the analysis in Sec. III D, we seek analytic constraints based on the convexity of our EOS model as $T\u2192\u221e$. From Eqs. (28) and (29) and $\Gamma | T \u2192 \u221e \u2261V q 0 $, we can write the following expansion at high temperatures:

*a posteriori*to determine the properties of a particular calibration.

#### 1. Atomic crystals with *a* = const.

## IV. ATOMISTIC SIMULATIONS

^{13}

^{,}

^{32}The latter is ignored for both single-crystal PETN and liquid nitromethane because both are excellent electrical insulators.

### A. Single-crystal pentaerythritol tetranitrate

Under ambient conditions, PETN (C $ 5 $H $ 8 $N $ 4 $O $ 12 $) adopts a tetragonal unit cell with space group $P 4 \xaf 2 1 c$ containing two molecules.^{33} The static compression curve for hydrostatic pressures less than 10 GPa was constructed from the room temperature isothermal compression data measured by Olinger, Cady, and Halleck.^{34} The DFT calculations used to extend the static compression curve to pressures in excess of 10 GPa and to compute the dependence of the vibrational normal modes on specific volume were described in detail in two earlier publications,^{16,35} but here we provide a concise summary for completeness. All of the DFT calculations for PETN were performed using the Gaussian and Plane Waves (GPW) method in the CP2k code^{36} with the generalized gradient approximation (GGA) exchange-correlation functional of Perdew, Burke, and Ernzerhof (PBE),^{37} the double zeta with polarization (DZVP) basis set, Goedecker–Teter–Hutter pseudopotentials,^{38} and the D3(BJ) dispersion correction.^{39} The plane wave and relative cut-off energies were set to 1200 and 60 Ry, respectively. The high pressure portion was calculated in a single unit cell under periodic boundary conditions with a $2\xd72\xd73$ Monkhorst–Pack k-point mesh.^{40} The phonon frequencies were calculated using a brute force (that is, without utilizing crystal symmetry) algorithm to build the force constant matrix via finite displacements of the atoms in a $2\xd72\xd72$ periodic supercell with a single k-point. The normal modes were calculated over compressions $1\u2265V/ V 0 \u22650.6$.

### B. Liquid nitromethane

The static lattice energy and dependencies of the vibrational normal modes on compression of liquid NM (CH $ 3 $NO $ 2 $) were computed using MD simulations rather than the static DFT methods because (i) the traditional force-constant matrix approach is ill-suited to liquids and (ii) the absence of long range order would make any DFT calculations prohibitively expensive. The MD simulations were performed using the LATTE quantum-MD code^{41} that represents the interatomic forces using semi-empirical density functional tight binding (DFTB) theory.^{42,43} We used the *lanl*31 parameter set^{44} for molecules and materials containing C, H, N, and O with the empirical dispersion correction described in Ref. 45 to improve the description of non-bonded interactions between molecules. All of the MD simulations used a time step for the integration of the equation of motion of 0.25 fs. The extended Lagrangian Born–Oppenheimer MD formalism for the propagation of the electronic degrees of freedom during DFTB-MD allowed us to use only one self-consistent field update per MD time step for the calculation of the set of Mulliken partial charges while maintaining no systematic drift in the total energy.^{46,47} The DFTB density matrix was computed using $ O ( N 3 )$ matrix diagonalization at zero electronic temperature.

The MD simulations used to calculate the static lattice energy and vibrational modes were performed in cubic simulation cells that contained 80 NM molecules (560 atoms) under three dimensional periodic boundary conditions. The initial atomic configuration was generated by placing NM molecules in the simulation cell at random positions, while avoiding close approaches and overlaps. The dimensions of this supercell were determined by requiring an initial density of 1.14 g cm $ \u2212 3 $. This configuration was then thermalized at 500 K and constant volume using a Langevin thermostat^{48} for 50 000 time steps to ensure a disordered, liquid-like state was achieved. The volumetric response of liquid NM was examined by generating a series of 61 supercells with compressions $1.05\u2265V/ V 0 \u22650.45$ by uniformly scaling internal coordinates and lattice vectors of the thermalized, ambient density configuration. Following a short geometry optimization to eliminate any highly distorted covalent bonds, MD trajectories were computed for each system in the canonical ensemble at 298 K using a Langevin thermostat. Information from the first 50 000 time steps from each trajectory was discarded. The time histories of the pressure and atomic velocities was recorded for post-processing from the subsequent 100 000 time steps.

The time average of the instantaneous pressures from the set of MD trajectories yields a 298 K isothermal compression curve up to $V/ V 0 =0.45$. A least-squares fit of Eq. (7) yields an equilibrium density of 1.18 g cm $ \u2212 3 $, which is slightly higher than the reported experimental values.^{49} The static compression curve is generated from the 298 K isotherm by removing the contribution to the pressure from $ F ion (V,T)$.

*defined*by Dickey and Paskin

^{21}as

*not*return the number of vibrational modes that are physically associated with those frequencies according to the number and types of bonds present in the molecule. We illustrate the shortcomings in the original Dickey and Paskin form of $\gamma (t)$ below and propose an alternative form for $\gamma (t)$ whose transform to the frequency domain yields a phonon DOS with the appropriate properties under integration.

We start by computing the normal mode frequencies of the nitromethane molecule in the gas phase using the dynamical matrix approach with the *lanl*31 DFTB model. The frequencies of the 21 normal modes are reported in Table I together with values derived from a DFT calculation at the B3LYP/cc-pVTZ level of theory^{53} and experiment.^{54} These results show that the *lanl*31 DFTB model yields normal mode frequencies that are in excellent agreement with experiment and more computationally expensive DFT calculations. These results also highlight that the frequencies and their groupings or degeneracies are connected to the bonding and structure of the molecule. For instance, we observe that the first six modes have imaginary frequencies in the gas phase—these correspond to trivial translations and rotations of the NM molecule. There are three modes with frequencies $ > 3000 c m \u2212 1 $ that correspond to the three C–H stretches in CH $ 3 $NO $ 2 $. The single C–N stretch corresponds to the mode at about $920 c m \u2212 1 $. Indeed, all of the normal modes can be assigned to specific harmonic distortions of the NM molecule.^{54} It is imperative for parameterization of the $ F ion $ term in the EOS that the phonon DOS captures accurately not only the frequencies of the normal modes but also their degeneracies based on the structure and symmetry of the molecule.

. | Frequency (cm^{−1})
. | ||
---|---|---|---|

Mode ID . | B3LYP/cc-pVTZ^{53}
. | lanl31
. | Experiment^{54}
. |

1-6 | … | … | … |

7 | 10 | 1.8 | … |

8 | 481 | 459 | 475 |

9 | 616 | 571 | 603 |

10 | 663 | 626 | 657 |

11 | 926 | 850 | 918 |

12 | 1110 | 1136 | 1096 |

13 | 1137 | 1137 | 1131 |

14 | 1400 | 1357 | 1380 |

15 | 1428 | 1415 | 1397 |

16 | 1465 | 1418 | 1410 |

17 | 1478 | 1444 | 1434 |

18 | 1633 | 1693 | 1583 |

19 | 3077 | 3189 | 2974 |

20 | 3164 | 3277 | 3045 |

21 | 3197 | 3302 | 3080 |

. | Frequency (cm^{−1})
. | ||
---|---|---|---|

Mode ID . | B3LYP/cc-pVTZ^{53}
. | lanl31
. | Experiment^{54}
. |

1-6 | … | … | … |

7 | 10 | 1.8 | … |

8 | 481 | 459 | 475 |

9 | 616 | 571 | 603 |

10 | 663 | 626 | 657 |

11 | 926 | 850 | 918 |

12 | 1110 | 1136 | 1096 |

13 | 1137 | 1137 | 1131 |

14 | 1400 | 1357 | 1380 |

15 | 1428 | 1415 | 1397 |

16 | 1465 | 1418 | 1410 |

17 | 1478 | 1444 | 1434 |

18 | 1633 | 1693 | 1583 |

19 | 3077 | 3189 | 2974 |

20 | 3164 | 3277 | 3045 |

21 | 3197 | 3302 | 3080 |

The cosine transform of Eq. (45) for the condensed-phase MD trajectory computed for the system with density $1.14 g c m \u2212 3 $ is presented in Fig. 4. For clarity, the result has been normalized using Eq. (44) such that $ N mode =21$. The frequencies associated with each peak in the Fig. 4 correspond clearly to those reported in Table I. However, upon integrating this spectrum, which is presented in Fig. 4, we observe that the number of modes associated with each peak is not in accord with the results in presented in Table I. For instance, we observe about six rather than three C–H stretches and four rather than seven modes at very low frequencies. Hence, the cosine transform of Eq. (45) does not correspond to a phonon DOS.

The phonon DOS, $g(f(V))$, for liquid NM obtained using Eq. (46) from our series of DFTB-MD trajectories over the compression ratios $1.05\u2265V/ V 0 \u22650.45$ are presented in Fig. 5. The set of DOS show the shift and broadening of the low frequency peak associated with the harmonic motion of the NM molecules as well as the mode-specific positive shifts in the frequencies of the intramolecular vibrational modes as compression is increased.

### C. Beryllium

The static lattice energy and phonon DOS of HCP Be were calculated using the plane wave DFT code VASP.^{55–58} All calculations used the PBE exchange-correlation functional, a plane wave cut-off energy of 400 eV with the PREC=ACCURATE setting, and the projector augmented wave pseudopotential for Be.^{59,60} The occupancies in the vicinity of the Fermi energy were smeared using the Methfessel–Paxton method with a 0.2 eV smearing width.^{61} The static compression curve was obtained by performing a series of lattice parameter optimizations under hydrostatic stresses between $\u2212$20 and 240 GPa with a $30\xd730\xd720$ $\Gamma $-centered Monkhorst–Pack k-point mesh. The lattice parameters predicted by the static DFT calculations at zero pressure, $a=2.265$ Å, $c=3.564$ Å, and $c/a=1.573$ are slightly smaller than experimental measurements at room temperature, $a=2.287$ Å, $c=3.583$ Å, and $c/a=1.567$.^{62} These results were encouraging because they imply that in the case of Be, the calculated static lattice energy is close to the true one and that a correction by the thermal pressure is not required.

The phonon DOS, $g(f(V))$, of HCP Be were calculated via the dynamical matrix method with a $6\xd76\xd76$ supercell over compression ratios $0.51\u2264V/ V 0 \u22641.0$ with a $5\xd75\xd73$ k-point mesh and are presented in Fig. 6. These calculations were managed using the phonopy package, which uses the space group symmetry of the crystal to construct the full dynamical matrix from a minimal number of atomic displacements.^{63} This reduces the computational burden of building the dynamical matrix from $ O ( N 4 )$ for a brute force algorithm to $ O ( N 3 )$. As compression increases the phonon DOS shift systematically to higher frequencies. The Be pressure-temperature phase diagram shows that HCP Be is thermodynamically stable at low temperature under hydrostatic compression to at least 250 GPa.^{64} The computed phonon DOS are fully consistent with this result because they exhibit no signs of any soft or imaginary phonon modes, which would indicate a structural instability in the HCP phase, at pressures up to 240 GPa.

^{65}the electronic contribution to the free energy is approximated by

*et al.*,

^{66}which yields $ \gamma el , 0 = 2.867 52 \xd7 10 \u2212 5 J g \u2212 1 K \u2212 2 $, $\kappa =0.59$, and $ V ref =0.539 c m 3 g \u2212 1 $ from the static lattice.

## V. CALIBRATION

In this section, we describe the global optimization steps for the thermal and cold components of the model, as well as the different features of each calibration.The devised EOS is intended to capture all three terms of the total free energy in Eq. (41) in the single phase $V$– $T$ range of interest. The calibration is conducted in two sequential steps corresponding to the thermal and cold components of the model. To overcome well-known shortcoming of DFT and MD with regard to its ability to accurately predict the ambient-pressure densities of materials, we constrain the ambient state ( $ P a =0$ GPa, $ T a =298$ K) during the optimization to match the experimental values of density $ \rho a $, thermal expansion coefficient $ \beta a $, and Grüneisen parameter $ \Gamma a $. The values and references for each thermodynamic property are included in Table II. Note that we use the subscript “ $a$” to refer to the ambient state—as opposed to the more common “ $0$”—to avoid confusion with the reference state at $ T 0 =0$ K.

Parameter . | PETN^{34,67}
. | Nitromethane^{49,68}
. | Beryllium^{69}
. |
---|---|---|---|

ρ_{a} (g cm^{-3}) | 1.774 | 1.131 | 1.842 |

β_{a} (K^{−1}) | 23.2 × 10^{−5} | 12.2 × 10^{−4} | 34.5 × 10^{−6} |

Γ_{a} (−) | 1.07 | 1.22 | 1.20 |

### A. Thermal model

^{70}to the entropy along each isotherm before evaluating the volume derivative in Eq. (48). The range and spacing of the temperature points used to generate the $ C V $ and $\Gamma $ surfaces can be chosen based on different criteria. For all three materials, we set an upper limit of $ T max =2000$ K, which approximately corresponds to the Von Neumann spike temperature of detonating PETN and nitromethane, as well as the hexagonal close-packed to body centered cubic (BCC) transition of beryllium at $P\u223c250$ GPa.

^{66}The lower temperature limit is $ T min =0$ K for PETN and nitromethane and $ T min =50$ K for beryllium. The latter is chosen based on where $\Gamma $ derived from the DFT calculations loses accuracy.

All three calibrations are conducted using the *basin-hopping* global optimization algorithm of Wales and Doye^{71} with the implementation included in the LMFIT^{72} Python package. The fits to the normalized specific heat for (a) single-crystal PETN, (b) liquid nitromethane, and (c) HCP beryllium are shown in the top row of Fig. 7 along with their root-mean-square (RMS) deviation. For beryllium, a similar plot for $\Gamma $ is given in Fig. 8. The markers for the thermal part indicate the DFT data plotted every 5 points in $V$ and $T$ and the solid lines represent the EOS model. The colormap represents different isochores ranging from low (blue) to high (red) compression ratios. Note that the results are only used to demonstrate the properties of the equation of state model and do not account for any uncertainties in the DFT or the experimental results.

### B. Cold curve

^{70}is used to obtain $ F st + F zpe $. The parameters $ e 0 $, $ V 0 $, $A$, $B$, and $C$ are then obtained by fitting Eq. (9) to $ F st + F zpe $ in a least squares sense. We impose additional constraints on $ V 0 $ and $A$ during the optimization such that $P( V a , T a )=0$ and $\beta ( V a , T a )= \beta a $ from Table II. Imposing such constraints involves solving the following non-linear root-finding problem for each guess of $ e 0 $, $B$, and $C$:

Parameter . | PETN . | Nitromethane . | Beryllium . |
---|---|---|---|

Thermal model
. | |||

$ C V \u221e ( kJ g \u2212 1 K \u2212 2 )$ [fixed] | 2.2881 × 10^{−3} | 2.5700 × 10^{−3} | 2.9000 × 10^{−3} |

Γ_{0} (−) [fixed] | 1/84 | 1/18 | 5/3 |

Γ_{∞} (−)[fixed] | 2/3 | 2/3 | 2/3 |

θ_{0} (K) | 218.95 | 245.59 | 216.10 |

a_{0} (−) | 1.9477 | 3.7262 | 9.3365 × 10^{−2} |

V_{∞} (cm^{3} g^{−1}) | $ 0.437 52 $ | $ 0.613 27 $ | $ 0.477 78 $ |

m (−) | 3.1674 | 1.9108 | 2.9117 |

n (−) | 3.1674 | 1.9108 | −5.2451 |

Cold curve | |||

e_{0} (kJ g^{−1}) | 1.5700 | 1.9709 | $ 0.966 22 $ |

V_{0} (cm^{3} g^{−1}) | $ 0.539 23 $ | $ 0.744 65 $ | $ 0.540 19 $ |

A (GPa) | $ 0.100 54 $ | 4.5067 × 10^{−2} | 4.6488 |

B (−) | $ 0.666 67 $ | 1.3189 | 1.1014 |

C (−) | 10.466 | 8.4919 | 3.9484 |

Parameter . | PETN . | Nitromethane . | Beryllium . |
---|---|---|---|

Thermal model
. | |||

$ C V \u221e ( kJ g \u2212 1 K \u2212 2 )$ [fixed] | 2.2881 × 10^{−3} | 2.5700 × 10^{−3} | 2.9000 × 10^{−3} |

Γ_{0} (−) [fixed] | 1/84 | 1/18 | 5/3 |

Γ_{∞} (−)[fixed] | 2/3 | 2/3 | 2/3 |

θ_{0} (K) | 218.95 | 245.59 | 216.10 |

a_{0} (−) | 1.9477 | 3.7262 | 9.3365 × 10^{−2} |

V_{∞} (cm^{3} g^{−1}) | $ 0.437 52 $ | $ 0.613 27 $ | $ 0.477 78 $ |

m (−) | 3.1674 | 1.9108 | 2.9117 |

n (−) | 3.1674 | 1.9108 | −5.2451 |

Cold curve | |||

e_{0} (kJ g^{−1}) | 1.5700 | 1.9709 | $ 0.966 22 $ |

V_{0} (cm^{3} g^{−1}) | $ 0.539 23 $ | $ 0.744 65 $ | $ 0.540 19 $ |

A (GPa) | $ 0.100 54 $ | 4.5067 × 10^{−2} | 4.6488 |

B (−) | $ 0.666 67 $ | 1.3189 | 1.1014 |

C (−) | 10.466 | 8.4919 | 3.9484 |

### C. Discussion

As shown in Fig. 9, one of the main improvements of the present EOS is having a Grüneisen parameter that can vary with both temperature and volume. We use the wireframe in Fig. 9 to indicate the $V$– $T$ range where atomistic calculations were used to calibrate the thermal component of the model. For the case of beryllium with $n<0$, the Grüneisen parameter is mostly a function of volume except below room temperature where it decreases. Figure 10 shows that $\Gamma $ becomes negative at sufficiently low temperatures ( $T\u227230$ K), which can lead to an anomalous wave behavior. Thermodynamics places no constraints on the sign of the Grüneisen parameter but $\Gamma <0$—thus $\beta <0$—physically implies that the material contracts upon heating at constant pressure. The entropy surface becomes multivalued in the $P$– $V$ plane when $\Gamma $ changes sign with the Hugoniot and the isentrope becoming tangent where $\Gamma =0$ (see Sec. V C of Menikoff and Plohr^{9}).

A more problematic region for the beryllium calibration occurs where $ K T <0$ at $V/ V a \u22732.5$ expansions along the principal isentrope. This is due to a relatively high $ \Gamma 0 $ value that causes the thermal pressure term of the EOS to increase with volume rather than decrease. Having regions of thermodynamic instability when $ \Gamma 0 \u223c O (1)$ is a deficiency of the model that would require a different functional form for the temperature scale $\theta (V)$. Despite the negative $ K T $ region, the isentropic bulk modulus $ K S $ remains positive in all $V$– $T$, which is needed for the non-thermally conducting compressible fluid flow equations to be hyperbolic. Both PETN and nitromethane calibrations are thermodynamically stable in all $V$– $T$ space and all three calibrations have $ G >0$ everywhere. As shown in Fig. 10, the principal isentrope (blue solid line) expands down to zero temperature for all three materials. In the limiting case of $ \Gamma 0 =0$, the expansion takes place down to a finite temperature because $\Gamma $ controls the slope of the insentropes in $log\u2061T$– $log\u2061V$ space. The principal Hugoniot (red solid line) has a local extremum of $V$ corresponding to the *strong shock limit* at $ V ^ min \u223c 1 2 \Gamma \u221e \Gamma \u221e + 1 $. This limit is physically understood as the balance between shock compression and thermal expansion due to shock heating.^{2}

## VI. DATA COMPARISON

### A. Hugoniots

The top row in Fig. 11 shows the shock velocity $ U s $ vs particle velocity $ u p $ while the bottom row shows pressure vs normalized specific volume. For PETN, the EOS calibration matches the overall data from the Marsh compendium^{74} and Armstrong *et al.*^{75} at medium and high compression, but it appears to worsen at low compression. Due to the lack of data at low particle velocities, we are unable to discern if the discrepancy is due to the EOS model or due to forcing a thermal expansion coefficient $ \beta a $ that is too high. Replacing $ \beta a $ in Table II with the calculated value reported by Perriot^{76} raises the ambient sound speed and improves the agreement with the data (not shown here). If accurate sound speed measurements were to become available, the calibration in Table III could be revisited. For liquid nitromethane, we found good agreement at low and high compression ratios against the data by Marsh^{74} and Lysne and Hardesty.^{68} The data by Lysne and Hardesty^{68} was corrected to $T=298$ K since the original measurements were collected at different initial temperatures. The correction is done using the expressions for the initial sound speed and specific volume provided by the authors. For beryllium, Fig. 11 shows the explosively-driven flyer experiments from Marsh^{74} and the magnetically-driven flyer experiments from McCoy *et al.*^{77} As shown in the phase diagram by Coe *et al.*,^{66} the $ HCP \u2192 BCC $ phase transition on the Hugoniot takes place at $P\u223c150$ GPa. The beryllium calibration in Table III is no longer expected to be accurate at higher pressures, but it remains thermodynamically well-behaved.

### B. Isentropic compression

Swift and co-workers^{73} report isentropic loading results of beryllium using the *Z* pulsed power facility at Sandia National Laboratories. In the experimental setup, a $0.6$ mm-thick aluminum 6061 liner is magnetically driven onto beryllium samples of different thicknesses. Lithium fluoride windows are mounted on the samples and the velocity history of the Be–LiF interface is recorded using laser Doppler velocimetry (point VISAR). The authors collected eight separate velocimetry measurements from *Z shot 843*: four records of the Al liner driving the LiF window without Be, two records for the $0.25$ mm sample of Be, and two records for the $0.50$ mm sample of Be. To assess the EOS model and the Be calibration, we conduct 1D hydrodynamic calculations corresponding to the Al liner driving a $0.25$ mm sample of Be and the LiF window.

The first step is to estimate the time-dependent drive at the aluminum free surface from the four Al-LiF VISAR measurements without beryllium. There are different techniques to unravel the drive but they all require assuming equations of state for the liner and window. We use a Mie–Grüneisen based on a linear $ U s $– $ u p $ and $\rho \Gamma = \rho a \Gamma a $. The aluminum EOS parameters are $ \rho Al =2.703 g c m \u2212 3 $, $ c Al =5.24 mm \mu s \u2212 1 $, $ s Al =1.40$, $ \Gamma Al =1.97$. The LiF EOS parameters are $ \rho LiF =2.638 g c m \u2212 3 $, $ c LiF =5.15 mm \mu s \u2212 1 $, $ s LiF =1.35$, $ \Gamma LiF =1.69$. The pressure at the boundary is obtained iteratively by assuming a $15$-point natural cubic spline and running wave calculations forward in time until the Al-LiF interface matches the experimental velocity record. This is done for all four VISAR records of the Al liner driving the LiF window since they show time differences of approximately $3$ ns. Additional details and the numerical values to reproduce each of the drives can be found in the supplementary material.

Once the time-dependent drive at the free surface is known, we conduct the full calculations with the Al liner, the LiF window, and a $0.25$-mm thick Be sample between the two. The Be is modeled using our EOS and the calibration parameters in Table III. Figure 12 (*left*) shows the velocity histories at the Be–LiF interface from the two VISAR records and the simulations with each of the four pressure drives. Good agreement is found between the experiment and simulations given the time uncertainty in the measurements of the drive. To better understand the dynamics of this problem, Fig. 12 (*middle*) depicts the magnitude of the density gradient corresponding to one of the pressure drives in position–time coordinates. Two acoustic disturbances appear at the Al–Be interface and Be–LiF interface that travel backward to the free surface and then reflect forward to impinge the Be–LiF interface at $ t ( 1 ) \u22480.335\mu s $ and $ t ( 2 ) \u22480.360\mu s $. If we use the particle velocity rather than the pressure at the boundary, the “(1)” disturbance interferes with the drive in a different manner. This is shown in Fig. 12 (*right*) where there is a split in the Be–LiF velocity profile between both types of presumed drives. The “(2)” disturbance can also be seen in both profiles as a small kink in the velocity. Using the pressure as the boundary condition is the common practice and gives better agreement with experiments.

All the calculations shown in Fig. 12 were conducted using a 1D Lagrangian research hydrocode^{78} with the PVRS approximate Riemann solver^{79} and a Courant number of $0.8$. The Be–LiF velocity profiles use an initial spatial resolution of $\Delta x 0 =5\mu m $ and a second-order *minmod* spatial reconstruction that matches the one applied in the forward calculations of the drives. For the density gradient in Fig. 12 (*middle*), we applied a first-order reconstruction instead and an initial resolution of $\Delta x 0 =1\mu m $.

## VII. SUMMARY AND CONCLUSIONS

There is a need for analytic equations of state that can accurately reproduce theoretical and experimental data at a wide range of temperature and specific volume. The devised model is a complete and thermodynamically consistent EOS with four distinct features: (1) a reference isotherm that remains thermodynamically stable, (2) a flexible specific heat model based on a fourth-order rational polynomial, (3) a Grüneisen parameter that depends on $V$ and $T$, and (4) pressure and internal energy functions that are analytically invertible in $T$. The model aims to improve the accuracy of existing equations of state while remaining computationally efficient. To demonstrate the features of the EOS, we include calibrations for single-crystal PETN, liquid nitromethane, and HCP beryllium. The calibrations are conducted in two sequential steps for the thermal model and the cold curve. The optimization uses the specific heat, Grüneisen parameter, and static compression curve obtained from density functional theory while enforcing the ambient state from experimental measurements. To limit the range of values for the fitting parameters, we derive analytic constraints using thermodynamic stability and convexity at the infinite temperature limit. Good agreement was found between the model and experimental Hugoniot data for the three example materials. The agreement is largely tied to the values of the experimental ambient state and the accuracy of the DFT in the range of interest. We also include 1D hydrodynamic simulations of isentropic compression experiments conducted for beryllium at the *Z* facility. The simulations match the velocimetry records within the time uncertainty of the data. Additional calibrations will be conducted in the future for other condensed phase materials—such as polymer-bonded explosives—where a porosity model will be used to augment the current EOS.

## SUPPLEMENTARY MATERIAL

The supplementary material includes a section dedicated to code developers interested in the implementation and verification of the equation of state as well as details of the Z compression experiment.

## ACKNOWLEDGMENTS

The authors would like to thank Ralph Menikoff, Sven Rudin, Romain Perriot, and Jeff Leiding for many insightful discussions. They also thank Damian Swift for sharing the experimental data of the beryllium isentropic compression test and Joseph Zaug for the ultrafast interferometry data of shocked PETN. This work was supported by the Advanced Simulation and Computing (ASC) Program, the Institutional Computing (IC) Program and the Conventional High Explosives (CHE) Grand Challenge project through Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Eduardo Lozano:** Conceptualization (equal); Formal analysis (lead); Methodology (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (lead). **Marc J. Cawkwell:** Data curation (equal); Resources (equal); Validation (equal); Writing – review & editing (equal). **Tariq D. Aslam:** Conceptualization (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material

### APPENDIX A: TEMPERATURE ROOT FROM QUARTIC

### APPENDIX B: DERIVATIVES OF THE TEMPERATURE SCALE AND LINEAR COEFFICIENT

## REFERENCES

*Equations of State for Solids in Geophysics and Ceramic Science*

*ShockWave Science and Technology Reference Library*(Springer, 2007), pp. 143–188.

*Sixth International Symposium on Detonation, Pasadena, CA*(Office of Naval Research, 1976).

*Dynamical Theory of Crystal Lattices*

*Interatomic Forces in Condensed Matter*

*Solid State Physics*

*Fourteenth International Symposium on Detonation, Coeur D’ Alene, Idaho*(Office of Naval Research, 2010).

*Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction*