Abstract

Sub-Neptunes and super-Earths, the most abundant types of planet in the galaxy, are unlike anything in the Solar System, with radii between those of Earth and Neptune1,2. Fundamental questions remain regarding their structure and origin. Although super-Earths have a rocky composition3, sub-Neptunes form a distinct population at larger radii and are thought to consist of a rocky core overlain by a hydrogen-rich envelope4,5. At the extreme conditions of the core–envelope interface (exceeding several gigapascals and several thousand kelvin4,6), reaction between core and envelope seems possible, but the nature and extent of these reactions are unknown. Here we use first-principles molecular dynamics driven by density functional theory to show that silicate and hydrogen are completely miscible over a wide range of plausible core–envelope pressure–temperature conditions. We find the origin of miscibility in extensive chemical reaction between hydrogen and silicate, producing silane, SiO and water species, which may be observable with ongoing or future missions. Core–envelope miscibility profoundly affects the evolution of sub-Neptunes and super-Earths, by dissolving a large fraction of the hydrogen of the planet in the core and driving exchange of hydrogen between core and envelope as the planet evolves.

Main

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 MgSiO3 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 MgSiO3–H2 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).

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 MgSiO3 stoichiometry, that is, MgSiO3 behaves as a component, and the system is nearly binary. We choose as our compositional variable the molar fraction of the H2 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 TC. 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).

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 TC. The negative volume of solution causes TC to decrease with increasing pressure: from 3,500 K at 2 GPa to 2,600 K at 10 GPa (Fig. 3).

Fig. 3: Phase diagram and core–envelope interface conditions.

The equilibrated conditions of our simulations: black circles represent the conditions at which two phases coexist; white circles represent the conditions at which a single homogeneous phase is stable (uncertainties in pressure and temperature are smaller than the size of the symbols). The black line is the critical curve TC(P) and is computed from the same asymmetric regular solution model used to generate the curves in Fig. 2 (see Methods); the shading provides an estimate of the uncertainty in TC. The black dashed line is the MgSiO3 liquidus59. Our results are compared with the conditions at the core–envelope interface in models of sub-Neptunes from the literature (red symbols and dashed lines): squares29 refer to rocky planets accreting by pebbles from 0.75M⊕ to 2.5M⊕ within a stellar disk at 1 au; down triangles15 represent Kepler-36c, outer core–envelope boundary at the beginning of the isolation stage (hotter) and at 7 Gyr; diamonds36 represent a sub-Neptune with M = 4M⊕ at the end of the accretion stage (hotter) and at the end of the mass loss stage when the envelope mass Menv/M = 0.025; red dashed line with no symbols30 represents rocky planets embedded in the disk and accreting by planetesimal impact from 0.1M⊕ (colder) to 3.8M⊕; up triangles16 denote a sub-Neptune with mass M = 4M⊕ and hydrogen envelope mass Menv/M = 0.02 on cooling from 1 Myr (hotter) to 10 Gyr in decade increments. The red solid line shows the conditions at the core–envelope interface in sub-Neptunes with M = 1–8M⊕, Menv/M = 0.02, equilibrium temperature Teq = 500 K and intrinsic temperature Tint = 70 K, with the crosses marking unit increments in M/M⊕, following existing methodology16,36.

The behaviour that we find was not anticipated by previous studies. Estimates of the solubility of hydrogen in silicate melt24,25,26,27 used an expression based on the solubility of noble gases (Henry’s law)28. This expression cannot capture complete miscibility because it assumes that hydrogen is inert. For example, one parameterization24 yields 20 mol% solubility of H2 in silicate melt at 2 GPa, 4,000 K, and the other two studies yield smaller solubilities, as compared with the unlimited solubility (that is, miscibility) that we find. Estimates of the solubility of silicate in hydrogen fluid15,29,30,31,32 assumed that silicate solubility is limited by congruent silicate vaporization and ideal mixing, leading to underestimated solubility: for example, 0.2 mol% of silicate in hydrogen at 2 GPa, 4,000 K, as compared with the unlimited solubility that we find. A previous study19 argued that TC should be greater than or equal to the congruent vapour–liquid critical temperature Tsilicate, and that by analogy with the water–hydrogen system, TC might increase with increasing pressure. However, we find TC (3,500 K at 2 GPa) to be substantially less than Tsilicate (6,600 K at 0.14 GPa) (ref. 20) and TC decreasing with increasing pressure, illustrating the importance of fully non-ideal interactions as in our simulations. Complete miscibility in the fluid silicate–water system is well known33,34,35, and partly by analogy, a recent study speculated that miscibility might occur in the silicate–hydrogen system as well36. However, the silicate–water system is fundamentally different as reaction does not require redox exchange, unlike the silicate–hydrogen system.

We probe the microscopic mechanism of reaction between the two phases by examining the evolution of species during equilibration (Fig. 4). Initially, phase coexistence simulations have all H as H2, all O bonded to Si and Si surrounded by four O on average with no SiO. As the simulation proceeds (towards a homogeneous state in the example shown), H2 diminishes at the expense of hydroxyl and molecular water species, and SiO and SiHn increase in abundance. These changes are captured by the redox reactions

$${\text{SiO}}_{2}+{\text{H}}_{2}\to \text{SiO}+{\text{H}}_{2}\text{O}$$

(2)

$${\text{SiO}}_{2}+3{\text{H}}_{2}\to {\text{SiH}}_{4}+2\text{OH}$$

(3)

producing reduced Si species SiO and SiH4, while oxidizing H. We confirm the presence of O–H and Si–H bonding in our simulations by computing the radial distribution function37 over the equilibrated, homogeneous portion of the simulation, which shows clear peaks at short distance (<2 Å) in the O–H and Si–H radial distribution functions (Fig. 4). For illustration, we have written equation (3) in terms of SiH4, although in the simulation shown in Fig. 4 SiH4 is subordinate to lower coordinate species SiHn (n < 4). However, we find that silane becomes more abundant than lower coordinate SiHn species at greater hydrogen concentration, for example, in the hydrogen-rich fluid in phase coexistence simulations performed at temperatures less than TC. Trace amounts of silane and water have been identified as reaction products in experiments in the MgSiO3–H2 system, albeit at much lower temperatures than our simulations and without the ability to quantify abundances38. Apart from the redox reactions, the relative amounts of OH and H2O are governed by H2O + O → 2OH, where the O is bonded to Si (ref. 39).

Fig. 4: Speciation in the miscible fluid.

Evolution of the speciation in a phase coexistence simulation at 5 GPa and 4,000 K that ultimately becomes homogeneous (inset) and the structure of the equilibrated homogeneous fluid (main figure). The main figure shows the H– and O– partial radial distribution functions and the inset shows the proportions of H atoms that are bound to one other H atom (H2) (black), the proportion of O atoms that are bound to one (OH) or two H (H2O) atoms (red), and the proportion of Si atoms that are bound to a single O atom (SiO) (blue), or to any number of H atoms (SiHn) (green).

The products of reactions (2) and (3) may be observable by the James Webb Space Telescope (JWST). H2O has been observed in the atmosphere of sub-Neptunes at concentrations greater than stellar metallicity5,40,41,42. Determining whether the observed water enrichment is endogenous, that is, due to core–envelope reactions, such as equation (2), or from accreted ice from beyond the snow line43,44,45, is central to understanding the origin of these planets46. Some previous planetary modelling studies have proposed that water is produced primarily by reaction of hydrogen with FeO (ref. 26), but reaction of hydrogen with SiO2 (equation (2)) may be a more important source of water27 because SiO2 is expected to be much more abundant in rocky cores than FeO (ref. 47). Previous studies have assumed ideal mixing between species from the low pressures to which observations are most sensitive (10−6 GPa) all the way down to the core–envelope boundary (5 GPa) (refs. 26,27,48), and it will be important in future to include the strongly non-ideal effects that we find to evaluate whether some or all of the observed water enrichment is endogenous. Silicon species have not yet been identified in the atmospheres of sub-Neptunes, but the possibility of observing SiO with JWST has been evaluated in detail49, and silane has been observed in a brown dwarf with JWST50. Planetary modelling studies have examined whether core–envelope reactions such as equations (2) and (3) can produce silicon species at detectable levels high in the atmosphere27,32,51,52. Among the important considerations for future studies will be to include the non-ideal thermodynamics that we find and the incongruent saturation and condensation of silicon species.

Our results indicate that the miscibility of core and envelope must be considered in models of the structure and thermal evolution of sub-Neptunes (Fig. 3). The reason is that probable pressure–temperature conditions at the core–envelope boundary lie well above the critical curve, in the regime of complete miscibility between silicate and hydrogen. For example, in a thermal evolution model of a sub-Neptune that assumes no reaction between core and envelope, the conditions at the core–envelope boundary lie above the critical curve throughout the first billion years of evolution (for example, 4.5 GPa, 3,600 K at 1 Gyr) (ref. 16). This model is, therefore, not in chemical equilibrium. Earlier states of evolution are likely to be even hotter and lie even farther above the critical curve. For example, the core–envelope boundary lies at 4 GPa and 13,000 K in a recent model at the end of the accretion phase36, and even hotter conditions are found during accretion15,30. Although the precise conditions at the core–envelope boundary depend on many factors, including the mass of the rocky core, how much gas the core accretes and how much it retains after disk dispersal, the pressure–temperature conditions that we have highlighted are typical of a wide variety of sub-Neptune models at various stages of evolution6,16. The only sub-Neptunes that do not exhibit H2–MgSiO3 miscibility are less massive (<4M⊕), older (≥10 Gyr) or with lower envelope mass fractions (<0.02) than most of the sub-Neptune population.

Our results have important implications for every stage of the evolution of sub-Neptunes, from the accretion phase to boil off, to present-day observations of atmospheric chemistry. Consider first the accretion phase in which a sub-Neptune forms by accreting a rocky core while embedded in the disk6,53. The rocky core gravitationally captures a portion of the disk, forming a hydrogen-rich envelope. Previous analyses of this stage assume limited dissolution of hydrogen in the core, or silicate dissolution in the envelope dictated by the congruent vapour pressure of the silicate, while neglecting hydrogen solubility in the silicate15,24,25,29,30,36. Our results indicate that these are not good assumptions. If the planet maintains chemical equilibrium at the interface between the envelope and the core, then silicate and hydrogen are completely miscible over a wide range of plausible evolutionary scenarios (Fig. 3).

As the disk disperses, the planet loses mass by evaporation. If the envelope and core are chemically distinct, most of the envelope mass is lost during this phase: for example, in one scenario54 85% of the envelope is lost to boil off. However, if much of the hydrogen envelope is dissolved in the core, then a greater fraction of hydrogen may be retained by the planet after disk dispersal. A more in-depth analysis of this problem would consider the complete miscibility that we find, along with variations of the pressure–temperature conditions at the core–envelope interface that occur as a result of boil-off, and the attendant changes in the phase equilibria.

After boil-off, the planet cools, altering the pressure and temperature conditions at the core–envelope interface. As the equilibrium between hydrogen-rich and hydrogen-poor phases evolves on cooling, partitioning of hydrogen between the core and envelope alters, leading to changes in the mass of the envelope and the planetary radius. Some sub-Neptunes may cross from earlier hotter states that are super-critical to later colder states that are sub-critical (Fig. 3), leading to changes in the composition of core and envelope that may involve silicate condensation, rain-out55, cloud formation or regions of convective inhibition19,56. Some sub-Neptunes may lose their envelopes because of photo-evaporation or core-powered mass loss, forming super-Earths8,9. The formation of super-Earths may also, therefore, be affected by silicate–hydrogen miscibility: if much of the hydrogen of sub-Neptunes is dissolved in the rocky core, the hydrogen may be less susceptible to loss24. Some super-Earths may have small low-molecular-weight atmospheres that originate in core–envelope reactions57, and even those that have completely lost their envelopes may still bear the imprint of chemical reaction between core and the envelope, including sequestration of hydrogen to a metallic core, and the production of endogenous water14,58.

Methods

Phase coexistence method simulations

Our simulations follow the phase coexistence method described in our previous work20,21,60. First-principles molecular dynamics simulations are based on density functional theory in the PBEsol approximation61 using the projector augmented wave method as implemented in VASP (Vienna Ab initio Simulation Package)62,63. We assume thermal equilibrium between ions and electrons using the Mermin functional64. Sampling the Brillouin zone at the Gamma point and a basis set energy cutoff of 500 eV converges the total energy and pressure to within 3 meV per atom and 0.2 GPa, respectively. The outermost core radii (in Bohr) and number of electrons treated as valence are as follows: H (1.1,1), Mg (2.0,2), Si (1.9,4) and O (1.52,4); the simulations are non-spin polarized.

We perform phase coexistence simulations in the canonical ensemble with a Nosé–Hoover thermostat65,66 with a time step of 0.1–1.0 fs for 10–70 ps, sufficient to achieve equilibration, and with drift of the conserved quantity not exceeding 1.5 meV atom−1 ps−1. For the longest time steps (1.0 fs), we used a simple form of hydrogen mass rescaling67, setting the mass of the hydrogen atom to 4 u or 16 u. We prepare the initial condition by first performing a homogeneous simulation of MgSiO3 fluid at temperature T in a cubic simulation cell containing 135 atoms, adjusting the dimension L until the pressure P matches the desired pressure (2 GPa, 5 GPa or 10 GPa). We then perform a homogeneous simulation of H fluid at T in a simulation cell of identical dimensions and adjust the number of H atoms until the pressure P matches that of the MgSiO3 simulation. The number of H atoms, therefore, varies depending on P and T, reflecting the contrast in compressibility and thermal expansivity between MgSiO3 and H fluid. We form the initial configuration of the phase coexistence simulation by combining equilibrated snapshots of MgSiO3 and H fluid, forming a simulation cell of dimensions L × L × 2L.

We quantify the compositions of the two coexisting fluids in two ways20,21. First, we use the one-dimensional density profile normal to the interface, which we find follows the expected hyperbolic tangent form[22](https://www.nature.com/articles/s41586-025-09970-4#ref-CR22 “Widom, B. Potential-distribution theory and the statistical-mechanics of