Hydrogen liquid-liquid transition from first principles and machine learning

Abstract

The molecular-to-atomic liquid-liquid transition (LLT) in high-pressure hydrogen is a fundamental topic touching domains from planetary science to materials modeling. Yet, the nature of the LLT is still under debate. To resolve it, numerical simulations must cover length and time scales spanning several orders of magnitude. We overcome these size and time limitations by constructing a fast and accurate machine-learning interatomic potential (MLIP) built on the MACE neural network architecture. The MLIP is trained on Perdew-Burke-Ernzerhof (PBE) density functional calculations and uses a modified loss function correcting for an energy bias in the molecular phase. Classical and path-integral molecular dynamics driven by this MLIP show that the LLT is always supercritical above the melting temperature. The position of the corresponding Widom line agrees with previous ab initio PBE calculations, which in contrast predicted a first-order LLT. According to our calculations, the crossover line becomes a first-order transition only inside the molecular crystal region. These results call for a reconsideration of the LLT picture previously drawn.

Figures (4)

• Figure 1: Results for MD simulations obtained with the MACE model for different temperatures: (a) $T=900$ K, (b) $T=950$ K, (c) $T=1000$ K, (d) $T=1100$ K. For each temperature, the left panel reports the EOS close to the LLT for the different system sizes. For $T=900$ K the red shaded area highlights the $r_s$ range for which hysteresis is observed at all values of $N$. The right panels show the average $\max_{\mathbf{k}} S(\mathbf{k}) /N$ along the trajectory (top right) and the average molecular fraction $m$ (bottom right) as a function of the system size $N$ for two different values of the Wigner-Seitz radius $r_s$.
• Figure 2: (a) Value of $\max_{\mathbf{k}}S(\mathbf{k}) / N$ as a function of the simulation time for two systems, with $N=256$ and $N=512$ atoms, respectively, at temperature $T= 950$ K and $r_s = 1.43$. Colored lines represent running averages over a time window of $8$ ps. (b) Results for $200$ ps long simulations, $N=128$, $T=1400$ K and corresponding confidence intervals, estimated from a running average with time window $\tau_{\textrm{run}} = 4$ ps. The AIMD result reported in Ref. Tirelli2022 is also shown. (c) Same as (b) but with $N=512$, $T=1000$ K and $\tau_{\textrm{run}} = 10$ ps.
• Figure 3: Specific heat per particle vs pressure along the isotherms, for $N = 2048$ obtained with the MACE model. The shaded areas indicate the uncertainty in the peak position. In the inset, we report the size scaling of the maximum of $c_v/k_B$ for four different temperatures (in log-log scale). The results confirm a first-order transition for $T=900$ K, where $c_v$ scales with a behavior close to linear, and a smooth crossover for higher temperatures in the thermodynamic limit.
• Figure 4: Classical PBE-LLT location as computed with different methods. AIMD results by Ref. Lorenzen2010 (light blue dashed line), Ref. Morales2010 (green dashed line), Ref. Karasiev2021 (blue dashed line). The results for the molecular-to-atomic crossover obtained with an NN MLIP by Ref. Cheng2020 are reported with orange and violet markers, corresponding to the maximum of the isobaric specific heat and density, respectively. The black markers and line indicate the recently proposed PBE melting line by Ref. Niu2023, obtained with an NN MLIP and the two-phase method. The MACE model results are indicated with red markers. The filled point at $T= 900$ K indicates the first-order character of the transition, while the empty points correspond to the location of the Widom line given by the $c_v$ maximum (see Fig. \ref{['fig: CV nvt']}).