Miscibility of rock and ice in the interiors of water worlds

Super-Earths and sub-Neptunes are the most common planet types in our galaxy. A subset of these planets is predicted to be water worlds, bodies that are rich in water and poor in hydrogen gas. The interior structures of water worlds have been assumed to consist of water surrounding a rocky mantle and iron core. In small planets, water and rock form distinct layers with limited incorporation of water into silicate phases, but these materials may interact differently during the growth and evolution of water worlds due to greater interior pressures and temperatures. Here, we use density functional molecular dynamics (DFT-MD) simulations to study the miscibility and interactions of enstatite (MgSiO3), a major end-member silicate phase, and water (H2O) at extreme conditions in water world interiors. We explore pressures ranging from 30 to 120 GPa and temperatures from 500 to 8000 K. Our results demonstrate that enstatite and water are miscible in all proportions if the temperature exceeds the melting point of MgSiO3. Furthermore, we performed smoothed particle hydrodynamics simulations to demonstrate that the conditions necessary for rock-water miscibility are reached during giant impacts between water-rich bodies of 0.7–4.7 Earth masses. Our simulations lead to water worlds that include a mixed layer of rock and water.

. Isochores for systems 1-3 (starting from left to right) are plotted in P-T space. The blue symbols represent conditions under which rock and water remained separate while the red symbols indicate a fully mixed rock-water fluid. Our simulations predict that mixing occurs as soon as the rock melts. The T-P conditions agree with the experimental melting line of MgSiO 3 59 . With increasing temperature, water transitions from a solid (blue hexagons) to a superionic phase (blue crosses) and finally into a liquid state (blue triangles). The solid-superionic transition 28 and H 2 O melt curve 81 are plotted. The pv-ppv phase transition was obtained from 82 . The transition to the ultra-high pressure phase α-Mg 2 SiO 5 H 2 from Ref. 38 is also plotted. The grey dashed line shows a model for a young Uranus planet at an age of 0. www.nature.com/scientificreports/ miscibility of high-pressure phases of MgSiO 3 and H 2 O has not been studied experimentally and the interactions of these materials, are assumed, but unknown at extreme P-T conditions. In earlier works, DFT-MD has been used to determine miscibility of planetary material such as water in hydrogen and helium 14 pertinent to Jupiter and of iron and MgO for terrestrial planets 13 . Therefore, we are motivated to use DFT-MD to evaluate a larger P-T landscape and extend petrological studies 6,18,43 of rock and water interactions providing an upper bound for the miscibility curve of these mixtures at high pressures.

Results
Results from atomistic simulations. We simulated MgSiO 3 -H 2 O mixtures at conditions relevant to rock-water interfaces within water worlds. Three systems were simulated at different temperatures while keeping their densities constant. System 1: we prepared a superlattice of bridgmanite (MgSiO 3 ) and ice VIII (H 2 O) at a density of 3.45 g cm −3 , corresponding to a pressure of 30 GPa at 0 K. System 2: we considered bridgmanite and ice X at a density of 3.85 g cm −3 (60 GPa); and system 3: post-perovskite and ice X at 4.38 g cm −3 (120 GPa) (see Fig. S5 in the supplementary material). Since the lattice constant mismatch between the MgSiO 3 and H 2 O crystals is unavoidable, we strained the lattice constants of both structures to fit them in a single, periodic supercell. The relative change in the lattice constants of each crystal was less than 2.5%. With our rock-water systems built, we step-wise increased the temperature from one DFT-MD simulation to the next in 500 K increments. Water transitioned through three phases during our calculations. Initially, it remained in the solid ice phase during the first couple temperature increments depicted by the blue hexagon symbols in Fig. 1. At temperatures beginning at ∼ 1500 K, water transitioned into the superionic phase where the oxygen atoms remain fixed in their lattice positions and the hydrogen atoms diffuse through that lattice 28 . We do not observe hydrogen atoms diffusing into the rock while the rock is solid and water is superionic (' × ' symbols in Fig. 1). At even higher temperatures, water melts completely while the rock remains solid and does not mix with water.
The temperature at which two materials become miscible can be determined with at least two different methods, in the same way, that melting temperatures of a material can be inferred: (1) Most simply, we can proceed in a similar manner to the heat-until-it-melts approach 44-47 by performing a series of MD simulations at a fixed volume while increasing the temperature step by step until we observe mixing, which we will call heat-until-itmixes approach. While the heat-until-melts approach is prone to overestimating the melting temperatures due to overheating of defect-free samples in computer simulations, our heat-until-it-mixes approach should not lead to a significant overestimation of this miscibility temperature because there are already a solid and a liquid phase present in our simulations that resembles the two-phase method 48 , which is regarded as an accurate method to obtain melting temperatures. Still, when the sample melts, the heat-until-it-melts method yields an upper bound to the melting point, and so it does our heat-until-it-mixes approach for the miscibility temperature. One should also note that there are cases when the heat-until-melts approach works very well. For example, the melting temperatures of beryllium predicted in this way are in good agreement with two-phase and free energy calculations 49 , and the solid-to-superionic transitions in water ice 28 can be determined very accurately, where the transition can be reversed within a small temperature interval. (2) Alternatively, the melting temperature can be determined by computing the Gibbs free energies of both phases for a given pressure and temperature 28,[49][50][51][52][53] . This is a much more complex process that typically involves some form of thermodynamic integration 13,14,54,55 . Similarly, one could determine the P-T conditions at which rock and water become miscible by comparing the Gibbs free energy of the mixture with that of the separate components.
In this paper, we determine the miscibility by pursuing the simpler heat-until-it-mixes approach. We expect this method to yield a reasonable upper bound of the true miscibility temperature because in our simulations, there are already a fluid and a solid phase present when the system transitions into fully mixed fluid state. We conducted our DFT-MD simulations step-wise from 500 to 8000 K in 500 K increments. If the two different materials spontaneously start to mix during our simulations, we know the materials are miscible at this temperature. With this comparatively simple method, we obtained novel results relevant to planetary science.
The most direct way to detect the melting and mixing of materials is by following the motion of the atoms throughout the trajectory. Once equilibrated, the mean squared displacement (MSD) was used to determine the onset of superionicity for water and distinguish between solid and liquid phases. However, while the MSD helps determine atomic diffusion, it does not show whether two materials are miscible because, while necessary, it is not a sufficient condition to determine homogeneous mixing. For example, once the water ice goes superionic, the MSD of hydrogen (H) atoms increases as a function of time which indicates diffusion, but does not indicate that there is mixing of rock and water. Similarly, if one layer melts, all of its atomic species will begin diffusing, the MSD values will increase over time, but the two fluid phases may remain separate. On the other hand, diffusion across the rock-water interface for all atomic species is a direct way to detect mixing ( Fig. 2A, B). Therefore, we followed the motions of each atom in the direction perpendicular to the interface illustrated in Fig. 2C (and Figs. S2 and S3) to determine the mixing of material across the rock-water interface. It is straightforward to identify the transition to the mixed-phase with this method, and we employed it to derive the thermodynamic conditions that define the boundary between the miscible and immiscible regimes in Fig. 1.
Another method used to demonstrate material mixing is plotting the atomic radial distribution function (RDF). Figure 2D and E show how the Mg and Si nuclei are correlated with the oxygen atoms in the rock and water phases. For example, in the unmixed state at 5000 K (dashed curves in Fig. 2D and E) the Mg and Si atoms are strongly correlated with oxygen atoms in MgSiO 3 , while they are less correlated with the oxygen atoms of H 2 O, shown by the blue dashed curve. These correlations are similar to those found in silicate melts 37 and pure MgO at much higher pressures and temperatures 56 , where a peak in the Mg-O g(r) function around 1.8 Å was also found for 2-fold compressed MgO and 5-fold compressed Mg 57 , which is larger than the Mg-O coordination However, under no conditions did we find a stable molecular bond between the Mg and O species when MgSiO 3 is liquid. In the mixed state at 7000 K, the distinction between the groups of oxygen atoms disappears, and both types of oxygen atoms are equally correlated with the Mg and Si nuclei in a fluid mixture (solid lines). We see in Fig. 1 at temperatures above the MgSiO 3 melt curve 59 , rock and water fully become miscible. Therefore, we find that the MgSiO 3 melt curve is a lower bound for rock-water miscibility, but because we used the heat-until-it-mixes method, we may be overestimating the temperatures for miscibility. The melting of water before the rock should also result in depressed melting temperatures for rock and therefore lower the temperature for mixing. Nevertheless, using the heat-until-it-mixes approach, we show that rock-water miscibility occurs when rock melts and allows us to determine the conditions at which mixing occurs deep within water worlds. Therefore, these results suggest that these chemical interactions should be considered in models of planetary interiors.
Results from impact simulations. During the rapid growth of planets, the combination of radiogenic heating and the heat of accretion leads to melting and separation into compositional layers. Nevertheless, these stratified layers can mix depending on the heat budget of the planet and the conditions reached within the interior, which can be influenced by giant impacts that deposit substantial energy deep inside a growing planet 60 . These major thermodynamic events transiently lead to a redistribution of material in an already differentiated body. Here, we consider whether collisions between water world planets could generate miscible layers of rock and water. We considered an exploratory set of accretionary giant impacts between 0.7 and 4.7 M ⊕ mass planets with varying water fractions using a smoothed particle hydrodynamics (SPH) code (Tables S2 and S3). The final bodies have masses (3.8-5.3 M ⊕ ) that fall near the observed transition in exoplanet mass-radius trends, around 1.5-2 Earth radii. In SPH codes, the continuum is represented by a large number of particles of fixed mass and composition, represented by an equation of state. Each body contained three distinct layers: metal core, silicate mantle, and pure water. As shown in Fig. 3, the material layers are shocked and mechanically mixed during a giant impact. Depending on the impact parameters and initial interior profiles, different amounts of each material are transiently mixed together. The treatment of materials in SPH does not enable consideration of miscibility, so during gravitational re-equilibration, the phases separate by buoyancy which tends to re-differentiate the merged www.nature.com/scientificreports/ body into strongly thermally stratified layers, as found in collisions between terrestrial planets with rock and iron layers 60 . In addition, some material is emplaced into orbit, and some escape the system. In general, during graze-and-merge collisions, the compositional layers are heated and partially mixed during the event (Movies S1 and S2). After gravitational equilibration, the densities of the most shock-heated material may be lowered to a point where it is stable within an outer material layer, leading to pressure levels with mixed composition in the postimpact body, as shown in Figs. 4A and S5A. Here, super-heated metal is mixed into the silicate layer, and superheated silicate is mixed into the water layer. The temperatures of these buoyant mantle particles are often much greater than the adjacent water particles, which correspond to the least shock-heated water (Figs. 4B, S5B). In an SPH calculation, this thermal stratification persists because of the absence of miscibility and heat transport. In a natural giant impact, thermal and chemical equilibration can occur within these mixed layers (brown shading in Fig. 4). When mixing occurs under conditions that exceed the miscibility boundary, the post-impact planet's chemical profile contains a new compositional zone with a rock-water mixture.
In general, during accretionary giant impacts, we find that the outer portion of the silicate mantle is heated far above the melt curve temperature (Table S3) and partially mixed into the overlying water layer. For the cases explored here, 6-22% of the mantle material is mixed into the water layer at the end of the event. In each case, a greater portion of silicate experienced transient mixing with water during the collision, and in some cases, these transient contacts exceeded the miscibility threshold but could not be captured in the end state due to the limitations of SPH. In our calculations, the lower portions of the water layer sometimes reach the point of melting superionic ice, but exceeding this melt curve depends on the specifics of the initial planet (thermal profile and water layer thickness) and the impact conditions. In general, our results support the idea that when giant impacts occur between water worlds, conditions at the top of the rocky mantle exceed the conditions for miscibility. Even though much of the impact energy is deposited into the outer water layers, a large fraction of the internal energy increase is deposited in the deeper mantle materials. We found the internal energy increase of the mantles to correspond to 16-44% of the initial kinetic energy of the event (Table S3). How much water  www.nature.com/scientificreports/ can dissolve into a mixed rock-water layer depends on the details of the event. Our results suggest that the mixed fraction of the rocky mantle can be substantial, on the order of 10% or more.

Discussion
Previous works, e.g. 15,17 , have considered the dissolution of water into magma oceans on exoplanets by primarily focusing on lower pressure regimes found closer to the surface and rock-water interactions with an atmosphere rich in H 2 O. Instead, we investigated rock-water miscibility at greater depths, or extreme pressures and temperatures, during the collisional growth of planets with a massive water layers overlaying a rocky interior. After giant impacts between water worlds, we expect their thermal evolution to be significantly influenced by the presence of a mixed, stably stratified compositional layer that evolves through secular cooling. In Fig. 1, we also plot the adiabatic profile for a model 'young' Uranus with an age of 0.1 Gyr 8 . A substantial fraction of the heat from formation is still retained at this age, and the predicted P-T profile passes through our miscibility regime for rock and water. If ice-giant-type planets formed with such high interior temperatures (or experienced giant collisions as shown in Fig. 3), they should include an extended, mixed layer of rock and water in their interiors. Only after sufficient cooling will the rocky and icy components separate. The mixed layer may therefore be stratified and stable against convection 61 , delaying the planet's cooling. Therefore, characterizing the miscibility properties of different planetary materials is vital in developing planetary evolution models across various planet types. As we discover and relearn how material interacts at extreme conditions, we need to shift our assumptions of fully differentiated planetary interiors to dynamic ones with 'fuzzy' interior layers where complex chemical interactions occur.
We studied the miscibility of rocky and icy material under conditions pertinent to water world interiors by performing DFT-MD simulations of bridgmanite and water mixtures at high pressure and temperatures. Under the conditions investigated, MgSiO 3 exists as either solid or liquid 59,62 while H 2 O can be solid, superionic, or liquid 28 . We find that rock and water are miscible and form a homogeneous mixture once MgSiO 3 melts. We computed the equation of state of the mixture and the validity of the AVL based on the equation of state of liquid MgSiO 3 and liquid H 2 O at the same pressure and temperature conditions. Our findings indicate the predicted densities obtained from AVL are in excellent agreement with the density of the fully interacting mixture. For example, at 7000 K the AVL overpredicts density by only 0.3%.
We also performed dynamic smoothed particle hydrodynamics simulations of two colliding planetary bodies with a three-layer structure of metal, silicate, and water. We show that all three materials are sufficiently heated during impact events and become dynamically mixed. Therefore we provide proof for the presence of a mixed rock-water layer that is stably stratified and provided this layer remains sufficiently hot; it may persist for long periods during secular cooling. Since giant impacts are a ubiquitous stage of planet formation, our work suggests that the formation of miscible rock-water layers may be common in water-rich exoplanets. Studies of the internal properties of water worlds should consider the effects of mixing and demixing on their thermal evolution, and assumptions on non-ideal mixing should be revisited 63 . This work ultimately provides theoretical evidence that many water-rich exoplanets have mixed mantles. , and H(ultrasoft). A canonical ensemble (constant number of particles N, volume V, and temperature T) was used and temperature was regulated with the the Nosé-Hoover thermostat 67,68 . Exchange-correlation effects were modeled by the Perdew, Burke, and Ernzerhof functional 69 . All molecular dynamics simulations used a 0.2 fs timestep. Electronic wave functions were expanded in a plane-wave basis with an energy cut-off of 1400 eV and we used the Ŵ-point to sample the Brillouin zone of our supercells for all simulations.
The supercells consisted of MgSiO 3 and H 2 O merged along the [001] direction (see Figs. 2A and S6). The crystal phases used for MgSiO 3 were bridgmanite (30 GPa and 60 GPa) and post-perovskite (120 GPa). The phases for water used were ice VIII (30 GPa) and ice X (60 and 120 GPa). System 1 (bridgmanite-ice VIII) refers to the initial 0 K configuration of the rock-water cell at 30 GPa, system 2 (bridgmanite-ice X) the 0 K configuration of the rock-water cell at 60 GPa, and finally system 3 (ppv-ice X) being the 0 K configuration of the rock-water cell at 120 GPa. The water to rock ratio is defined as the number of H 2 O formula units over the total number of MgSiO 3 and H 2 O formula units in the system. Therefore, system 1 has 32 MgSiO 3 units and 64 H 2 O units resulting in an water to rock ratio of 0.26. System 2 and 3 have 32 MgSiO 3 units and 72 H 2 O units leading to a water to rock ratio of 0.29. System size convergence tests were performed to verify that these cells provide sufficiently well-converged results. The details are provided in the supplementary material.
In system 3, the orthorhombic primitive cell for MgSiO 3 was replicated 4 × 1 × 1 , resulting in cell dimensions of 9.9098 × 8.1253 × 6.1493 Å 3 . For the ice X phase, the cubic primitive cell was replicated 4 × 3 × 3 , leading cell lengths of 10.4738 × 7.8554 × 7.8554 Å 3 . Since the lattice constant mismatch between the MgSiO 3 and H 2 O crystals is unavoidable, we strained the lattice constants of both structures to fit them in a single, periodic supercell. The cell dimensions for system 3 resulted in a superlattice of dimensions 10.1920 × 7.9900 × 14.6046 Å 3 . Since the lattice constants were strained, we performed an additional ionic relaxation within the supercell to remove residual stresses. The lattice constants of the supercell were relaxed again to achieve the hydrostatic conditions at a pressure of 120 GPa. The dimensions of this supercell were 9.8191 × 8.0286 × 15.1232 Å. This was used as the initial configuration for the DFT-MD simulations, for system 3 in particular, and we equilibrated it at temperatures ranging between 500 and 8000 K, in 500 K degree steps. www.nature.com/scientificreports/ When the MD trajectories were equilibrated, analysis to determine miscibility was performed. The mean squared displacement (MSD) in the z-direction (perpendicular to the interface), given by where <> i,α denotes a statistical average over all particles of type α and D α the diffusion coefficient, was used to determine the onset of superionicity in the ice and distinguish a solid from liquid. When D α > 0 for all species, the sample is liquid, while D H > 0 with D α = 0 for all other species represents superionic ice. We found that this method helped determine atomic diffusion but could not tell us whether the materials were miscible, as diffusion is a necessary but not sufficient condition to generate a homogeneous mixture.
Therefore, we also rely on the analysis on the individual trajectories of the atoms to confirm that the atomic species diffuse through the system instead of the two separate liquid compounds remaining. We calculated the radial distribution function for each atomic species using, where we calculate the distances between two atomic species of interest α and β throughout a trajectory. δ(r) is the Dirac delta function. The prefactor guarantees that the function equals 1 for noninteracting particles.
SPH. The GADGET-2 smoothed particle hydrodynamics (SPH) code 70 , modified to include tabulated equations of state (EOS) 71 , was used to calculate the structure of isolated bodies and to simulate giant impacts. The giant impact simulations were conducted in a similar manner as in previous work 72 , and the source code is available in the supplementary information in 72 and https:// github. com/ PhilJ Carter/ gadge t2-plane tary. Here, we used an updated version of the Analytic Equations of State (ANEOS) code package that produces a better fit to shock and post-shock temperatures 73,74 . The source code for ANEOS is available in 75 .
For this initial study, we considered water worlds with three distinct compositional layers: metal, silicate, and water. As the properties of water world interiors are largely unknown, the pre-impact bodies were assumed to be analogous to terrestrial bodies with a solid silicate mantle and liquid iron-alloy core. The thermal profile of the mantle was chosen to fall near the solidus of peridotite at the mantle-water boundary, which is a conservatively cool initial condition during the giant impact stage of planet formation. All had a metal core to silicate mass ratio comparable to Earth's value of 0.32. The planets had an exploratory range of initial water entropies, and the total water mass fractions ranged from 7wt% to 50wt%. The metal, silicate, and water layers were modeled using ANEOS models for an iron alloy (Fe 85 Si 15 ) 76 73 but with improvements to the liquid phase. The water ANEOS model is a new parameter set that takes advantage of the triatomic formulation for the gas and improves the model vapor curve compared to previous versions. The ANEOS formulation treats each material as a single component with solid, liquid, vapor and plasma phases. Each of these natural materials has a more complicated phase diagram with multiple solid structures and, in the case of pyrolite, a multi-component composition. The EOS models are developed to capture the bulk response of each material. The ANEOS material parameters are given in Table S1 and further documentation are available at https:// github. com/ stste wart. Non-rotating, isolated bodies were initialized with a specified mass and entropy for each material layer (Table S2). The bodies were then gravitationally equilibrated for 24 hours. At each time step, the entropy of each layer was imposed and the particle velocity was damped, reaching residual velocities of 50-60 cm s −1 .
Here, we considered giant impacts that are typical of accretionary events between protoplanets 78,79 . The most common growth events were graze-and-merge collisions with intermediate impact angles. More grazing impact angles (e.g., 45 • ) lead to hit-and-run events. Simulations were typically run for 48 hours to achieve gravitational re-equilibration, but some required longer times. In most cases, the graze-and-merge events created synestias 80 with massive water-rich disk regions. The impact parameters and post-impact properties for the water world giant impacts are presented in Table S3. For all cases considered, a substantial fraction of the mantle (6-20%) is super-heated and mixed into the outer water layer. Generally, the temperatures of the silicate at the pressures of the mantle-water interface substantially exceed the melt boundary and promotes miscibility with water.

Data availibility
All DFT-MD data are included in the manuscript and in the SI Appendix. SPH simulation initial and final timestep files will be deposited in Dataverse upon acceptance of this manuscript and linked here.