The fate of planetary cores in giant and ice-giant planets

Context. The Juno probe that currently orbits Jupiter measures its gravitational moments with great accuracy. Preliminary results suggest that the core of the planet may be eroded. While great attention has been paid to the material properties of elements constituting the envelope, little is known about those that constitute the core. This situation clutters our interpretation the Juno data and modeling of giant planets and exoplanets in general. Aims. We calculate the high-pressure melting temperatures of three potential components of the cores of giant planets, water, iron, and a simple silicate, MgSiO3, to investigate the state of the deep inner core. Methods. We used ab initio molecular dynamics simulations to calculate the high-pressure melting temperatures of the three potential core components. The planetary adiabats were obtained by solving the hydrostatic equations in a three-layer model adjusted to reproduce the measured gravitational moments. Recently developed ab initio equations of state were used for the envelope and the core. Results. We ﬁnd that the cores of the giant and ice-giant planets of the solar system di ﬀ er because the pressure–temperature conditions encountered in each object correspond to di ﬀ erent regions of the phase diagrams. For Jupiter and Saturn, the results are compatible with a di ﬀ use core and mixing of a signiﬁcant fraction of metallic elements in the envelope, leading to a convective and / or a double-di ﬀ usion regime. We also ﬁnd that their solid cores vary in nature and size throughout the lifetimes of these planets. The solid cores of the two giant planets are not primordial and nucleate and grow as the planets cool. We estimate that the solid core of Jupiter is 3Gyr old and that of Saturn is 1.5Gyr old. The situation is less extreme for Uranus and Neptune, whose cores are only partially melted. Conclusions. To model Jupiter, the time evolution of the interior structure of the giant planets and exoplanets in general, their luminosity, and the evolution of the tidal e ﬀ ects over their lifetimes, the core should be considered as crystallizing and growing rather than gradually mixing into the envelope due to the solubility of its components.


Introduction
The core-accretion model, which assumes that giant planets form by accretion of hydrogen and helium around a solid core, is a reference model in planetary modeling. It is used to explain the rapid formation of large giant planets of several hundred Earth masses by the rapid runaway accretion of gaseous hydrogenhelium material from the planetary nebula if the core size is beyond a critical size (Pollack et al. 1996). In conjunction with basic hydrodynamics arguments, it brings a simplified picture of the interior structures of these planets as two adiabatic layers of varying densities in hydrostatic equilibrium, one for the hydrogen-helium envelope and a second corresponding to a primordial core that is enriched in heavy elements. This simple picture carries over to evolutionary models where the time variation of the luminosity is obtained by integrating the energy dissipated by these layered structures backward in time. The composition of these structures is often assumed to be fixed during the planet lifetime (Guillot et al. 1995;Guillot & Gautier 2014).
The measurements of gravitational moments of the giant planets of the solar system, their luminosity, atmospheric composition, and the mass-radius relationships obtained for the thousands of exoplanets that are now detected provide mounting evidence that this model of planetary interiors needs to be improved (Baraffe et al. 2014;Helled & Guillot 2018). For the giant planets of the solar system, the gravitational moments measured for either Jupiter or Saturn are not compatible with a pure and homogeneous hydrogen-helium envelope (Nettelmann et al. 2012(Nettelmann et al. , 2013. This suggests that in addition to a varying helium concentration, metallic elements such as water or silicates may be present deep in the envelope and close to the core. This view of the Jupiter interior was re-enforced by a recent analysis of the Juno measurements, where a diffuse core was invoked to explain the measured gravitational moments (Wahl et al. 2017;Debras & Chabrier 2019).
The latest formation models, which are more in line with global simulations of the formation of planetary systems, further suggest that late accretion of planetesimals may also explain the varying metallicity that is observed in the atmospheres of giant planets (Zhou & Lin 2007;Alibert et al. 2018). While this approach also leads to a non-uniform density in the envelope, the fate of the primordial core during the planet lifetime currently remains a major unknown in deciding how the underlying hypothesis of the core accretion model might be relaxed. This is by extension also true of the two-or three-layer models. This mostly stems from the lack of accurate physical properties for the elements that may constitute the core at the extreme pressuretemperature conditions that are encountered deep within the envelope and at the core-envelope boundary (Baraffe et al. 2014).
A first step to address this open question was recently obtained by considering the miscibility of potential core constituents within a pure hydrogen plasma (Wilson & Militzer 2012;González-Cataldo et al. 2014;Soubiran & Militzer 2015). This assumes that a solid core is slowly dissolving into the envelope as the planet cools down. These calculations further consider that miscibility for a given element in a pure hydrogen plasma is representative of the miscibility in a planetary envelope consisting of a hydrogen-helium mixture with impurities. This approach consequently neglects the effect of potential additional elements that are also dissolved in the envelope in a significant fraction that may even become dominating as we approach the core. To proceed beyond this step becomes rapidly untractable because the miscibility of all the potential constituents in varying concentrations in the core and envelope needs to be considered. We here suggest that considering the high-pressure melting properties of potential elements that constitute the core is an alternative way to build quantitative interior structure models that are more consistent with their cooling history.

High-pressure melting properties of potential core materials
We considered three basic components that might constitute the primordial core: a simple silicate, MgSiO 3 , water, H 2 O, and iron, Fe. We included iron as a possible core component because dynamical simulations indicate that migration of the planet below the ice line may cause a non-negligible amount of iron to accumulate that in turn may finally rest within the core. We first turn to the high-pressure properties of MgSiO 3 , which presents many polymorphs at low pressures. Because we are mainly concerned with the conditions that are encountered in the core of giant planets, we focus on the post-perovskite (PPV) phase of MgSiO 3 , which is considered as the stable one in the 1−10 Mbar range. Figure 1a shows the equation-of-state (EOS) points calculated using molecular dynamics simulations (see Appendix A). We obtained the pressure dependence of the melting temperature we calculated, T m , by adjusting the semi-empirical Simon law (Poirier 2004) to the ab initio results (see Appendix A). The melting temperature steadily increases with pressure and reaches 18 000 K at 20 Mbar. The liquid or solid states were obtained by considering the mean square displacement once the simulation was equilibrated. When we compare our results with previous estimates that are available up to 4 Mbar (Belonoshko et al. 2005;de Koker & Stixrude 2009), we see that this approach leads to a significant overestimation. This overheating effect is well known. It comes from simulations that are performed at fixed volume while using a limited number of particles. The melting can be overestimated by up to 30%. Because more refined approaches are beyond reach for this system, we applied a conservative coefficient of 0.7 to our estimate (noted 0.7 * T m in Fig. 1a). Figure 1a shows that this allows us to extend the results obtained by Belonoshko et al. (2005) while at the same time satisfying the pressure dependence found in our simulations. Our results further rule out the predictions of de Koker & Stixrude (2009).
Direct inspection of the stress tensor obtained from the simulations also shows that the off-diagonal components become non-negligible for pressures beyond 10 Mbar. This indicates that the PPV phase becomes unstable beyond this pressure range. This result is in line with previous work (Umemoto et al. 2006) that found no stable post-PPV phase for MgSiO 3 and dissociation into simpler compounds beyond 10 Mbar. The current understanding of the MgSiO 3 phase diagram is displayed in Fig. 1b. Since the pioneer work of Umemoto et al. (2006), the dissociation pathway has been refined to include intermediate compounds found in the 10-30 Mbar range. This includes Mg 2 SiO 4 and MgSi 2 O 5 (Umemoto & Wentzcovitch 2011;Wu et al. 2013). Because the details of this dissociation pathway are likely a second-order effect for planetary modeling, we approximated the melting temperature of MgSiO 3 as the combination of the high-pressure melting temperature obtained from our simulations, corrected for overheating to match the results of Belonoshko et al. (2005) at 4 Mbar, with the high-pressure melting temperature of the simpler components, MgO and SiO 2 , beyond 20 Mbar that we obtained previously (Mazevet et al. 2015;Musella et al. 2019).
In Fig. 2a we show the results obtained for dense water. This work complements previous studies performed by French et al. (2009) and Wilson et al. (2013), who reported for the super-ionic phase that the oxygen atoms either lie in FCC or BCC structures. Because different melting temperatures were obtained, we first revisited these calculations by performing direct melting simulations and considered both structures as initial states (see Appendix A). Figure 2a shows that our estimate of the BCC structure agrees well with the result of French et al. (2009). In contrast with the findings of Wilson et al. (2013), we find that the FCC structure is not more stable than the BCC structure above 1 Mbar because the high-pressure melting temperature in the FCC phase is equal to or lower than in the BCC phase. Like French et al. (2016), we also find that the FCC structure appears more stable around Mbar and below this pressure.
This shows that the stability of the BCC structure is a good estimate for planetary modeling. We also find that the BCC structure is unstable just below melting, which indicates that another crystalline structure may be present under these conditions. To further investigate whether super-heating is an issue for the super-ionic state of dense water, we performed two-phase simulations in the BCC structure. We find that simulation of direct melting tends to overestimate the stability of the superionic phase by a few thousand Kelvin. This result is consistent with the fact that the volume change is small through this transition (French et al. 2016). Figure 2a shows our new estimate of the pressure-temperature domain where the super-ionic state should be considered. This reaches 13 500 K at 100 Mbar and includes the FCC phase below 1 Mbar.
We finally turn to the behavior of the last potential element that may be present in planetary embryos: iron. To estimate the stability of the solid phase of iron in the pressure range considered here, we extended our previous calculations (Morard et al. 2011;Bouchet et al. 2013) up to 100 Mbar (see Appendix A). We further performed two-phase simulations by considering the stable phases of iron as predicted by linear response theory (Stixrude 2012). Figure 2b shows that solid iron is found in an HCP state up to the 30-60 Mbar range. The FCC crystalline structure is predicted to be the most stable one beyond this pressure and up to 200 Mbar. Beyond 200 Mbar, the BCC structure is predicted to be the most stable (Stixrude 2012).
Our two-phase simulations confirm this overall result. At 60 Mbar, the high-pressure melting temperature of the FCC L4, page 2 of 6 phase is higher than that of BCC by 5000 K, indicating that it is the most stable phase when both temperature and anharmonic effects are taken into account. We further point out that we used fewer atoms in the simulation cell in the BCC than in the FCC phase. We previously found at lower pressures that this can lead to an underestimation of the melting temperature by a few thousand Kelvin. It is thus likely that the difference between the BCC and FCC melting temperatures reported here is slightly smaller.
Beyond 100 Mbar, we switched to the computationally more efficient Thomas Fermi molecular dynamics approach to converge in the number of particles used. Probably due to the influence of the 3 s state, Fig. 2b shows that the Thomas Fermi regime is not yet fully reached at 100 Mbar because pressures in the BCC phase are overestimated by 15%. Despite this limitation, the method allows for an additional estimate of T m . Figure 2b shows that this method predicts that the BCC phase is more stable with T m up to 3000 K higher than for the FCC phase. This tends to confirm the results obtained using linear response theory. The ab initio results obtained with a smaller number of atoms in the BCC phase do not confirm this finding and tend to suggest that the FCC phase may still be the stable phase. Resolving this question is beyond the scope of this paper, therefore we considered the FCC melting temperature from 30 to 100 Mbar and the Thomas Fermi result beyond 100 Mbar, as indicated in Fig. 2b, to extend the pure iron melting temperature beyond 15 Mbar.

Implications for giant planet cores
We now use these high-pressure melting temperatures we calculated to estimate the state of the core in giant and ice-giant planets. Figure 3 shows the temperature profile obtained for Jupiter when two different models based on different ab initio EOSs (Nettelmann et al. 2012;Militzer & Hubbard 2013) for the envelope and the core were used. The origin of this difference has been well documented elsewhere (Miguel et al. 2016;Mazevet et al. 2019a). We concentrated on the region corresponding to the core. This corresponds to a plateau in the temperature profile because the core is described by an isothermal profile. When we consider Jupiter, the comparison between the predicted interior profiles and the melting temperatures calculated previously indicates that if they are present within the core, MgO and Fe are clearly in a solid state. Figure 3  convect and mix with the envelope and/or establish a regime of layered convection with a significant density gradient due to gravitational settling (Leconte & Chabrier 2013). Figure 3 also shows that the SiO 2 melting temperature is in between the predictions of Nettelmann et al. (2012) and Militzer & Hubbard (2013) for the interior profile of Jupiter. This indicates that depending on the interior model that is used, SiO 2 should either be completely melted up to the center or be present with a solid fraction close to the core center. In both cases, Fig. 3 shows that SiO 2 is melted at the core-envelope interface and can also mix into the envelope and/or participate in a layered convection regime. Considering the low-temperature profile, Fig. 3 shows that the current state of the Jupiter core should be considered as made of solid MgO, possibly including some solid iron, surrounded by solid SiO 2 , some fraction of SiO 2 present at the time of formation potentially mixing into the envelope, and H2O, which can fully mix within the envelope. We point out that this picture provides a stronger physical basis to the latest interpretation of the JUNO data, where a partially diluted core was invoked to reproduce the gravitational moments to high order (Wahl et al. 2017;Debras & Chabrier 2019). This result is obtained by only considering the melting temperature of possible core components and does not rely on their solubility.
For Saturn, our calculations suggest a core containing solid MgSiO 3 together with iron that is either dissolved in silicate or is present as metal. The high-pressure melting temperature predicted from our simulations is however rather close to the isotherm calculated for the core by Guillot & Gautier (2014). The difference between the two is only about 2000 K. This points to uncertainties on the presence of solid MgSiO 3 in the core of Saturn. The situation is clearer for the two ice giants Uranus and Neptune. We find that MgSiO 3 is expected to be in a solid state in their cores. The calculated core conditions are in this case well below the high-pressure melting temperatures. The situation differs for super-ionic water, for which we find that the core conditions for the two ice-giants reported by Guillot & Gautier (2014) lie very close to the melting temperature. It is thus likely that water is also partially mixed into the envelope of ice-giant planets. Overall, the calculations of the melting temperatures performed in this work clearly suggest that the cores of the giant and  Fig. 4. Comparison between the high-pressure melting temperatures obtained with density-temperature profiles of Jupiter (blue lines) and Saturn (green lines) where the surface temperatures at 1 bar, T atm , are increased from current conditions to 300 K. The hotter adiabats correspond to higher surface temperatures. The rectangles indicate the conditions for the core we showed in Fig. 3.
ice-giant planets of the solar system are currently partially mixed with the envelope. They also suggest that a double-diffusive or layered convection regime involving a mixture of the various elements that constitute the initial core is likely taking place for the giant planets and that it does not have a clear interface between a hydrogen-helium envelope and a solid core. To further explore this assumption, we now turn to the time evolution of the cores of the two giant planets as they cool down. We show in Fig. 4 adiabats for Jupiter and Saturn that we obtained by increasing the atmospheric temperature T atm at 1 bar. The calculations were performed using our newly developed ab initio EOSs for hydrogen, helium, water, and MgO Mazevet et al. 2019b;Musella et al. 2019). We used a three-layer model that consists of an H-He envelope, and a core made of H 2 O, and MgO or MgSiO 3 for Jupiter and Saturn, respectively. We solved the standard hydrostatic equations and optimized the fraction of each layer using the theory of figures to the third order to match the radius and the first three gravitational moments. By considering T atm measured by the Galileo probe at a pressure P atm = 1 bar and T atm = 165 K (von Zahn et al. 1998), we obtain a profile for Jupiter in which the temperature of the core is slightly higher than the temperature obtained by Militzer & Hubbard (2013), that is, favoring a low-temperature profile for Jupiter. Additional details can be found in (Mazevet et al. 2019a). Figure 4 shows the interior profiles obtained for Jupiter by increasing T atm from current conditions to 200 K, 250 K, and 300 K while keeping all the other parameters fixed. At T atm = 200 K, SiO 2 is completely melted up to the center of the core. At temperature T atm greater than 250 K, Fig. 4 indicates that the MgO core starts to be melted at the core-envelope interface. The MgO core is no longer completely solid and can mix into the envelope or participate in a double-diffusive or layered convection regime. We assume here that viscosity of the solid at these extreme conditions prevents it from participating in a convective or layered convective mode.
The situation is similar for Saturn. We first point out that at T atm = 140 K, the measured temperature at 1 bar, our threelayer model adjusted to reproduce the gravitational moments predicts a core temperature 2000 K lower than the calculations L4, page 4 of 6 of Guillot & Gautier (2014). This comes from the different H-He EOSs that were used and translates for the case of Saturn the sensitivity to the H-He EOSs found for Jupiter (Mazevet et al. 2019a). Figure 4 shows that under current conditions, the core of Saturn contains solid MgSiO 3 . At the coreenvelope boundary, H 2 O is in a liquid state and can thus mix into the envelope or enter a double-diffusive or layered convection regime. A fraction of H 2 O is in a super-ionic state deeper within the center of the core. This suggests that the core is currently made of super-ionic water and MgSiO 3 that is probably dissociated into solid MgO and solid SiO 2 deep within the core. When T atm is increased to 200 K, MgSiO 3 is no longer in a solid state at the core-envelope boundary, while water can no longer be found in a super-ionic state. As T atm is increased to 250 K and 300 K, Fig. 4 shows that the size of the solid core reduces as the core isotherm crosses the MgSiO 3 high-pressure melting curve at successively higher temperatures.
Using standard evolution models for the giant planets (Guillot et al. 1995), we can directly relate the interior profiles calculated at different T atm to earlier states in the history of the planets. Using the time evolution of T atm calculated in Guillot et al. (1995), we estimate that T atm = 300 K corresponds to Jupiter at 1 Gyr after formation, that is, 3.5 Gyr ago, while T atm = 250 K corresponds to the state of Saturn 1.5 Gyr ago. Because T atm further increases earlier on for both planets, this indicates that the cores were likely completely melted during the early stages after formation, in the same way as for HD209458b, for example. This suggests that soon after formation, the cores of both planets were completely melted and likely entered a doublediffusive or layered convection regime involving all the primary constituents of the core, with a likely significant density gradient due to gravitational settling. Figure 4 also shows that Jupiter and Saturn not only have different cores at the present time, but also followed different evolutionary tracks. As the planets cooled down, MgO first precipitated in the core of Jupiter and no longer contributed to a double-diffusive or layered convection regime, but SiO 2 is only recently doing so. Conversely, the core of Saturn likely formed more recently, experienced first precipitation of MgSiO 3 , and is experiencing super-ionic H 2 O precipitating at the present time.

Discussions and outlook
Our calculations of the high-pressure melting temperatures of several components that might be present in the core of giant planets indicate that the solid cores of these objects are not primordial. After an initial regime where all the constituents likely convected and/or entered a double-diffusive or layered convection regime, the solid fraction evolved throughout the lifetime of the planet with a time-dependent solidification of its components. Although less spectacular, a similar situation with an only partially melted core is suggested for the ice giants. Our calculations thus bring a different perspective to the diffuse core approximation that has been invoked to reproduce the gravitational moments of Jupiter that were measured by the Juno probe. They further suggest that the solubility of the various core constituents in a pure hydrogen plasma may not be relevant. In effect, no clear interface between the core and the H-He interface exists throughout the lifetime of the planets of the solar system, with water acting as a buffer between the core forming and the H-He envelope. We further point out that if melted, layered convection or double-diffusive convection of the core constituents can take place, which might affect the energetics and the evolutionary track of each planet differently. We suggest that this may be particularly the case for Saturn, for which current models struggle to explain the current luminosity, and for giant exoplanets in general. Our calculations indicate that the time evolution of their luminosity should consistently take the evolving state of the core into account. Finally, our results suggest that the varying nature of planetary cores should also be considered when tidal effects of the giant planets of the solar system on their satellites are estimated.