Equation of state of warm dense deuterium and its isotopes from density-functional theory molecular dynamics

Of the two approaches of density-functional theory molecular dynamics, quantum molecular dynamics is limited at high temperature by computational cost whereas orbital-free molecular dynamics, based on an approximation of the kinetic electronic free energy, can be implemented in this domain. In the case of deuterium, it is shown how orbital-free molecular dynamics can be regarded as the limit of quantum molecular dynamics at high temperature for the calculation of the equation of state. To this end, accurate quantum molecular dynamics calculations are performed up to 20 eV at mass densities as low as $0.5\mathrm{g}/{\mathrm{cm}}^3$ and up to 10 eV at mass densities as low as $0.2\mathrm{g}/{\mathrm{cm}}^3$. As a result, the limitation in temperature so far attributed to quantum molecular dynamics is overcome and an approach combining quantum and orbital-free molecular dynamics is used to construct an equation of state of deuterium. The thermodynamic domain addressed is that of the fluid phase above 1 eV and $0.2\mathrm{g}/{\mathrm{cm}}^3$. Both pressure and internal energy are calculated as functions of temperature and mass density, and various exchange-correlation contributions are compared. The generalized gradient approximation of the exchange-correlation functional, corrected to approximately include the influence of temperature, is retained and the results obtained are compared to other approaches and to experimental shock data; in parts of the thermodynamic domain addressed, these results significantly differ from those obtained in other first-principles investigations which themselves disagree. The equations of state of hydrogen and tritium above 1 eV and above, respectively, $0.1\mathrm{g}/{\mathrm{cm}}^3$ and $0.3\mathrm{g}/{\mathrm{cm}}^3$, can be simply obtained by mass density scaling from the results found for deuterium. This ab initio approach allows one to consistently cover a very large domain of temperature on the domain of mass density outlined above.

Figure 1

Relative differences between the pressures $P_{{}_{\mathrm{QFCC}}},P_{{}_{\mathrm{VAAQP}}},P_{{}_{\mathrm{INFERNO}}},P_{{}_{\mathrm{OFAA}}}$, and $P_{{}_{\mathrm{OFWAA}}}$, with the latter taken as the reference, at (a) $0.2\mathrm{g}/{\mathrm{cm}}^3$, (b) 2.$5\mathrm{g}/{\mathrm{cm}}^3$, and (c) $20\mathrm{g}/{\mathrm{cm}}^3$ for deuterium. The results have been obtained with the LDA exchange-correlation functional.

Figure 2

Ionic excess pressure of deuterium, calculated with QMD or OFWMD, at the mass densities (a) $0.2\mathrm{g}/{\mathrm{cm}}^3$, (b) $0.5\mathrm{g}/{\mathrm{cm}}^3$, (c) $1\mathrm{g}/{\mathrm{cm}}^3$, (d) 2.$5\mathrm{g}/{\mathrm{cm}}^3$, (e) $8\mathrm{g}/{\mathrm{cm}}^3$, and (f) $20\mathrm{g}/{\mathrm{cm}}^3$. The simulations have been performed with the LDA exchange-correlation functional.

Figure 3

Ratio of ionic excess pressure to total pressure of deuterium at $T=1$, 5, 10, and 20 eV in the domain 0.2–$20\mathrm{g}/{\mathrm{cm}}^3$. The simulations have been performed with QMD and the LDA exchange-correlation functional. The error bars are not displayed when they are smaller than the character size.

Figure 4

Ratio of the ionic internal energy per atom of deuterium, calculated with QMD or OFWMD, to $E_{{}_{iK}}=$ $1.5kT$ at the mass densities (a) $0.2\mathrm{g}/{\mathrm{cm}}^3$, (b) $0.5\mathrm{g}/{\mathrm{cm}}^3$, (c) $1\mathrm{g}/{\mathrm{cm}}^3$, (d) 2.$5\mathrm{g}/{\mathrm{cm}}^3$, (e) $8\mathrm{g}/{\mathrm{cm}}^3$, and (f) $20\mathrm{g}/{\mathrm{cm}}^3$. $k$ designates the Boltzmann constant. The simulations have been performed with the LDA exchange-correlation functional.

Figure 5

Total pressure of deuterium calculated with OFWMD, OFWHMD, and QMD at (a) $1\mathrm{g}/{\mathrm{cm}}^3$, (b) $5\mathrm{g}/{\mathrm{cm}}^3$, and (c) $20\mathrm{g}/{\mathrm{cm}}^3$. The simulations have been performed with the LDA exchange-correlation functional.

Figure 6

Sensitivity of the pressure of deuterium calculated with QOFWH to exchange correlation at (a) $0.2\mathrm{g}/{\mathrm{cm}}^3$, (b) $1\mathrm{g}/{\mathrm{cm}}^3$, and (c) $20\mathrm{g}/{\mathrm{cm}}^3$. The pressure calculated with the LDA exchange-correlation functional is chosen as the reference. In (c), the curves GGAT and LDAT are superimposed.

Figure 7

Comparison of the pressure of deuterium calculated with various approaches along the isochores (a) $0.199561$ $\mathrm{g}/{\mathrm{cm}}^3$, (b) $0.389768$ $\mathrm{g}/{\mathrm{cm}}^3$, (c) $1.15688$ $\mathrm{g}/{\mathrm{cm}}^3$, (d) $1.96361$ $\mathrm{g}/{\mathrm{cm}}^3$, (e) $4.04819$ $\mathrm{g}/{\mathrm{cm}}^3$, and (f) $15.7089$ $\mathrm{g}/{\mathrm{cm}}^3$. The approaches considered are designated as follows: PIMC [1], Caillabet [15], Kerley03 [53], and Holst [14]. QOFWH implemented with the GGAT exchange-correlation functional is taken as the reference.

Figure 8

Hugoniot of deuterium, calculated with the initial conditions ${\rho }_0=0.171$ $\mathrm{g}/{\mathrm{cm}}^3$ and $E_0=-15.886$ eV/at, in the domain $T\ge 1$ eV and $P\le 60$ Mbar. The approaches considered are designated as follows: PIMC2011 [1], PIMC2000 [7], Caillabet [15], and Kerley03 [53]. Our approach, QOFWH, is implemented with the GGAT exchange-correlation functional.

Figure 9

Hugoniot of deuterium, calculated with the initial conditions ${\rho }_0=0.171$ $\mathrm{g}/{\mathrm{cm}}^3$ and $E_0=-15.886$ eV/at, in the domain $T\ge 1$ eV and $P\le 6$ Mbar. The approaches considered are designated as follows: PIMC2011 [1], PIMC2000 [7], Caillabet [15], Lenosky [10], and Kerley03 [53]. Our approach, QOFWH, is implemented with the GGAT exchange-correlation functional. The experimental results found in Refs. [54] and [55] are respectively designated Knudson and Hicks.