Both sub-Neptunes and super-Earths may have begun with a hydrogen-rich envelope gravitationally captured from the surrounding stellar disk during accretion. Sub-Neptunes retained some portion of this primary atmosphere to the present day, whereas in the case of super-Earths, the envelope was subsequently lost7. The population of super-Earths and sub-Neptunes, and the origin of the radius valley that separates these two classes of planets, is best explained by cores that are made of an Earth-like composition without a substantial amount of accreted ice8,9,10,11. For sub-Neptunes, the hydrogen-rich envelope overlies the rocky core for billions of years, whereas for super-Earths, the envelope may be retained for about 100 Myr (refs. 2,8,12). During this time, the hydrogen envelope and the rocky core may have reacted chemically with each other, altering the composition of the atmosphere and interior. Our own Earth may have also once had a primary atmosphere that left its mark on the present-day composition of its interior13,14.
The structure and evolution of sub-Neptunes and super-Earths challenge our understanding of fundamental properties of materials at extreme conditions. For example6, at the core–envelope boundary of a typical sub-Neptune, the temperature may be 5,000 K, whereas the pressure may be 5 GPa. Both core and envelope are fluid at the conditions of the interface, which exceed the silicate melting temperature: the rocky core forms a magma ocean6,15,16. This part of the pressure–temperature space is not much explored by experiments or theory. Similar temperatures are found in the interiors of giant planets, but only at much higher pressures, at which dynamic compression experiments are able to probe17. By the same token, similar pressures are found in the interiors of rocky planets, but only at much lower temperature conditions that are readily accessible to static compression experiments18. Previous modelling studies have emphasized the importance of a better understanding of fluid silicate–hydrogen phase equilibria for super-Earth and sub-Neptune structure and evolution19. Experiments have demonstrated that small concentrations of hydrogen are soluble in silicate liquid, albeit at pressure–temperature conditions much less extreme than those we consider here18.
We use the phase coexistence method20,21 to study the reaction between hydrogen fluid and a liquid of MgSiO 3 composition, taken to be representative of the rocky portion of planets (Fig. 1 and Methods). The initial condition consists of equal-sized domains of pure hydrogen and pure silicate composition, pre-equilibrated to the same pressure and temperature. The Hellmann–Feynman forces are computed in density functional theory. As the molecular dynamics simulation proceeds, a reaction occurs, and a dynamic equilibrium is established. At temperatures below the critical temperature, two phases with distinct compositions coexist in equilibrium. By quantifying the compositions of the two coexisting phases, we map out the phase diagram and identify the critical temperature above which a single homogeneous fluid is stable. To gain additional insight into our results, we also perform a series of simulations of homogeneous fluids with compositions on the MgSiO 3 –H 2 join, from which we extract the energetics of solution. Further details are given in the Methods.
Fig. 1: Illustration of the phase coexistence method. a,b, Snapshots from a phase coexistence simulation at 10 GPa and 2,500 K: the initial condition (a) and from the equilibrated portion of the trajectory (b). H, white; Mg, yellow; Si, blue; O, red; and Gibbs dividing surfaces, light blue. c, The one-dimensional density profile from the equilibrated portion of the trajectory (symbols with colours corresponding to atom types); lines are fits to equation (4). Full size image
Over a range of conditions, we find two phases coexisting in dynamic equilibrium, one hydrogen-rich and the other hydrogen-poor (Fig. 1). In both phases, the silicate fraction nearly maintains MgSiO 3 stoichiometry, that is, MgSiO 3 behaves as a component, and the system is nearly binary. We choose as our compositional variable the molar fraction of the H 2 component
$$x=\frac{{N}_{{\text{H}}_{2}}}{{N}_{{\text{H}}_{2}}+{N}_{\mathrm{MgO}}+{N}_{{\text{SiO}}_{2}}}=\frac{{N}_{\text{H}}/2}{{N}_{\text{H}}/2+2({N}_{\text{Mg}}+{N}_{\text{Si}}+{N}_{\text{O}})/5}$$ (1)
where N i is the number of molecules or atoms of type i. We quantify the composition of the two coexisting phases based on (1) one-dimensional density profiles using a Widom-like expression22 and (2) a three-dimensional coarse-grained density field23 (see Methods for further details).
The compositions of the two coexisting phases approach each other with increasing temperature (Fig. 2a and Extended Data Table 1) and merge into a single homogeneous phase at temperatures greater than the critical temperature T C . We find that the compositions of the two phases are asymmetric: the hydrogen-poor phase has a greater amount of hydrogen than the amount of silicate in the hydrogen-rich phase. An asymmetric regular solution model captures our results (see Methods for further details).
Fig. 2: Phase diagram and energetics of solution. a, Composition of the two coexisting phases from our simulations (symbols) near 2 GPa (red), 5 GPa (grey) and 10 GPa (blue) with uncertainties in composition indicated (Extended Data Table 1). Curves with corresponding colours are our computed phase equilibria with estimated uncertainty indicated by shading. Two experimental measurements of the hydrogen solubility in mafic melt18 at 2 GPa (red squares with error bars) are also shown. b, Excess properties of solution at 5 GPa and 4,000 K from our homogeneous simulations: enthalpy (kJ mol−1) (solid line and black symbols), entropy (J mol−1 K−1) (long dashed line and dark grey symbols) and volume (cm3 mol−1) (short dashed line and light grey symbols). The total Gibbs free energy of solution (thin solid line) is also shown. Curves in a and b are computed from the same asymmetric regular solution model (Methods). Full size image
Homogeneous simulations show that immiscibility between hydrogen-poor and hydrogen-rich fluids originates in a positive enthalpy of solution (Fig. 2b and Extended Data Table 2). The excess entropy of solution far exceeds the ideal entropy of mixing (5.8 J mol−1 K−1), producing solvus curves (binodes) that are flat-topped: that is, the compositions of the two coexisting phases diverge rapidly from each other on cooling below T C . The negative volume of solution causes T C to decrease with increasing pressure: from 3,500 K at 2 GPa to 2,600 K at 10 GPa (Fig. 3).
... continue reading