Ammonia-Hydrogen Phase Equilibria in Planetary Envelopes from First-Principles-Based Machine-Learning Molecular Dynamics
Authors: Akash Gupta, Austin Guo, Jie Deng
Organizations: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA · Department of Geosciences, Princeton University, Princeton, NJ 08544, USA
Ammonia and hydrogen are major constituents of planets within and beyond the solar system, yet how they interact in planetary envelopes is poorly known. We perform large-scale two-phase coexistence simulations of the ammonia-hydrogen system from 500 to 4000 K and 0.1 to 100 GPa, driven by machine-learning interatomic potentials trained on density functional theory calculations. From these we map the compositions of coexisting hydrogen-rich and ammonia-rich phases, and the critical curve that separates conditions of complete miscibility from those of phase separation. The critical temperature increases with pressure, from 500 K at 0.1 GPa to 2000 K near 30 GPa, and we find no phase separation at 3000 and 4000 K up to 100 GPa. The slope of the critical curve decreases substantially between 0.3 and 30 GPa, and the curve nearly flattens above about 30 GPa and 2000 K, where molecules begin to give way to atomic species. The coexistence curves are asymmetric with an ammonia-rich critical composition. For the expected envelope compositions, the phase boundary lies close to, and intersects, published thermal profiles of Uranus and Neptune, implying ammonia may be raining out in these planets today. In contrast, in warm sub-Neptunes such as TOI-270 d, ammonia and hydrogen are likely miscible in all proportions. Together with previous studies, our results suggest that silicates, water and ammonia separate from hydrogen in sequence as planets cool, with implications for their compositional, thermal and structural evolution.
Figures & tables
Figure 1: A two-phase coexistence simulation at 1000 K and 7 GPa, containing 4096 H 2 and 2048 NH 3 molecules (16,384 atoms). Blue spheres are nitrogen atoms and gray spheres are hydrogen atoms. Top: the initial condition, in which slabs of pure ammonia and pure hydrogen are joined (left), followed by a perspective view and a side view of the system after 2 ns. The brown surfaces are the Gibbs dividing surfaces between the ammonia-rich and hydrogen-rich phases [ 123 ] . The rightmost panel shows the mass density along the long axis of the cell after 2 ns; symbols are simulation results and the solid curve is the theoretically motivated hyperbolic-tangent profile expected for coexisting phases [ 17 , 120 ] . Dashed lines mark the determined interface positions and shaded bands reflect the interface widths. Bottom: total energy and composition of each phase as a function of time; the final 0.5 ns (darker) is the window over which the compositions are averaged.
Figure 2: Compositions of coexisting phases as a function of pressure at six temperatures. Filled circles are the measured compositions of each phase in simulations that result in two phases, with bands indicating their 2 σ variation in time, and open circles show the compositions of simulations that remain a single homogeneous fluid. The solid curves are the posterior-median coexistence curves from our thermodynamic model, with shaded 2 σ envelopes (see Methods ). The simulations have bulk compositions of xH2= 0.67 and 0.33, and contain between 16,384 and 70,304 atoms. No simulation at 3000 or 4000 K separates into two phases at any pressure that we sampled and thus the coexistence curve shown at 3000 K is only partially constrained.
Figure 3: Phase diagram of the ammonia (NH 3 )-hydrogen (H 2 ) system in pressure and temperature space. The bold black line with its gray 2 σ band is the critical curve from our thermodynamic model delineating the regimes of phase stability (turquoise shading). Above 2000 K the critical curve is only partially constrained because none of the simulations at 3000-4000 K phase-separate up to 100 GPa. The bold dashed line below 500 K is physically motivated by previous work on hydrogen-water fluids [ 97 , 34 ] and joins our results to the critical point of pure ammonia (pink square) and its condensation curve [ 29 ] . The gray solid line, with its 2 σ band, denotes the isopleth of the coexistence surface at xH2∼0.9 , as expected in the deep envelopes of Uranus and Neptune [ 72 , 93 ] . It is extended to experimental results from Reamer and Sage [88] and Heise [37] . Also shown are model P-T profiles for Uranus (orange) and Neptune (purple) from Nettelmann et al. [72] (solid lines) and Tejada Arevalo [111] (dashed lines). The red lines are P-T profiles of TOI-270 d from Rigby and Madhusudhan [89] for internal temperatures of 25 K (dark red) and 50 K (light red). The pressure axis is logarithmic at low pressure and linear at high pressure. Comparison of the critical curve and this isopleth with the model profiles of Uranus, Neptune, and TOI-270 d suggests that while ammonia and hydrogen are likely fully miscible in the envelopes of warm sub-Neptunes, Uranus and Neptune may be undergoing ammonia rain today, at depths that depend on assumptions about the initial conditions and evolution of these planets.
Figure 4: Variation with pressure of the slope of the critical curve (top) and of the speciation of the homogeneous fluid (bottom). In the top panel, the line shows the slope dTc/dPc of the critical curve, and the shaded band its 2 σ uncertainty. The bottom panel gives the fraction of each species in a homogeneous fluid with xH2=0.67 along the critical curve up to 3000 K and then at 4000 K at the highest pressures. The gray dashed curve is a Fermi-Dirac-like function tracing the change in atomic H fractions, showcasing the molecular-to-atomic transition of the fluid at about 30-100 GPa and 2000-4000 K. Species are identified from atomic neighbors that persist for 100 fs. The pressure axis is logarithmic below 15 GPa and linear above.
Appendix figures & tables4 assets
Supplementary material from the paper’s appendix.
Appendix
Figure S1: Training-set performance of the final machine-learning potential (MLP) across the sampled pressure and temperature space. The panels show the mean absolute energy error per atom and the root-mean-square errors of the force and stress components for 68,994 training configurations, containing 16,128,230 atoms in total. Configurations are binned by the midpoint temperature of their training zone (Methods) and by pressure, which includes the ionic kinetic contribution at that temperature. Values in parentheses mark bins with fewer than 300 configurations, and white cells contain none. The energy and force color bars are linear, and the stress color bar is logarithmic.
Figure S2: Comparison of the machine-learning potential (MLP) predictions with density functional theory (DFT) calculations, and the dependence of our results on system size. The top panels compare the energies, forces and static stresses predicted by the MLP with those from DFT for 61,562 atoms, corresponding to 258 energies, 184,686 force components and 1,548 stress components at 500-4000 K; the dashed lines indicate perfect agreement. The root-mean-square errors are 19 meV atom -1 in energy, 0.39 eV Å -1 in force and 0.96 GPa in stress. The bottom left panel compares the H-H, H-N and N-N radial distribution functions from the MLP (solid lines) and from DFT (open circles) at 1000 K and 30 GPa. The bottom right panel shows the compositions of the hydrogen-rich phase (squares) and the ammonia-rich phase (circles) at 1000 K and 10 GPa as a function of the number of atoms in the simulation; the shaded bands show the frame-to-frame fluctuations, and the filled symbols mark the chosen 16,384-atom system.
Figure S3: Temperature dependence of the mixing parameters of the thermodynamic model (Equation A2 ): the dimensionless asymmetry parameter r (left), the zero-pressure excess Gibbs energy scale u in kJ mol -1 (middle), and the excess volume scale v in kJ mol -1 GPa -1 (right), on a logarithmic temperature axis. Solid curves show pointwise posterior medians, and shaded bands 2 σ intervals, conditional on the model. The selected model has a temperature-independent r , while u is quadratic in t=T/(1000K) and v contains a term proportional to t−4 .
Figure S4: Coexistence (binodal) curves predicted by the thermodynamic model at fixed pressures of 0.3-100 GPa, shown as temperature against composition on a logarithmic temperature axis. Below each curve a mixture of that composition separates into two phases. Solid curves show posterior medians, and shaded bands 2 σ intervals. Portions below 500 K and above 2000 K extend beyond the temperatures with observed two-phase states.
Simulating microstructure evolution requires quantum-mechanical accuracy and mesoscopic reach in length and time scales, a combination that no current method achieves. Classical phase-field models provide this reach, but their accuracy is limited by phenomenological free energies and mobilities. Here we develop a framework for learning ab initio phase-field models, where the mesoscopic equation is not postulated but derived from a Mori-Zwanzig projection of molecular dynamics onto species-density fields under explicit assumptions. The nonlocal free energy and mobility left unspecified by this equation are parametrized by neural networks and learned from short molecular dynamics trajectories generated with machine-learning interatomic potentials of ab initio accuracy. We demonstrate the framework on an iron-boron melt and on hydrogen-helium mixtures under planetary conditions. For iron-boron, the model shows that the melt at the FeB4 composition is spinodally unstable at ambient pressure but stabilized at 10 GPa, offering a thermodynamic rationale for why FeB4 has been synthesized only under high pressure. For hydrogen-helium, the model predicts the immiscibility boundary and captures droplet nucleation and growth in helium-rain simulations of a column corresponding to 2.2 million atoms, far beyond the scale of atomistic modeling at comparable accuracy. Trained across compositions and conditions, such models could provide a mesoscopic counterpart to ab initio molecular dynamics.
Mengyi Chen, Peichen Zhong, Zihan Zhang +1
Department of Mathematics, National University of Singapore, Singapore · Department of Materials Science and Engineering, National University of Singapore, Singapore · Institute for Functional Intelligent Materials, National University of Singapore, Singapore
Accurate prediction of phase equilibria remains a central challenge in chemical engineering. Physics-consistent machine learning methods that incorporate thermodynamic structure into neural networks have recently shown strong performance for activity-coefficient modeling. However, extending such approaches to equilibrium data arising from an extremum principle, such as liquid-liquid equilibria, remains difficult. Here we present DISCOMAX, a differentiable algorithm for phase-equilibrium calculation that guarantees thermodynamic consistency at both training and inference, only subject to a user-specified discretization. The method combines discrete enumeration of feasible phase states with masked softmax aggregation in the backward pass, with the propagation of the true equilibrium state in the forward pass, using a straight-through gradient estimator to enable physics-consistent end-to-end learning of neural \gls{gE}-models. We show that this approach bears analogy to statistical thermodynamics, and we evaluate it on binary liquid-liquid equilibrium data where it outperforms existing surrogate-based methods, while offering a general framework for learning from different kinds of equilibrium data.
Karim K. Ben Hicham, Moreno Ascani, Jan G. Rittig +1
RWTH Aachen University, Process Systems Engineering (AVT.SVT), 52074 Aachen, Germany · Forschungszentrum Jülich GmbH, Institute of Climate and Energy Systems ICE-1: Energy Systems Engineering, 52425 Jülich, Germany · JARA-ENERGY, 52056 Aachen, Germany
The crystal structure of high-pressure solid hydrogen remains a fundamental open problem. Although the research frontier has mostly shifted toward ultra-high pressure phases above 400 GPa, we show that even the broken symmetry phase observed around 130~GPa requires revisiting due to its intricate coupling of electronic and nuclear degrees of freedom. Here, we develop a first principle quantum Monte Carlo framework based on a deep neural network wave function that treats both electrons and nuclei quantum mechanically within the constant pressure ensemble. Our calculations reveal an unreported ground-state structure candidate for the broken symmetry phase with Cmcm space group symmetry, and we test its stability up to 96 atoms. The predicted structure quantitatively matches the experimental equation of state and gives the closest x-ray diffraction peak-position match among the tested candidates. Furthermore, our group-theoretical analysis provides a symmetry-counting compatibility check between the Cmcm structure and existing Raman and infrared spectroscopic data. Crucially, static density functional theory calculation reveals the Cmcm structure as a dynamically unstable saddle point on the Born-Oppenheimer potential energy surface, demonstrating that a full quantum many-body treatment of the problem is necessary. These results shed new light on the phase diagram of high-pressure hydrogen and call for further experimental verifications.
Shengdu Chai, Chen Lin, Xinyang Dong +4
Interdisciplinary Center for Theoretical Physics and Information Sciences (ICTPIS), Fudan University, Shanghai 200433, China · Shanghai Artificial Intelligence Laboratory, Shanghai 200232, China · Department of Engineering, University of Oxford, Oxford OX1 4BH, UK +4