*3.2. Elastic Coefficients*

The calculated elastic coefficients and corresponding bulk and shear moduli on the 300, 500, 700, 900, and 1100 K isotherms are listed, respectively, in Tables 1–5. There are two clear trends in the results. For a given temperature, magnitudes of the elastic coefficients increase substantially with increasing pressure. For a given pressure, magnitudes of the elastic coefficients decrease with increasing temperature. Both trends are physically plausible: As the crystal is compressed at a given temperature, it is expected to become more stiff and so the elastic coefficients become larger. On the other hand, as the crystal temperature increases at constant pressure, the crystal expands and in doing so becomes less stiff.

**Table 1.** Pressure-dependent elastic coefficients and isotropic moduli of *β*-HMX at 300 K in units of GPa.


**Table 2.** Pressure-dependent elastic coefficients and isotropic moduli of *β*-HMX at 500 K in units of GPa.



**Table 3.** Pressure-dependent elastic coefficients and isotropic moduli of *β*-HMX at 700 K in units of GPa.




**Table 5.** Pressure-dependent elastic coefficients and isotropic moduli of *β*-HMX at 1100 K in units of GPa.

As an example, Figure 4 shows the pressure dependence of the elastic coefficients and isotropic moduli at 300 K (using data from Table 1). As the pressure changes from 1 atm to 30 GPa, the magnitudes of the elastic coefficients increase by a minimum of about four-fold for *C*<sup>66</sup> to about 60-fold for *C*15. The Voigt–Reuss–Hill bulk and shear moduli increase by about 12- and 6-fold, respectively. Similar behavior is observed at the higher temperatures (not shown).

**Figure 4.** Components of the elastic stiffness tensor at 300 K as functions of hydrostatic pressure. (**a**) Diagonal components and the Voigt–Reuss–Hill bulk modulus. (**b**) Off-diagonal components and the Voigt–Reuss–Hill shear modulus. Lines are added to guide the eye.

The typical temperature dependence of the elastic coefficients is shown in Figure 5, where the results for pressure 5 GPa are presented. This is the lowest pressure among the ones we studied for which *β*-HMX remains crystalline for all temperatures on the time scale considered. The elastic coefficients decrease with increasing temperature approximately linearly. The rate at which the magnitude of the elastic coefficients decreases with temperature ranges from about 2.5 GPa per 100 K for *C*<sup>11</sup> to about 0.07 GPa per 100 K for *C*15. The Voigt–Reuss–Hill bulk and shear moduli decrease with temperature at a rate of about 1 GPa per 100 K and 0.5 GPa per 100 K, respectively. Qualitatively similar temperature dependence of the elastic coefficients was observed at higher pressures but the rates at which the coefficients decrease become smaller.

**Figure 5.** Components of the elastic stiffness tensor at 5 GPa as functions of temperature. (**a**) Diagonal components and the Voigt–Reuss–Hill bulk modulus. (**b**) Off-diagonal components and the Voigt–Reuss–Hill shear modulus. Data for 0 K are taken from Mathew and Sewell [15], who used practically the same force field as here and a similar methodology. Lines are added to guide the eye.

In Figure 6, we compare our results at 1 atm and 300 K to the published theoretical and experimental data for the *β*-HMX elastic coefficients variously at ambient conditions [7–9,14] or zero temperature and pressure [15,33,34]. Our results are, in general, consistent with the MD simulations of Sewell et al. [14] for which practically the same force field was used but with all the covalent bonds frozen. As a result of the bond constraints, most of the elastic coefficients in [14] are slightly larger in magnitude compared to our values. Similarly, there is an overall consistency with the MD/energy-minimization results of Mathew and Sewell [15], obtained with practically the same force field used here but at zero temperature. As can be expected, their values are higher than ours because the crystal stiffens with decreasing temperature, as discussed above. The DFT results in [33,34] give higher values of elastic coefficients and elastic moduli compared to ours. This is consistent with the fact that those results correspond to zero temperature and pressure. Note that, although DFT calculations have been used to calculate lattice parameters and unit-cell volumes on the cold curve (see, e.g., [33]), to our knowledge, they have not been used to predict elastic coefficients at elevated pressures. Moreover, incorporating temperature dependence into such calculations using, for example, the quasi-harmonic approximation would be unreliable due to the need to account for (anisotropic) thermal expansion across hundreds of kelvins; and explicit simulations analogous to, and on the scale of, those studied here are practically infeasible. There is a significant disparity among the various experimental results, which has been attributed to variations in measurement techniques and sample purity. Reanalysis [6,14] of the experimental data has shown that this can also result from lack of redundancy in the acoustic velocity measurements and sensitivity of the numerical solution to initial conditions of the multivariate minimization used to extract the elastic coefficients. For example, it is known that only five of the thirteen nonzero elastic coefficients reported in [7] were accurately

determined [14]. Our results agree reasonably well with the experimental data of Sun et al. [9], with more pronounced differences compared to the other experimental data.

**Figure 6.** Components of the elastic stiffness tensor at 1 atm and 300 K from the present study, previous theoretical results [14,15,33,34], and experimental values [7–9]. (**a**) Diagonal components and the Voigt–Reuss–Hill bulk modulus. (**b**) Off-diagonal components and the Voigt–Reuss–Hill shear modulus. The purple and red arrows emphasize the results from this work and the experimental results due to Sun et al. [9], respectively.

To the best of our knowledge, there are no experimental results for the elastic coefficients of *β*-HMX at higher pressures even for 300 K. However, the room-temperature bulk modulus as a function of pressure can be obtained from fits of experimental isothermal compression data to an equation-of-state fitting form (here, the third-order Birch–Murnaghan equation of state), which yields pressure *P* as a function of volume *V*. The volume-dependent bulk modulus *K*(*V*) can be calculated as *K*(*V*) = −*V*(*∂P*/*∂V*)*T*. The pressure-dependent bulk modulus *K*(*P*) can then be obtained by expressing the volume in *K*(*V*) as a function of pressure using the equation of state. This approach is less precise than obtaining *K* directly from the elastic tensor at a given *T* and *P* due both to the assumed form of the equation of state and the typically small numbers of data points on the isotherm available for fitting. Table 6 lists the room-temperature bulk moduli calculated at five pressures using this approach, for three sets of experimental data [29–31] and our calculated 300 K isotherm in Figure 1. The table also includes the pressure-dependent Voigt–Reuss–Hill bulk moduli reported in Table 1. Note that the results for [30,31] at 10, 20, and 30 GPa represent extrapolations of the fitting form to pressures higher than those for which the corresponding isotherms were measured.



In the table, one can see that the MD-based results (top two rows) are self-consistent in that the bulk moduli computed directly from the elastic tensor and via the equation-of-state fit are in good agreement for all pressures considered. Our values obtained by the direct approach are slightly larger than those obtained using the data from Yoo and Cynn [29] for all pressures except 30 GPa, for which they are about the same. Our estimates based on the data from [30,31] are less reliable because, as noted above, the pressure states considered in those references were below 8 GPa. Nevertheless, there is good agreement between our results and those obtained based on the work of Gump and Peiris [30] with the exception of *P* = 1 atm. Our values are lower than those obtained from the data in [31] for all pressures except *P* = 1 atm, for which our value is larger.
