Thermodynamics of the insulator-metal transition in dense liquid deuterium

Recent dynamic compression experiments [M. D. Knudson et al., Science 348, 1455 (2015); P. M. Celliers et al., Science 361, 677 (2018)] have observed the insulator-metal transition in dense liquid deuterium, but with an approximately 95 GPa difference in the quoted pressures for the transition at comparable estimated temperatures. It was claimed in the latter of these two papers that a very large latent heat effect on the temperature was overlooked in the first, requiring correction of those temperatures downward by a factor of two, thereby putting both experiments on the same theoretical phase boundary and reconciling the pressure discrepancy. We have performed extensive path-integral molecular dynamics calculations with density functional theory to directly calculate the isentropic temperature drop due to latent heat in the insulator-metal transition for dense liquid deuterium and show that this large temperature drop is not consistent with the underlying thermodynamics.


I. INTRODUCTION
The long-standing quest to produce metallic hydrogen in the laboratory has progressed rapidly in the last few years.Recent improvements in static compression techniques, using diamond anvil cells (DAC) coupled with pulsed laser heating, have enabled investigation of the insulator-metal transition both in the cold, dense fluid [1][2][3][4][5] and in the solid 6 .Of the various attempts to achieve metallization in hydrogen, two experimental approaches have used dynamic compression of liquid deuterium to produce multi-megabar pressures 7,8 .Similar in concept, these experiments were designed to probe the metallization of dense liquid deuterium by inducing a small shock followed by nominally isentropic compression.Through this combination of shock and ramp compression, the experiments, by design, avoid the solid phases of deuterium entirely while still maintaining a compression path well below the theoretical critical point temperature on the first-order insulator-metal phase transition boundary.
Experiments on the Sandia Z machine, as reported in Knudson et al. 7 observed indications of a first-order insulator-metal transition in dense liquid deuterium at around 280 to 305 GPa.Subsequent dynamic compression experiments at the National Ignition Facility (NIF), reported in Celliers et al. 8 also found indications of an insulator-metal transition, but with a reported pressure of around 200 GPa.In the Z experiments, the pressure of the transition was deemed to be marked by the rapid drop in reflectivity upon release from a saturated high reflectivity state; it was argued that this provided the clearest signature of the transition, as transients in the system (particularly due to effects of thermal condition) would have ample time to damp out during the several tens of nanoseconds during which the system was in the metallic state.Conversely, in the NIF experiments the transition was denoted to correspond to a reflectivity value of 30% during the initial rise in reflectivity of the compressed deuterium and did not provide data on the release from high pressure.In this context, it is worth noting that the Z experiments indicate both a rise in reflectivity with increasing pressure as well as the abrupt fall with decreasing pressure.In principle, one can argue that the Z experiments also express the full extent of the phase boundary on the return back to low reflectivity.
Direct measurements of the temperature were not made in either of these experiments.Instead, the temperature was estimated using an equation of state (EOS) model for the shock ring-up, followed by estimates of the temperature increase for the quasi-isentropic ramp compression path using either (i) first-principles calculations (Z experiments), or (ii) an average of three tabulated EOS models (NIF experiments).We note that two of these three EOS models are PBE-based global models.It is well known that PBE 9 systematically underestimates the pressure conditions necessary for dissociation, and thus will predict isentropes that exhibit regions of strong −dT /dP (due to latent heat of the transition) at pressures well below the actual metallization boundary; any global model that builds in latent heat well before the transition will underestimate the experimental temperature, perhaps by as much as several hundred K.With this caveat, the resulting estimates of the temperature on encountering the phase boundary were comparable for both sets of experiments, thus presenting ∼95 GPa difference in the location of the metallization boundary at comparable temperatures.
In an effort to reconcile this significant difference in metallization pressure, Celliers et al. argue the NIF experiments probe the entrance of a given isentrope to the coexistence region while the Z experiments probe the exit from the coexistence region.This interpretation, as it stands, is consistent with the difference in the indicators used for identification of the phase boundary in the two sets of experiments: the rise of reflectivity with increasing pressure in the NIF experiments and the drop of reflectivity on the descent from high pressure in the Z experiments.However, it is further argued in Ref. (8) that a factor of two downward correction to the inferred temperatures reported by Knudson et al. 7 is required due to a very large latent heat effect.This reanalysis, formally based on the Clausius-Clapeyron relation, but equivalent to requiring that the two well-separated pressure points (∆P ∼ 95 GPa) lie on the same phase boundary, would require that the temperature estimates in the Z experiments decrease by ∼600 K to nearly 900 K for the lowest and highest temperature loading paths, respectively.Subject to this reinterpretation both sets of experimental data would be in quite good agreement with recent coupled electron-ion Monte Carlo (CEIMC) calculations for the phase boundary reported by Pierleoni et al. 10 , and prior density functional calculations for the vdW-DF1 functional 7 .While this result would be quite appealing from a theoretical point of view, we will show here that under more careful consideration this reinterpretation is in marked disagreement with the thermodynamics of the transition within a first-principles framework.
In particular, the Clausius-Clapeyron analysis presented in Ref. (8), by itself, is incomplete, as it does not include any constraints imposed by the thermodynamics of the phase transition.For a given phase boundary the pressure difference for the entrance and exit of an isentrope, the effective specific heat along the phase boundary, and the latent heat are not arbitrary, but rather dictated by thermodynamics.Here we demonstrate that while the various first-principles frameworks disagree on the location of the first-order transition boundary in (P, T ) space, they are in quite good agreement on the thermodynamics of the transition (i.e. the specific heats, Grüneisen gamma, bulk modulus, etc.).Furthermore, we show that subject to the thermodynamics of the transition, the large temperature drop suggested by Celliers et al. 8 is not thermodynamically consistent; such a large temperature drop would require an anomalously low specific heat for the metallic hydrogen phase.An initial critique of these temperature corrections on thermodynamic grounds was presented in Desjarlais et al. 11 ; counterarguments were presented in Celliers et al. 12 .Here we provide a more detailed analysis, including nuclear quantum effects.

II. PHASE BOUNDARY THERMODYNAMICS
To systematically address the question of the temperature drop resulting from the latent heat of the insulatormetal transition, we start with the assumption, common to both papers cited above, that the experimental path follows an isentrope.Figure 1 is a schematic for the phase boundary and representative isentropes in (P, T ) space.As shown in Fig. 1, the isentropes exhibit a negative slope close to the transition, consistent with a negative Grüneisen γ.See, for example, the isentropes obtained by thermodynamic integration in Fig. 1 of Knudson et al. 7 .For an isentrope S that enters the coexistence region at a point (P I 1 , T 1 ) and a given latent heat ∆H 1 = T 1 ∆S 1 we can readily compute the entropy at Here the superscripts I and M refer to points on the insulator and metallic sides of the coexistence region respectively.To obtain the temperature T 2 at which the isentrope S exits the coexistence region at (P M 2 , T 2 ) we need only compute the temperature drop required to remove the excess entropy ∆S 1 in going from (P M 1 , T 1 ) to (P M 2 , T 2 ).Note that this path is a convenient integration path and should not be confused with the actual experimental thermodynamic path, as was done in Celliers et al. 12 .
We derive in the following the exact expression for change in entropy along a line in (P, T ) space in terms of quantities readily calculated within an N V T firstprinciples framework.We start with the total derivative of S(P, T ): We can specify an arbitrary line in (P, T ) space by imposing the linear constraint This results in where we have used the definition where C P is the specific heat at constant pressure.We can further reduce this in terms of quantities that are readily calculable with the Maxwell relation Given the definitions where B T is the isothermal bulk modulus, C V is the specific heat at constant volume, and γ is the Grüneisen γ, we arrive at With calculations of C V , γ, and B T , C P is readily obtained via In Celliers et al. 12 this expression ( 7) was noted to be equivalent to Eq. (4.19) in Reichl 13 , but with the erroneous interpretation that this expression is only to be applied at constant volume.As is clear from the derivation, this expression for the change in entropy is completely general for a given line in (P, T ) space.For the specific application here, this expression is applied along a line adjacent to, but on the metallic side of the phase transition, as indicated by the red arrow in Fig. 1.For the slope of the coexistence boundary in the vicinity of 1300K to 1400 K, we use dP/dT coex = -0.12GPa/K as an average, which agrees well with the published coexistence line for vdW-DF1, including nuclear quantum effects, around that temperature range 7 .C V , γ, and B T are determined from first-principles calculations.

III. FIRST-PRINCIPLES CALCULATIONS WITH NUCLEAR QUANTUM EFFECTS
To obtain quantitative values for C V , γ, and B T , including nuclear quantum effects, we have performed several path integral molecular dynamics (PIMD) calculations with finite-temperature density functional theory as implemented in the VASP 5.3.5 [14][15][16] , using the PIMD scheme of Alfè and Gillan 17 .The DFT-PIMD calculations were performed with the vdW-DF1 18 exchangecorrelation functional.Note the choice of vdW-DF1 for these calculations is motivated by arguments in Celliers et al. 8 that (i ) the onset of the phase transition is in close agreement with the phase boundary predicted by either CEIMC or DFT calculations with vdW-DF1 and nuclear quantum effects, and (ii ) that the Z experiments are indicative of the exit from that boundary.Therefore vdW-DF1 is a logical choice for directly addressing the thermodynamics of their argument in that region of phase space in a density functional framework.However, FIG. 2. Energy versus temperature in the vicinity of the insulator-metal transition at specific volumes of 1.80 (green circles, dotted fit) and 1.78 (black circles, solid fit) Å3 /atom.
CV is given by the slope.
as discussed in Knudson et al. 7 the phase boundary suggested by the Z experiments is in better agreement with the vdW-DF2 functional and the effect of the latent heat was implicitly computed in Knudson et al. 7 through direct calculation of the isentropes on both sides of the phase boundary by thermodynamic integration.Each DFT-PIMD calculation consists of 8 path integral molecular dynamics images for a Trotter time step ≤ 9.5 × 10 −5 K −1 , with each image containing 256 deuterium atoms represented with a PAW potential 19,20 .The Brillouin zone was sampled at the Baldereschi mean value point 21 and the plane wave cutoff energy was 700 eV.Individual DFT-PIMD runs consisted of 4000 time steps of 0.25 fs, for a total simulated time of 1 ps.An Andersen thermostat 22 was employed to regulate the temperature and approximate a canonical ensemble.Due to the relatively short nature of each individual run, we observed statistical variation in the temperature about the target temperature, along with expected statistical variation in the thermodynamic quantities.To minimize the statistical variance in the thermodynamic averages, generalized virial estimators 23 were used for the pressure and energy.
The results of these calculations for specific volumes of 1.80 Å3 /atom (green circles, dotted fits) and 1.78 Å3 /atom (black circles, solid fits) are illustrated in Fig. 2 as the total energy versus temperature, and in Fig. 3 as pressure × volume versus the energy, providing C V and γ, respectively.The lower temperature bound of the data set for each fit was chosen to correspond with a sharp drop in the dimer peak of the pair correlation function, as computed from the path centroid positions and indicative of the molecular to atomic transition.By this construction we are performing calculations adjacent to the metallic side of the phase transition in the region traversed by the path from (P M 1 , T 1 ) to (P M 2 , T 2 ) suggested by the red arrow in Fig. 1.For a specific volume of 1.80 Å3 /atom we find a transition at ∼ 1345 K and P = 204 GPa.For a specific volume of 1.78 Å3 /atom we To determine B T adjacent to the phase boundary, we have computed P versus T over a wide range of volumes and temperatures as shown in Fig. 4. From the individual fits to the pressure, we generate P versus V at 1345 K, as shown by the corresponding second order polynomial fit in Fig. 5. Differencing the pressure fit at V = 1.80 Å3 /atom and using the second definition in Eq. ( 6) yields B T = 272 GPa.We find that B T increases rapidly over the course of several 10s of GPa, reaching 540 GPa within 50 GPa at 1345 K, however it is this lower local value that is relevant to the calculation along the phase boundary and which slightly depresses the effective specific heat along the coexistence line. From Å3 /atom, we find C P = 6.07 k B /atom.Combining terms and inserting in Eq. ( 7), we find with C coex = 2.64 k B /atom.Integrating Eq. ( 8) with the constraint 2 1 dS = −∆S 1 , and assuming C coex is weakly varying or represents an average value, gives

IV. ISENTROPIC TEMPERATURE DROP DUE TO LATENT HEAT
There are previous calculations of the latent heat for this insulator-metal transition in the literature 10,24 .We find, consistent with the Pierleoni et al. results for quantum protons, a latent heat of ∼0.05 ± 0.005 eV/atom (or 522 K to 638 K), over the temperature range of 1300−1350 K. Using ∆H 1 = 0.05 eV/atom for the latent heat, and T 1 = 1345 K, sets ∆S 1 =0.431 k B /atom.For C coex = 2.64 k B /atom as calculated above at V = 1.80 Å3 /atom, T 2 /T 1 = 0.85, according to Eq. ( 9).We do not find these calculations very sensitive to the 0.050 ± 0.005 eV/atom extrema quoted for the latent heats, producing only a ∓1.5% change in T 2 /T 1 .
These final results are only very weakly dependent on the explicit T dependence within the square brackets of Eq. ( 8), so for the purposes of computations above we have treated T as a constant equal to T 1 .Because of the weak dependence, iteratively replacing T by T = (T 1 + T 2 )/2 converges very rapidly.As a specific quantitive example of this, consider a case where the cal-culated C coex = 2.64 k B /atom at V = 1.80 Å3 /atom is assumed to represent a midpoint or average value C coex within the temperatures spanned from T 1 to T 2 , where T 2 /T 1 = exp −∆S 1 /C coex .Iterating the above equations as described, with the target T = (T 1 + T 2 )/2 = 1345 K, results in a solution with T 1 = 1447 K, ∆S 1 = 0.401 k B /atom, T 2 = 1243 K, and T 2 /T 1 = 0.86.For a phase boundary slope of dP/dT coex = -0.12GPa/K, this calculated ∆T = 202 K along the coexistence boundary would correspond to ∆P = 24 GPa.Comparing these results to the green curve in Fig. 1 of Knudson et al. 7 reveals very good agreement with the latent heat induced temperature drop obtained there by thermodynamic integration using the vdW-dF2 functional 25 .Furthermore, Figure 4c (green curve) of Knudson et al. 7 , which shows reflectivity as a function of pressure and corresponds to these same approximate conditions, is in good agreement with this ∆P estimate on the fast rising linear portion of the reflectivity history.
It is interesting to note that despite the relatively high value of C P = 6.07 k B /atom found on the metallic side of the phase boundary, the negative slope of the phase boundary and the negative slope of the isentrope near the boundary (γ < 0) combine to give an effective specific heat through Eq. ( 9) that is very close to what we find well away from the phase boundary in the region where γ ∼ 0. In that region, at a specific volume of 1.60 Å3 /atom and 1345 K (around 254 GPa in Fig. 4), the specific heat C P ≈ C V = 2.61 k B /atom.As was noted in Desjarlais et al. 11 , values of specific heat of this magnitude are expected for liquid alkali metals; C P = 3.5 k B /atom for liquid lithium metal 26 .

V. DISCUSSION
Given the close agreement between (i ) the results presented here using vdW-DF1, (ii ) the estimates provided in Desjarlais et al. 11 obtained using interpolations on CEIMC data from Pierleoni et al. 10 , and (iii ) the direct isentrope calculations presented in Knudson et al. 7 obtained through thermodynamic integration using vDW-DF2, we wish to comment now on the apparent discrepancy between these results and the estimated temperature reduction factor T 2 /T 1 = 0.66 presented by Celliers et al. 12 .Their estimate was obtained from analysis of their Eqs.( 1) and (2).However, in computing their Eq.( 2), γ was incorrectly treated as a constant over the entire 95 GPa range of the integral.As is clear from the temperature minimum of the isentropes obtained by thermodynamic integration in Fig. 1 of Knudson, et al. 7 and suggested also by Fig. S32 of the Supplementary Materials for Celliers et al. 8 along with Fig. 1 of Celliers et al. 12 , the Grüneisen γ → 0 and eventually turns positive over the course of several 10s of GPa beyond the phase boundary.This follows directly from (∂S/∂P ) T = −γC V /B T and is a direct consequence of exhausting the molecular to atomic transition.That γ → 0 within approximately 50 GPa is illustrated in Fig. 4 where (∂P/∂T ) V is essentially zero by 255 GPa for temperatures between 1200 K and 1400 K. Correcting the integration in Eq. ( 2) of Ref. 12, and accounting for the consequences in subsequent steps in that analysis, brings those estimates in line with the results obtained here.
A refinement of the finite-difference analysis outlined in Ref. 8 is presented in Ref. 12.This analysis combines the Clausius-Clapeyron relation with the assumption that the NIF and Z pressures are indicative of the entrance (P I 1 , T 1 ) and exit (P M 2 , T 2 ) of the coexistence region.By effectively constraining the 95 GPa pressure difference to the coexistence line, they arrive at a temperature reduction factor due to latent heat of T 2 /T 1 = 0.58.This value of 0.58, when combined with Eq. ( 9), and T 1 = 1447 K, implies an anomalously low C coex = 0.74 k B /atom, in gross disagreement with the C coex = 2.64 k B /atom obtained from our first-principles calculations.We emphasize that this larger value for C coex is consistent across vdW-DF1, vdW-DF2, and CEIMC; even though these different first-principles frameworks disagree on the precise location of the insulator-metal phase transition boundary, they are in quite good agreement on the thermodynamics of the transition (i.e. the specific heats, Grüneisen gamma, bulk modulus, etc.).
As noted in the Introduction, the criteria for identifying the location of the phase boundary are substantively different between the Z and NIF experiments.In the Z experiments the phase boundary was associated with the abrupt drop from a high reflectivity phase upon pressure release.The NIF experiments only probed the transition upon compression; there the phase boundary was associated with the reflectivity rising above a threshold of 30%.Both approaches are valid; however, comparing one to the other on unequal footing exacerbates the apparent discrepancy.Associating the phase boundary in the NIF experiments with a higher value of reflectivity, commensurate with completion of the transition, would result in a higher inferred value for the transition pressure.As a quantitative illustration, consider the calculations of Rillo et al. 27 , which suggest a reflectivity in excess of 40% at the completion of the transition at a temperature of 1500 K. Applying this criteria to the 1450 K (closest available and specifically N150914-2) reflectivity versus pressure trace in Fig. 2 of Celliers et al. 8 suggests a phase boundary exit pressure of 240 GPa.
Evaluating the completion of the transition would also result in a lower inferred temperature, in accordance with the latent heat considerations explored here.However, regarding the estimated temperature, as noted in the Introduction two of the three equation of state (EOS) models used to infer the temperature in the NIF experiments include PBE latent heat contributions at pressures below the observed transition pressure (in either the Z or NIF interpretations) and therefore result in lower predicted temperatures.The 2003 deuterium EOS of Kerley 28 , the same EOS used for estimating temperatures for the Z experiments prior to switching to DFT calculations of the isentrope, is much closer to, and slightly exceeds, the upper bound of the temperature estimates provided in Fig. 3 of Celliers et al. 8 .Calculation of the temperature path for this reported 1450 K transition case with the Kerley 2003 EOS is shown in Fig. S18 of the Supplemental Material for Celliers et al. 8 ; subtracting 200 K for latent heat at 240 GPa, as indicated by the calculations in this paper, suggests a temperature closer to 1625 K.We note that a transition at 240 GPa and 1625 K is much closer to the vdW-DF2 phase boundary than that of vdW-DF1.Interpreting this NIF experiment on the same footing as the Z experiments reduces the apparent discrepancy between the two experiments by about a factor of two, to approximately 42 GPa.The remaining differences are not insignificant, but will require future experiments and analyses to reconcile.

VI. CONCLUSION
We have performed an extensive study of the thermodynamics of the insulator-metal transition in dense liquid deuterium within a first-principles framework to assist in interpreting recent dynamic compression experiments.Specifically, we used density functional theory, including nuclear quantum effects, to directly calculate the temperature drop for an isentrope that traverses the first-order insulator-metal transition.This was accomplished by evaluating an exact expression for change in entropy along a line in (P, T ) space in terms of quantities readily calculated within an N V T first-principles framework.An extensive set of path-integral molecular dynamics calculations with the vdW-DF1 functional were performed to obtain quantitative values for C V , γ, and B T adjacent to the metallic side of the phase transition in the region traversed in recent dynamic compression experiments.The resulting temperature drops were found to be consistent with previous direct isentrope calculations 7 obtained through thermodynamic integration and estimates 11 based on interpolations of CEIMC data 10 .Furthermore, these temperature drops are in stark disagreement with a recent reinterpretation presented by Celliers et al. 8,12 The arguments presented in Celliers et al. 8,12 have their root in constraining both the Z and NIF experiments to a given theoretical phase boundary and approximate isentrope, with the NIF experiments marking the entrance to the coexistence region and the Z experiments probing the exit.For the quoted pressures in the two sets of experiments, forcing the Z experiments to the phase boundary would require a factor of two correction downward in the temperature estimates for the Z experiments.However, as is clear from the analysis presented here, any given pressure difference between the two experiments would, by this enforced temperature constraint, result in a different effective specific heat along the phase boundary, generally in conflict with the underlying thermodynamics.The large temperature corrections required for the estimated 95 GPa pressure difference would necessitate an anomalously low specific heat that is in gross disagreement with first-principles calculations presented here using vdW-DF1, those in Pierleoni et al. 10 using CEIMC as argued in Desjarlais et al. 11 , and our earlier direct calculations of the isentropes 7 obtained through thermodynamic integration using vdW-DF2.While interpreting the 1450 K NIF experiment in a manner analogous to that used for the Z experiments suggests a substantially smaller discrepancy between the two sets of experiments, an approximately 42 GPa difference remains.What is clear from the analysis presented here is that the pressure difference cannot be explained through the supposition of a very large latent heat effect between the onset and completion of the phase transition.

ACKNOWLEDGMENTS
MPD would like to thank Dario Alfè for providing many of the path-integral algorithms.Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA0003525.This paper describes objective technical results and analysis.Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.RR also thanks the Deutsche Forschungsgemeinschaft (DFG) for support via the SFB 652 and FOR 2440.

SS + DS 1 SFIG. 1 .
FIG.1.Notional schematic of the P T phase diagram for deuterium around the insulator-metal transition with representative isentropes.The black dashed line represents the phase boundary and terminates at the critical point.The isentropes are depicted with slight −dT /dP near the phase boundary, resulting from a negative Grüneisen γ, in accordance with firstprinciples calculations.The red arrow suggests a convenient integration path from S + ∆S1 back to S.

FIG. 3 .
FIG.3.Pressure × volume versus energy in the vicinity of the insulator-metal transition at specific volumes of 1.80 (green circles, dotted fit) and 1.78 (black circles, solid fit) Å3 /atom.The Grüneisen γ is given by the slope.

FIG. 4 .FIG. 5 .
FIG. 4. Pressure versus temperature for several specific volumes.Labels for each fit indicate the specific volume in Å3 /atom.