Elucidation of the subcritical character of the liquid--liquid transition in dense hydrogen
aa r X i v : . [ c ond - m a t . o t h e r] D ec Elucidation of the subcritical character of the liquid–liquid transition in densehydrogen
Valentin V. Karasiev, ∗ Joshua Hinz, S. X. Hu, and S.B. Trickey Laboratory for Laser Energetics, University of Rochester,250 East River Road, Rochester, New York 14623 USA Quantum Theory Project, Dept. of Physics, University of Florida, Gainesville FL 32603 USA (Dated: Dec. 26, 2020)
ARISING FROM Bingqing Cheng et al.
Na-ture https://doi.org/10.1038/s41586-020-2677-y(2020)
Determining the liquid–liquid phase transition (LLPT)in high-pressure hydrogen is a longstanding challengewith notable variation in experimental and calculated re-sults. See Refs. 1–4 and works cited therein for bothcalculational and experimental developments. Until re-cently, the computational consensus was for a first-ordertransition. Calculated values differed but, for example,our results on 700 ≤ T ≤ ≥ P ≥
70 GPa [2]. Driven by molecular H dissoci-ation, transition signatures include density jumps, qual-itative and sharp changes in ionic pair correlation func-tions (PCFs), and abrupt dc conductivity and reflectivitychanges. Coupled-electron ion Monte Carlo (CEIMC) [5]results concur at least roughly with those from ab-initio molecular dynamics driven by consistent density func-tional theory (MD-DFT) [2] and show reasonable agree-ment with experiment also.In marked contrast, Cheng et al. [6] found a continu-ous transformation from a molecular to an atomic liquidthat goes supercritical above P ≈
350 GPa, T ≈ N V T dynamics tending to increased defectconcentration compared to that from the
N P T ensem-ble. The other is much shorter simulation times in theMD-DFT and CEIMC calculations than in the MD-MLPones.Those diagnoses implicate other issues. Almost allof the MLP training was on small systems ( N ≤ N P T
MD simulations driven byDFT forces with PBE exchange-correlation [7]. (PBE was used in Ref. [6] to train the MLP.) We used systemsizes from 256 through 2048 atoms per cell. Brillouinzone sampling used the Baldereschi mean value point forthe simple cubic crystal structure k = ( , , ) [8]. Vasp [9, 10] was used for 1024 and 2048 atom systems, whilethe i-PI interface [11] with
Quantum Espresso [12] wasused for 256 and 512 atoms. Consistent results from thetwo shows that the MD code and technical choices (ther-mostat, barostat, etc.) are inconsequential.Our new large-system MD-DFT results agree with pre-vious DFT-based and CEIMC simulations [2, 3, 13]:there is a sharp molecular-to-atomic transition. Thequalitatively different character compared with whatcomes from the MD-MLP is shown in Fig. 1. The left-column panels show density profiles ρ H ( T ) along isobars.At 350 and 300 GPa, the large-scale MD-DFT ρ H ( T )values jump ≈
1% near T = 650 K. At 300 GPa, this isabove the melting temperature T m [14]. In contrast, the300-GPa MD-MLP isobar has a steep density increasenear T = 500 K (in the stable solid phase) [6], but passessmoothly through both that melt line and the LLPT. Ex-cept for a systematic offset, the MD-MLP ρ H ( T ) matchesthe MD-DFT ρ H ( T ) in the atomic fluid region.Figure 1 also shows unequivocally that there are noimportant finite-size effects on the calculated LLPT. Thedensity profiles on each of the isobars ( P = 250, 200, 150,and 100 GPa) are almost identical irrespective of atomcount (256, 512, 1024, or 2048). The transition characteris insensitive to system size and specific technical choicesof the MD code used, while the transition temperature T LLPT is affected only modestly. At P = 200 GPa, forexample, going from 256 to 2048 atoms decreases T LLPT by less than 100 K; ρ H values jump ≈
3% in MD-DFTsimulations for all system sizes. A 512 atom system seemsadequate to eliminate any major finite-size effects. Thisoutcome agrees with Ref. [15]. Those authors foundthat four well-defined molecular shells in the PCF of a3456-atom system were captured quite well in a 500-atomsupercell calculation.The molar heat capacity from MD-DFT as a functionof T is shown in Fig. 1, middle column. All the isobarsexhibit divergent heat capacity character across the tran-sition. They confirm that finite-size effects on T LLPT aresmall and do not modify that character.
FIG. 1. Comparison of MD results from the PBE exchange-correlation-based machine-learning potential (MLP) and ab initio
MD-DFT (DFT)
NP T simulations. Left column panels (a): Hydrogen density as function of T along six isobars. Meltingtemperature T m for each isobar is shown by a vertical dashed line [14]. Middle column panels (b): Molar heat capacity as afunction of T along the isobars. Right column panels (c): Pair correlation function (PCF) for each isobar for two temperaturesbelow the density jump and two above. Figure 1 right-hand column shows the PCF on each iso-bar at pairs of temperatures below and above the densityjump. Above, the first PCF peak virtually disappears,confirmation of the density jump being in conjunctionwith the molecular dissociation [2].To test possible long simulation duration effects on T LLPT or its character, we did up to six sequential MD-DFT runs of roughly 1.8-ps duration each for a total of ≈ N P T simulations beginningwith the atomic fluid at 200 GPa. Starting at 950 K,we cooled the system in sequential runs to 899, 849, and824 K with simulation durations around 8 ps for eachtemperature. If the nanosecond timescale were to yielda smooth transition, the hydrogen density during such afast cooling curve would not drop sharply below the hy-pothetical smooth long-duration curve. But, as evident
FIG. 2. The LLPT boundary from the present large-scale MD-DFT (DFT/PBE) simulations compared to MLP(MLP/PBE) C max P and ρ max curves. in the Fig. 1 density plot at 200 GPa (left column), thecooling curve (thin blue curve, circles), is almost identi-cal to the one from MD simulations when the molecularfluid T is increased gradually (sharp transition shown bythe solid orange curve).Figure 2 shows the LLPT curves associated with den-sity jumps, heat capacity peaks, and PCF peak disap-pearance. For the new large-scale MD-DFT calculations,those three criteria give one curve (virtually identical P, T values), shown in red with squares at data points.Two MD-MLP curves emerge from the analysis, how-ever, one for the location of molar heat capacity maxima C max P , and another for the maximum density, ρ max . Con-sistent with the foregoing discussion, there are strikingdifferences. The MLP C max P curve lies well below theMD-DFT curve. The MLP ρ max curve is flatter thanthe MD-DFT reference curve and lies close to it only atabout P = 70 GPa, T = 2800 K and then again for P between about 170 and 300 GPa.Given that neither the finite-size effect nor simulationduration diagnosis advanced by Cheng et al. [6] is sus-tained by our direct exploration, the remaining plausiblecause of the different physics they found must be in theMLP. The detailed origin of that different physics is a bitobscure. However, as discussed in our Supplemental In-formation, documentation in the Supplemental Informa-tion to Ref. 6 confirms that the MLP does not reproducethe behavior (be it physical or not) of several MD-DFTcalculations. Those differences, in addition to the starkLLPT differences discussed here, confirm that the MLP isnot systematically related to the physics of a well-definedBorn–Oppenheimer electronic structure treatment of theH system. The MD-MLP results instead are consistent,at least, with the MLP being a single interpolative, ap- proximate representation of the electronic structure oftwo chemically distinct regimes (molecular, atomic) ofthe hydrogen liquid.We conclude that the numerical evidence for super-critical behavior of high-pressure liquid hydrogen basedon the approximate MLP simulations is unsupported byMD-DFT simulations on much larger systems for signif-icantly longer durations. The diagnosis of the differencebetween MD-MLP and MD-DFT calculations as beingfrom size and duration effects is mistaken. Rather, thesupercritical behavior found in the MD-MLP calculationsseems plausibly to be an artifact of a disconnect of theMLP from underlying electronic structure differences in-herent in the chemistry of the LLPT. Data availability
The data that support the findings shown in the figuresare available from the corresponding author uponreasonable request. ∗ [email protected][1] Gregoryanz, E. et al. Everything you always wanted toknow about metallic hydrogen but were afraid to ask. Matter Radiat. Extremes , 038101 (2020).[2] Hinz, J. et al. Fully consistent density functional theorydetermination of the insulator-metal transition bound-ary in warm dense hydrogen. Phys. Rev. Research ,032065(R) (2020).[3] Rillo, G., Morales, M. A., Ceperley, D. M. & Pierleoni, C.Optical properties of high-pressure fluid hydrogen acrossmolecular dissociation. Proc. Natl. Acad. Sci. , 9770-9774 (2019).[4] Lu, B., Kang, D., Wang, D., Gao, T. & Dai, J. To-wards the same line of liquid–liquid phase transitionof dense hydrogen from various theoretical predictions.
Chin. Phys. Lett. , 103102 (2019).[5] Pierleoni, C., Morales, M. A., Rillo, G., Holzmann, M.& Ceperley, D. M. Liquid–liquid phase transition in hy-drogen by coupled electron–ion Monte Carlo simulations. Proc. Natl. Acad. Sci. , 4953-4957, (2016).[6] Cheng, B., Mazzola, G., Pickard, C. J. & Ceriotti, M. Ev-idence for supercritical behaviour of high-pressure liquidhydrogen.
Nature , 217-220 (2020).[7] Perdew, J. P., Burke, K. & Ernzerhof, M. Generalizedgradient approximation made simple.
Phys. Rev. Lett. , 3865-3868 (1996); , 1396(E) (1997).[8] Baldereschi, A. Mean-value point in the Brillouin zone. Phys. Rev. B , 5212-5215 (1973).[9] Kresse, G. & Furthm¨uller, J. Efficient iterative schemesfor ab initio total-energy calculations using a plane-wavebasis set. Phys. Rev. B , 11169-11186 (1996).[10] Kresse, G. & Joubert, D. From ultrasoft pseudopotentialsto the projector augmented-wave method. Phys. Rev. B , 1758-1775 (1999).[11] Kapil, V. et al. I-PI 2.0: A universal force engine for ad- vanced molecular simulations. Comput. Phys. Commun. , 214-223 (2019).[12] Giannozzi, P. et al. Advanced capabilities for materialsmodelling with quantum espresso.
J. Phys.: Condens.Matter , 465901 (2017).[13] Lorenzen, W., Holst, B. & Redmer, R. First-order liquid-liquid phase transition in dense hydrogen. Phys. Rev. B , 195107 (2010).[14] Zha, C. S., Liu, H., Tse, J. S. & Hemley, R. J. Meltingand high P–T transitions of hydrogen up to 300 GPa. Phys. Rev.Lett. , 075302 (2017).[15] Geng, H. Y., Wu, Q., Marqu´es, M. & Ackland, G.J. Thermodynamic anomalies and three distinct liquid-liquid transitions in warm dense liquid hydrogen.
Phys.Rev. B , 134109 (2019).
Acknowledgements
This report was prepared as an account of work spon-sored by an agency of the U.S. Government. Neitherthe U.S. Government nor any agency thereof, nor any oftheir employees, makes any warranty, express or implied,or assumes any legal liability or responsibility for the ac-curacy, completeness, or usefulness of any information,apparatus, product, or process disclosed, or representsthat its use would not infringe privately owned rights.Reference herein to any specific commercial product, pro-cess, or service by trade name, trademark, manufacturer,or otherwise does not necessarily constitute or imply itsendorsement, recommendation, or favoring by the U.S.Government or any agency thereof. The views and opin-ions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agencythereof.V.V.K., J.H., and S.X.H. were supported by theDepartment of Energy National Nuclear Security Ad-ministration Award Number DE-NA0003856 and USNational Science Foundation PHY Grant No. 1802964.S.B.T. was supported by Department of Energy GrantDE-SC 0002139. This research used resources of theNational Energy Research Scientific Computing Center,a DOE Office of Science User Facility supported bythe Office of Science of the U.S. Department of Energyunder Contract No. DE-AC02-05CH11231. Part of thecomputations were performed on the Laboratory forLaser Energetics HPC systems.
Author contributions
V.V.K. conceived the projectinitially and designed the study. V.V.K. and J.H. per-formed the MD-DFT simulations and post-processed thedata. V.V.K. wrote the initial manuscript. S.B.T. re-vised the conception and scope. V.V.K. and S.B.Trewrote the manuscript. All authors discussed the re-sults and revised the paper extensively.
Conflict of interests
The authors declare that theyhave no conflicts of interest.
Additional informationSupplementary information is available for this paperat https://doi.org/10.1038/xx