SubNeptune planets are a very common type of planet. They are inferred to harbor a primordial HHe envelope on top of a rocky core which dominates the mass. Here we investigate the longterm consequences of the core properties on the planet massradius relation. We consider the role of various core energy sources resulting from core formation its differentiation its solidification latent heat core contraction and radioactive decay. We divide the evolution of the rocky core into three phases the formation phase which sets the initial conditions the magma ocean phase characterized by rapid heat transport and the solidstate phase where cooling is inefficient. We find that for typical subNeptune planets of 210 M SUBSUB and envelope mass fractions of 0.510 the magma ocean phase lasts several gigayears much longer than for terrestrial planets. The magma ocean phase effectively erases any signs of the initial core thermodynamic state. After solidification the reduced heat flux from the rocky core causes a significant drop in the rocky core surface temperature but its effect on the planet radius is limited. In the long run radioactive heating is the most significant core energy source in our model. Overall the longterm radius uncertainty by core thermal effects is up to 15.
['methods: numerical', 'planets and satellites: composition', 'planets and satellites: interiors', 'planets and satellites: physical evolution', 'Astrophysics - Earth and Planetary Astrophysics']
Contribution of the Core to the Thermal Evolution of SubNeptunes
{'Contribution of the core to the thermal evolution of sub-Neptunes': 'A. Vazan 1 , C. W. Ormel 1 , L. Noack 2 , C. Dominik 1 1  \nAstronomical Institute Anton Pannekoek, University of Amsterdam, The Netherlands 2 Department of Earth Sciences, Institute of Geological Sciences, FU Berlin, Germany', 'ABSTRACT': 'Sub-Neptune planets are a very common type of planets. They are inferred to harbour a primordial (H/He) envelope, on top of a (rocky) core, which dominates the mass. Here, we investigate the long-term consequences of the core properties on the planet mass-radius relation. We consider the role of various core energy sources resulting from core formation, its differentiation, its solidification (latent heat), core contraction and radioactive decay. We divide the evolution of the rocky core into three phases: the formation phase, which sets the initial conditions, the magma ocean phase, characterized by rapid heat transport, and the solid state phase, where cooling is inefficient. We find that for typical subNeptune planets of ∼ 2-10 M ⊕ and envelope mass fractions of 0.5-10% the magma ocean phase lasts several Gyrs, much longer than for terrestrial planets. The magma ocean phase effectively erases any signs of the initial core thermodynamic state. After solidification, the reduced heat flux from the rocky core causes a significant drop in the rocky core surface temperature, but its effect on the planet radius is limited. In the long run, radioactive heating is the most significant core energy source in our model. Overall, the long term radius uncertainty by core thermal effects is up to 15%.  \nSubject headings: Methods: numerical, Planetary systems, Planets and satellites: composition, Planets and satellites: interiors, Planets and satellites: physical evolution', '1. Introduction': "Exoplanets in the range of several Earth masses are very common in our galaxy (Howard et al. 2012; Batalha et al. 2013; Coughlin et al. 2016). Observed mass-radius relation for some of  \nthese planets suggest that some of these planets are bare rocky planets with no envelope, and some are inferred to contain some amount of hydrogen and helium (sub-Neptunes) on top of rock/iron core (hereafter core ). The inferred envelope masses for sub-Neptunes are typically of several percent (Lopez et al. 2012; Jontof-Hutter et al. 2016; Wolfgang & Lopez 2015), i.e., most of the planet mass is in the core. In such a planet the core can perform as an energy reservoir for the envelope.  \nIn contrast to envelope (gases) cooling and contraction, the core radius is not expected to change much during the planetary evolution (Rogers & Seager 2010). However, the heat flux from the core can affect the thermal properties of the low-mass envelope by heating it from below. High enough heat flux from the core can lead to envelope mass loss (Ginzburg et al. 2016, 2017), but even a moderate flux from the core can change the envelope thermal properties. Hence, the thermal evolution of the core indirectly affect the planet radius.  \nPrevious astrophysical studies have accounted for some of the core properties, like the decay of radioactive materials, when modeling the thermal evolution of the envelope (e.g., Lopez & Fortney 2014; Howe & Burrows 2015; Chen & Rogers 2016), assuming an isothermal core that cools in the cooling rate of the envelope. However, in the mass range of Neptune-like planets, the core thermal properties can be significant (Baraffe et al. 2008). A key factor in the core-envelope thermal evolution is thus the timescale on which the core releases its heat; if the cooling is on Gyrs timescale, the observed radius is affected by it (Vazan et al. 2018). In order to study the timescale for the core cooling, its thermal evolution should be modeled explicitly.  \nIn geophysical studies core heat transport is explicitly modelled (Turcotte & Oxburgh 1967; Stevenson et al. 1983). The cooling rate is determined by the properties of high viscosity convection and the resulting conductive boundary layers between convective regions. In these models the surface temperature of the planet is fixed as there is no thick envelope. In contrast, in the case of sub-Neptune planets the surface temperature changes in time due to the envelope's cooling and contraction. The thick, gaseous envelope keeps the temperature (and the pressure) at the core surface higher than in an Earth-like case, affecting the viscosity and the state of the core (solid/liquid). Thus, for sub-Neptunes the heat transport properties of the core depends on the envelope properties. Therefore, both envelope and core should be modeled simultaneously.  \nIn this work we model the thermal evolution of the planet as a whole, center to surface, and study the effects of the underlying thermal parameters on the state of the core, and on the derived planet radius. In section 2 we discuss our model, and define different phases of the core thermal evolution. In section 3 we show our results for sub-Neptune planets of ∼ 2 - 10 M ⊕ with envelope mass fractions of 0.1% - 20%, and examine thermal properties  \nthat lie within the uncertainty of geophysical models. We discuss the sub-Neptunes thermal evolution perspective in section 4, and draw our conclusions in section 5.", '2. The model': 'We calculate the thermal evolution of sub-Neptune planets with a rocky core and a hydrogen-helium envelope on a single structure grid (see Section 2.4). As illustrated in Fig. 1, we divide the evolution of the planet into three phases:  \n- 1. Formation - the conditions of the core and the gaseous envelope as derived from estimates of core formation (see Section 2.1).\n- 2. Magma ocean - vigorous (liquid) core convection and efficient heat transport from the core to the gaseous envelope.\n- 3. Solid state - solidification of the core surface and conductive (reduced) heat transfer from the core to the envelope.  \nThe formation phase, in which the planet assembles in the presence of a gas-rich disk, is not explicitly modeled, but provides the initial conditions for the subsequent disk-free phases. The transition from phase 2 to phase 3 occurs when the core surface temperature drops below the solidification temperature for the surface pressure (see Section 2.2). This transition is approximated as instantaneous. The timing of the transition depends on the cooling of the core and the envelope, and on the core energy sources (see Section 2.3). After solidification (phase 3) core heat transport slows down significantly. We distinguish two end-member scenarios: (3a) conductive core; (3b) convective core with conductive surface boundary layer.  \nIn contrast to previous works, we determine the core thermal evolution from the perspective of the core thermal properties, and not from the envelope evolution. The change in core properties in time is simulated by the three phases. The evolution path is continuous between the phases: at the end of each phase the structural parameters (radius, temperature, density, luminosity and composition for each mass layer) are being used for the first step of the next phase.  \nFig. 1.- The three phases of the core evolution: 1. Formation - the core is hot due to conversion of binding energy to heat; 2. Magma ocean phases - vigorous (liquid) core convection; 3. Solid state - we consider either an entirely conductive core (3a), or a convective core with conductive core surface (3b).  \n<!-- image -->', '2.1. Initial conditions': 'The initial energy content of the planet is determined by its formation. During the formation of a core the gravitational energy from the accumulation of the solid materials partially transforms to thermal energy. The speed of the accumulation process, as well as the material properties, determine which fraction of the gravitational binding energy is locked up in the planet in the form of thermal energy. A fast-forming core may retain large fraction of the initial binding energy and reach high initial temperatures. On the other hand, for a slow core formation a substantial amount of the binding energy is released already during the formation process (by radiation), which results in lower initial temperatures.  \nThe maximal temperature for the core formation can be estimated from the gravitational binding energy:  \nE binding = 3 GM c 2 5 R c (1)  \nThe maximal temperature is achieved when all the binding energy is converted to heat, i.e., E binding = C p M c T max and thus  \nT max ∼ GM c R c C p = 4 . 8 × 10 4 K ( M c 1 M ⊕ )( R c 1 R ⊕ ) -1 ( C p 1 kJ (kg K) -1 ) -1 (2)  \nwhere M c and R c are the core mass and radius, G is the gravitational constant, and C p is rock heat capacity (Guillot et al. 1995). For a 4.5 M ⊕ rocky planet of 1.65 R ⊕ we get T max = 1 . 3 × 10 5 K . This is rather high, since it assumes no radiative or advective losses. But even when the initial temperature will be a fraction of T max , substantial amount of  \nheat is locked in the core, and the rocky part will radiate at a high luminosity. Since the fraction of gravitational binding energy that is left in the core after formation is unknown, we consider here a wide range of fractions between 5%-50% of E binding .  \nThe core is embedded in an initially diluted adiabatic envelope. We assume that the planet envelope is formed by gas accretion during the disk phase (Pollack et al. 1996; Ikoma & Hori 2012), and therefore is composed of hydrogen and helium at solar ratio. The envelope extends initially to the planet Hill radius 2 .', '2.2. Flux from the core': "Heat transport by convection in the core is determined by the Rayleigh number  \nRa = ραg ∆ TD 3 κη , (3)  \nwhich depends on structure properties as gravity ( g ), density ( ρ ), layer thickness ( D ) and the temperature difference within the layer (∆ T ), but also on material properties such as thermal expansion coefficient ( α ), thermal diffusivity ( κ ) and kinematic viscosity ( η ). Convection occurs for Rayleigh numbers greater than a critical value (here we use Ra crit = 450, Morschhauser et al. (2011)). According to Eq. 3 convection is less likely and less vigorous if viscosity η and thermal conductivity κ are high. Between these two, the viscosity is the key factor in regulating the convection, since it is an exponential function of the local physical conditions. Specifically, we approximated the viscosity by (e.g., Karato & Wu 1993; Poirier 2000; Noack et al. 2016) 3 :  \nη ( T ) = η 0 exp [ E 0 R g ( 1 T -1 T 0 )] , (4)  \nwhere reference values are extrapolated from Earth-like composition: viscosity η 0 = 10 21 Pas, temperature T 0 = 1600 K, and activation energy E 0 = 240 kJ/mol; R g is the gas constant.  \nIn the core-envelope-boundary (hereafter CEB) heat is transported by conduction. The thickness of the conductive layer δ is determined by the vigor of the underlying convection  \n(Stevenson et al. 1983):  \nδ = D ( Ra crit Ra ) 1 / 3 (5)  \nAs the core cools and the viscosity increases, the Rayleigh number is lower and thus the conductive boundary layer becomes thicker. Since heat transport in the boundary layer operates by conduction - a diffusive process - core cooling depends on the conductive timescale, which is in the order of τ cond = δ 2 /κ .  \nInitially, the core is in a magma ocean phase - a fluid phase of very low viscosity - due to its high temperature from formation. In this phase the conductive boundary layer is very thin, and thus the conductive timescale is short ( τ cond /lessorequalslant 10 6 yr ). In this phase we assume efficient heat transport from the core to the envelope, and in the core, i.e., we simply ignore the conductive CEB and model the core as adiabatic.  \nThe magma ocean phase continues until at some point the core is cold enough to allow solidification of the core surface. To find this point we compare the CEB temperature with the critical melting temperature for Earth mantle composition, as described in Appendix A1. The magma ocean solidifies from the bottom upward, since the melting temperature of rock increases more strongly with pressure than the increase in temperature of the adiabatic core structure. Therefore, CEB solidification implies that the core can no longer be modeled by low viscosity convection.  \nAfter solidification, the high viscosity of the rock slows down the heat transport in the core. The change in the core heat transport after solidification is a gradual process, which depends on uncertain structure and material properties 4 . For simplicity, we consider here two extreme limits:  \n- (3a) conductive core, where we assume that convection will be suppressed due to the high interior pressure (e.g., van den Berg et al. 2010; Stamenkovi'c et al. 2012).\n- (3b) convecting core, with conductive thermal CEB layer (including the crust) on top of it. We model the transition between the magma ocean phase (phase 2) and the solid state phase (phase 3) as instantaneous.  \nIn the case of an entirely conductive core we take the heat transport in all core layers to be conductive. We set the actual temperature gradient (Eq. 8-9 in Section 2.4 below) to conductive heat transport, using Earth mantle conductivity (Table 2). In the case of conductive CEB layer we assume a conductive CEB layer on top of a convecting core. The thickness of the conductive layer is calculated according to Eq. 5 at the CEB solidification  \nTable 1: Core energy sources  \npoint. We take the heat transport only in that layer to be conductive, while the rest of the core remains convective.", '2.3. Core energy sources': 'Our model includes several energy sources. The magnitude of the involved energies and the timescale on which they operate affect the state of the core and the envelope. Below are the sources included in our model:  \n- I Formation As discussed in Section 2.1, the exact fraction of the accretion energy to be locked in the core depends on formation and thus unknown. Therefore, we calculate evolutionary tracks for planets with different fractions of the binding energy from core formation, in the range of 0.05-0.5 E binding . The lower bound of this range is calculated from the minimal energy of early Earth-like geophysical models (Noack et al. 2017), and the upper bound is an overestimation, assuming half of the impact energy is left in the core. We calculate models for values of 0.05, 0.1 0.2 0.3, 0.5 E binding .\n- II Differentiation The core is initially of low viscosity (molten). Under this condition iron sinks to the center of the planet very efficiently (Stevenson 1990), i.e., during the early evolution. The released gravitational energy further heats the interior. We use the  \nformalism of Solomon (1979) to add the differentiation energy to the core. We assume the iron-to-rock ratio to be Earth-like, and find the differentiation energy to be of a few percent of the binding energy. For simplicity, we add this energy at once (i.e., at t = 0), to the initial energy content of the core 5 .  \n- III Radioactive heating The decay of radioactive elements is an important heat source in rocky planetary interiors (Valencia et al. 2007). The dominant elements with half lives in the Gyrs regime are 238 U , 235 U , 232 Th and 40 K (Nettelmann et al. 2011; Anders & Grevesse 1989). We apply the radiogenic luminosity to Eq. 6 (see below), by using the values of Nettelmann et al. (2011). The abundances of radioactive elements for exoplanets are unknown, but can range between 0.5 and 2.5 times the Earth ratio in solar analogue stars (Unterborn et al. 2015). Therefore, we calculate for planets in this abundance range.\n- IV Solidification (latent heat) As the planet cools the core changes from liquid to solid. The solidification process releases latent heat. We include the latent heat release in our model by adding it to each planetary layer that cools below the solidification temperature. We use the melting curve as provided in Appendix A1. Rock latent heat of 6 × 10 5 J/ kg (Morschhauser et al. 2011) is added to the luminosity (Eq. 6) of each solidified layer.\n- V Core contraction The pressures in sub-Neptune interiors can reach GPa levels, where core compression may heat the core and affect the thermal evolution (Mordasini et al. 2012). In our calculation this effect is automatically included since the heavy element core is part of the evolution structure matrix, and pressure-temperature-density relation are derived from the rock EOS (Vazan et al. 2013). Estimate of the energy from core contraction (by p ∆ V ) is found to be on the order of 10 4 -10 5 J/ kg for sub-Neptune mass range.  \nEstimates of the above energy sources and their release timescales, for a 5 M ⊕ planet with 10% envelope, appear in Table 1. The fits for the energy flux are implemented in the model by using the method described in Vazan et al. (2018). Collecting these effects, the luminosity in the core is taken to be:  \nL core = M c ( c v dT c dt + E radio τ r e ( -t/τ r ) + E solid ∆ t δ ( T -T solid ) ) (6)  \nwhere c v is the specific heat capacity, M c the core mass, and t is time. dT c /dt describes the release of initial energy from formation and differentiation, E radio and τ r are adjusted to fit the heat production by radioactive decay as in Nettelmann et al. (2011), and E solid is the solidification (latent) energy (Morschhauser et al. 2011) released on time interval ∆ t , when the temperature reaches the solidification temperature (T=T solid ). The parameter values appear in Table 2. External core energy sources, such as late planetesimal capture (Chatterjee & Chen 2018), are not included in the model.', '2.4. Thermal evolution': 'The evolution is modelled by a 1D hydrostatic planetary code that solves the structure and evolution equations for the entire planet on one grid (see Kovetz et al. 2009; Vazan et al. 2013, 2015, for details). We use an equation of state (EOS) for hydrogen, helium (Saumon et al. 1995) and rock, as described in Vazan et al. (2013). The thermodynamic properties of the core (such as density, entropy, etc.) are modeled by the EOS of one material (SiO 2 ), as a simplification for the full mineralogy of a core.  \nThe energy balance during the evolution is described by  \n∂u ∂t + p ∂ ∂t 1 ρ = q -∂L ∂m (7)  \nwhere the symbols ρ, p, u are density, pressure and specific energy, respectively, q is the contribution by the additional energy sources, and m,L are the planetary mass and luminosity.  \nThe temperature profile is determined by the heat transport rate according to  \n∂ ln T ∂m = ∇ ∂ ln p ∂m (8)  \nIn convective regions the actual temperature gradient ∇ is the adiabatic temperature gradient ∇ A ; otherwise, heat transfers via radiation (in the envelope) and conduction (in the core), i.e. ∇ = ∇ R , where  \n∇ R = κ op L 4 πcGm p 4 p R . (9)  \nκ op is the harmonic mean of radiative and conductive opacities, and p R is the radiation pressure.', '2.4.1. Atmospheric conditions': "The atmosphere opacity regulates the planet luminosity and the contraction of the en-  \nvelope. We use Rosseland mean of the radiative opacity by Sharp & Burrows (2007) for solar metallicity grain-free atmosphere. Envelopes of sub-Neptunes may have an enhanced atmospheric metallicity from their formation (Fortney et al. 2013; Thorngren et al. 2016). Moreover, extended clouds in sub-Neptune atmospheres (e.g., Bean et al. 2011; D'esert et al. 2011) indicate large atmospheric metallicities and a large opacity (Morley et al. 2013). Because of the important role of atmospheric opacity on the cooling rate of planets, we tested the sensitivity of our results to higher (30 × solar) atmospheric metallicity.  \nAt the atmospheric boundary, which is taken to be the planetary photosphere, the outer boundary condition is κ op p = τ s g where g = GM/R 2 is the gravitational acceleration, and τ s is the optical depth of the photosphere. The envelope mass is constant during the evolution, i.e., no evaporation or gas accretion are included in the model.", '2.4.2. Irradiation': "We assume the temperature distribution in a gray, plane-parallel atmosphere, with constant net outward flux F and irradiation temperature T irr , to be  \nσT 4 ( τ ) = ( 3 4 τ + 1 4 ) F + g ( τ ) σT 4 irr . (10)  \nwhere g ( τ ) = 3 2 (1 -1 2 e -τ ) (Kovetz et al. 1988) and σ is Stefan-Boltzmann's constant. The temperature distribution is calculated for a vertical (maximal) irradiation flux, with no angle dependency of the incident flux (Guillot 2010). At the photosphere, where τ = τ S = 1 , the net outward luminosity of the planet is  \nL = 4 πR 2 F = 4 πR 2 σ [ T 4 -g ( τ S ) T 4 irr ] . (11)  \nThe irradiation temperature as a function of the distance from the star is  \nT irr = ( L /star (1 -A ) 16 πσd 2 ) 1 / 4 (12)  \nwhere L /star is the stellar luminosity (we use L /star = L /circledot ), d is the distance of the planet from the star, and A is the albedo. In order to separate the core effects from the environmental thermal effects, we first take d=1 AU as our standard model. This corresponds to the outer edge of Kepler's detection region. Next, we test closer-in cases of d=0.3 AU and examine the effect on the results. We avoid planets at d < 0.3 AU, from which photoevaporation can become significant, since photoevaporation by the parent star removes (part of) the gaseous envelope (Owen & Wu 2013, 2016), and thus changes the envelope mass in time.  \nTable 2: Model parameters  \nNote. - Set of parameters we use in the model. References: 1 Nettelmann et al. (2011), 2 Morschhauser et al. (2011), 3 Guenther et al. (1992), 4 Sharp & Burrows (2007), 5 Stevenson et al. (1983).  \nThe albedo strongly depends on the atmospheric composition and is an unknown parameter for sub-Neptunes. Here we take all the irradiation to be absorbed by the planet ( A = 0), which is the upper bound of the irradiation effect. Thus, a non-zero albedo planet should be located closer-in for an equivalent irradiation flux.", '3.1. Magma ocean phase': "3.1.1. Role of core energy sources  \nIn Fig. 2 we show the evolution of 5 M ⊕ planets with 90% (4.5 M ⊕ ) core surrounded by 10% (0.5 M ⊕ ) of hydrogen-helium envelope, located at 1 AU from a Sun-like star. The radius (left) and core surface (CEB) temperature (right) are shown in the figure, where the different curves are for different approximations for core energy sources. The CEB solidification occurs when the CEB temperature crosses the dashed horizontal melting line. At this time point (which we call t MO ) the magma ocean phase ends. At first, all cases in the figure are modeled as in the magma ocean phase, i.e., efficient convection, even for t > t MO . Thus, the cooling after the solidification point (thin curves) is overestimated. In Section 3.2 below we show  \nthe effect of the solidification on the results.  \nThe lowest core energy case (blue curve) is achieved for the lower bound of 0.5 Earth radioactive abundance and minimal initial energy from formation ( E acc =0 . 05 E binding ). Then, we increase the initial energy content of the core to 20% of the binding energy ( E acc =0 . 2 E binding ) and include the differentiation energy (red-dotted curve). As is shown in the figure, the long term evolution is the same as the previous case, since this energy is being released efficiently during the magma ocean phase on a timescale of ∼ 10 7 yr. Although the initial core energy content doesn't affect the long term evolution, this initial energy can expand the envelope until it isn't bound to the planet. We test this idea (see appendix A2), and find that for an initial energy of more than 30% of the binding energy, part of the envelope is being lost, for the planets in Fig 2.  \nNext, we increase the contribution from radioactive element decay to levels similar to Earth (dashed-yellow). Since the radioactive energy is being released on Gyr timescale, the long term radius and CEB temperature are higher for this case. In the next case (purple) we added the latent heat by core solidification. The latent heat release by core solidification occurs from inside-out. This is a gradual process during the magma ocean phase, which ends when the CEB solidifies. The latent heat release delays somewhat the CEB solidification by up to 1 Gyr, and slightly increases the radius. This case (purple curve) of a core with Earth abundance of radioactive content, and a latent heat release at solidification is our standard model.  \nFinally, we considered a model with radioactive levels enhanced by a factor 2.5 over Earth. This high radioactive level is motivated by the maximal measured abundance around sun-like stars (Unterborn et al. 2015). This enhanced radioactive heating (dashed-green) in the core is found to be the most important energy contribution to the long term evolution. Variation in the fraction of radioactive elements in the observed range of 0.5-2.5 × Earth abundance, results in up to 5% radius change for planets with solar-metallicity envelope. In general, fraction of radioactive elements and solidification (latent) heat are the most significant core energy sources for the long term. In addition to inflating the radius, the excess heat delays the CEB solidification time (from 4 Gyr to 8 Gyr in this case), and thus prolongs the magma ocean phase. For these cases, many observed sub-Neptunes are likely to be in a magma ocean phase.", '3.1.2. Dependence on the envelope mass': 'In Fig. 3 we follow the change in surface (CEB) temperature during planet evolution, of planets with same core mass surrounded by envelopes of different masses. We use the 4.5 M ⊕ standard core model (purple curve in Fig. 2), and vary the mass of the envelope between 0.1%-20% (0.005-1 M ⊕ ). The horizontal curves in the figure show the point when the CEB solidifies, for the surface pressure and temperature of this planet, as described in appendix A1. In these runs, as before, magma ocean phase conditions are assumed for the entire evolution, thus the cooling after solidification (thin curves) is overestimated.  \nAs is shown in the figure, the magma ocean phase (thick part of each curve) for planets with significant envelopes lasts much longer than in the case of Earth. As the atmospheric mass increases, the temperature on the surface of the core is higher and, as a result, the surface stays molten for longer time. Therefore, for most of the cases in Fig. 3 the duration of the magma ocean phase increases with envelope mass.  \nHowever, with increasing envelope mass, the pressure on the planet surface also increases and so does the melting temperature. For a thick envelope of about 1 M ⊕ (red curve) the melting temperature is above 3000 K and therefore the CEB solidifies earlier than for lower mass envelopes. As a result, the magma ocean duration is limited to up to several Gyrs, since for high mass envelopes ( > 0.5 M ⊕ ) the required higher melting temperatures shorten the magma ocean phase. We find that sub-Neptune planets with envelope masses between 0.01 M ⊕ < M env < 1 M ⊕ make the transition from magma ocean to solidified state during the time we observe them (1-7 Gyr). In the next section we estimate the effect of such a transition on the properties of the planet.', '3.2. Solid state phase': 'In Fig. 4 we present the planetary radius (left) and the CEB temperature (right) for a 4.5 M ⊕ core with 0.5 M ⊕ (10%) and 0.05 M ⊕ (1%) envelopes. The solid curves represent the efficient core cooling as in the magma ocean phase (same as the blue and green curves in Fig. 3). After the CEB solidifies (dashed horizontal temperature line) we calculate for the two scenarios described in Section 2.2: thermal evolution with entirely conductive core (dashed-dotted), and thermal evolution with conductive CEB layer (dashed). The thickness of the conductive layer, according to Eq. 5, is about 100 km, and the timescale for cooling by conduction through this layer is on the order of 10 8 yr.  \nAs is shown in the figure, the CEB temperature rapidly changes at the point where the model (suddenly) assumes a conductive structure. The conductive CEB scenario moderates  \nthe flux from the core but keeps the long term surface temperature similar to the magma ocean case. The conductive boundary layer acts as a bottle-neck for the cooling, but it is not thick enough to slow the cooling substantially. In the conductive core scenario the cooling is much slower; heat is locked in the core while the surface temperature drops.  \nThe effect of solidification on the thermal evolution depends on the envelope mass. In the case of the 0.5 M ⊕ envelope (upper panels) the CEB solidification occurs at a later time than in the 0.05 M ⊕ envelope case. At a later time the core energy budget is smaller (mainly due to radioactive decay). In addition, the cooling at the bottom of the thicker envelope (the core interface) is slower. Therefore, the drop in CEB temperature is more moderate.  \nWhile the effect of the core properties is significant for some of the cases, the effect on the radius is more limited. The maximum radius change for sub-Neptune planets in our sample as resulting from solidification is about 6% percent, or 0.1 R ⊕ . This maximum is achieved for a planet with 0.05 M ⊕ envelope mass (bottom left in Fig. 4). For planets with lower mass envelopes, the contribution of the envelope to the total radius is small. On the other hand, planets with high envelopes masses, such as the 0.5 M ⊕ envelope case (upper left in Fig. 4), have lower core-to-envelope mass ratios and smaller energy contents at solidification, which limit the effect on the radius.', '3.3. Envelope conditions and planet mass': "Since the overall thermal evolution of the core depends also on the thermal evolution of the envelope, the thermal properties of the envelope are expected to change the results. Therefore, we calculate for enhanced (30 × solar) atmospheric opacity, and higher irradiation by the star (d=0.3 AU). In figure 5 we show the radius (left) and CEB temperature (right) in time for identical planets with different atmospheric conditions. The standard case (blue) is compared to a case with enhances opacity (dashed-dotted) and to a case with a stronger irradiation (dashed). Additional representative results appear in Table 3. As is shown in the figure, irradiation has a significant effect on the envelope radius 6 , within the range of previous works (Lopez & Fortney 2014; Howe & Burrows 2015). However, irradiation doesn't change the surface conditions for significant envelope masses ( > 0.01 M ⊕ ). As the outer layers of the envelope expands by irradiation, the conditions at the bottom of the envelope remain similar to the standard case. Thus, the location of the planet (for d > 0.3AU) has a small effect on the core evolution in the presence of a significant envelope (photo-evaporation is  \nnot included in the model).  \nAtmospheric opacity, on the other hand, changes both the radius and the core surface temperature substantially. In the enhanced opacity case the envelope cooling is slower. Thus, the envelope traps the heat from the core and delays the core cooling. As is shown in figure 5, metallicity of 30 × solar keeps the core surface hotter (molten) for much longer time than in the standard case. In the case of lower mass envelope, like for example a 4.5 M ⊕ core with 0.1% envelope, the enhanced atmospheric metallicity (30 × solar) delays the surface solidification from less than 0.5 Gyr to more than 2 Gyr. Moreover, the effect of core thermal properties on the radius become more important for the high atmospheric opacity. Uncertainty in radioactive elements, in the range of 0.5-2.5 Earth-ratio, results in up to 15% radius change for planets with enhanced envelope metallicity, in comparison to only 5% for similar planets with solar-metallicity.  \nWe also vary the mass of the core by a factor of two (2.25 M ⊕ and 9 M ⊕ ) for different envelope masses (examples in Table 3). We find that cores of different masses with the same envelope mass result in similar cooling rates of the core surface 7 , i.e., the envelope mass, and not the envelope fraction, is the key parameter driving the core thermal evolution. The reason is that during the magma ocean phase the pressure-temperature conditions at the bottom of the envelope are determined by the adiabatic structure of the envelope. Hence, the same envelope mass leads to similar pressure-temperature conditions at the core surface. In figure 6 we show this trend for two planets with different total masses, but the same envelope mass: the difference in magma ocean duration (CEB temperature) between a 10 M ⊕ planet with 5% envelope and a 5 M ⊕ planet with 10% envelope (same envelope mass) is smaller than the uncertainty within the core thermal properties for each planet. While the two cores feature a similar thermal evolution until solidification, the solidification point is somewhat different. The 10 M ⊕ planet reaches the CEB solidification earlier than the 5 M ⊕ planet, because of higher surface pressure which increases the melting (solidification) temperature. The radii of the 5 M ⊕ planet is larger, because of irradiation effect on cores with lower gravity (e.g., Lopez & Fortney 2014). Conversely, the higher (lower) solidification temperature for more (less) massive cores, due to surface pressures, and the sensitivity to thermal effects by the higher (lower) gravity, are found to have only small effects on the radius variation by the core properties. Also here, high atmospheric opacity intensifies the core thermal effects.", '4. Discussion': "Our results indicate that many of the observed sub-Neptune planets are in the magma ocean phase for several Gyrs. This is very different than the short magma ocean phase of planets in our solar system (e.g., Elkins-Tanton 2012). Because of its prolonged existence, it is likely that core and envelope interact strongly; not only regarding their thermal properties (as modeled here) but also regarding their composition. For example, solubility of the envelope hydrogen in the silicate melt (Hirschmann et al. 2012) can be significant for high temperature and pressure conditions as in sub-Neptunes (Chachan & Stevenson 2018). Consequently, sub-Neptune planets with less than 1 M ⊕ envelope may contain significant volatiles fraction in their cores, which affects the rock thermal and physical properties.  \nConversely, the long interaction of the molten convective core with the hydrogen envelope may also enrich the envelope with metals from the molten core. Currently, there is a lack of knowledge (experimental data as well as modelling) of rock-envelope interaction for sub-Neptune conditions. Most of the current knowledge pertains to processes on Earth-like terrestrial planets (e.g., Abe & Matsui 1986; Hirschmann 2012). Studies of rock-envelope interaction for sub-Neptune pressure-temperature conditions are necessary to improve our understanding of the structure and thermal properties of these objects.  \nAs we show in this work, the envelope mass determines the state and the cooling rate of the core. The Kepler's data reveals a statistical dip in planetary radius between two picks of 1.3 R ⊕ and 2.6 R ⊕ (Fulton et al. 2017). This valley, which divides the close-in planets into two populations, can be a result of envelope mass loss by photoevaporation (Owen & Wu 2017; Jin & Mordasini 2018). Envelope mass loss during the evolution shorten the magma ocean phase (as long as t evap < t MO ) and thus affects the thermal evolution of the core. Thus, the valley actually divides the planets also into two thermal populations, where the radius pick of 1.3 R ⊕ is for bare (solid state) cores, and the 2.6 R ⊕ pick is for magma ocean phase cores (M env > 0.05 M ⊕ ). Future studies of magma-envelope chemical interaction can provide atmospheric markers to be detected, in order to distinguish between the two populations.  \nThis study is focused on sub-Neptune planets around sun-like star (G-type). SubNeptune planets around M-type stars, which appear to be common (e.g., Mulders et al. 2015), may have different formation environment (e.g., Kennedy & Kenyon 2008; Ormel et al. 2017), and thus different core thermal contribution. For example, different core formation timescale in a lower mass disk will change the formation energy left in the core. Moreover, low mass disk may have different composition. As a result the differentiation energy, which is derived from the iron-to-rock ratio will change. The radioactive element abundances and the rock latent heat, which we find to contribute the most for the planet long term evolution, are derived from the rock mineralogy, which depends on metals abundances and thermal  \nconditions of the building blocks in the disk. As a result, the (Earth-like) values for core energy sources we used in this work (Secion 2.3) will change for M-type stars.  \nThis work consider rocky cores, without any fraction of ice. Ice-rich cores have substantial effects on both radius-mass relation and the atmospheric properties (Chen & Rogers 2016). However, during the solidification of the molten core water would be expelled from the core (e.g., Elkins-Tanton 2008). For the planetary and envelope mass range of sub-Neptunes, water is found to stay in the envelope in vapor form; the surface pressure-temperature conditions don't allow for liquid water for the parameter range we studied here. As we show in figure 3, envelope mass larger than 0.02 M ⊕ keeps the surface temperature above 1000 K for more than 10 Gyr. Thus, if the core contains some fraction of ice, the envelope is presumably saturated with vapor and thus denser than the H, He envelope in our model. In this case, the increase in core radius by the icy (lighter) materials diminishes by the decrease in envelope radius by the higher mean molecular weight. Moreover, ice-rich cores are expected to contain less radioactive elements, which are scaled with the rocky part. Therefore, the core energy source of radioactive heating is reduced.  \nClose-in planets usually have low envelope masses, and thus exhibits limited radius change by core thermal contribution, as is discussed in Section 3.2. However, for close-in planets additional mechanisms, that are not included in our model, can affect the contribution of the core to the overall thermal evolution: (1) photoevaporation by the parent star. The early released core energy (formation and differentiation) extends the envelope radius and hence accelerates photoevaporation. Photoevaporation is expected to reduce the envelope mass, and therefore shorten the magma ocean phase duration. Including photoevaporation simulation (e.g., Murray-Clay et al. 2009) in the core-envelope evolution model is essential for modeling of close-in planets, and we would like to address it in a future work. (2) close-in planets experience tidal forces by the parent star. Tidal heating is a continuous energy source in the core, and may be significant for super-Earth cores 8 (Efroimsky 2012). (3) late planeteseimal capture can add energy during the planet evolution, as well as increase the atmospheric metallicity. This contribution is relevant in particular for planets less massive than 10 M ⊕ and with envelope mass fractions less than 10% (Chatterjee & Chen 2018).  \nThe core model appears in this work is a first order estimate for core thermal evolution. In detailed (envelope-less) geophysical models (e.g., van den Berg et al. 2019), the transfer between the different core thermal phases is continuous, and the heat transport depends on thermal and physical properties of the core minerals. As we show here the core-envelope  \nthermal effects are mutual, i.e., the envelope changes the core thermal state in time. Thus, in order to improve the existing sub-Neptune models one should link geophysical core model to envelope thermal evolution. Such a combined model is challenging, because of its selfdependent nature, but is necessary in order to better understand interiors of sub-Neptune planets.", '5. Conclusion': 'We have modeled the thermal evolution of sub-Neptune planets with core and hydrogenhelium envelope on one structure grid. Our model divides the evolution of the core into three phases: initial (formation phase), efficient cooling (magma ocean phase), and inefficient cooling (solid state phase). We have examined the contribution of the core energy sources to the thermal evolution of the planet as a whole. In particular, we followed the mutual core-envelope thermal effects on the core solidification and the planet radius evolution.  \nWe summarize our main conclusions below:  \n- 1. Most of the observed sub-Neptune planets are in the magma ocean phase (molten surface). We find that the duration of the magma ocean phase for planets with envelope masses between 0.01 M ⊕ - 1 M ⊕ lies between 1 and 7 Gyrs.\n- 2. Because of its efficient cooling, the magma ocean phase renders the evolution insensitive to the initial conditions. For this reason, the initial thermodynamic state of the core (heat of formation and iron differentiation) does not influence the radius evolution for more than several 10 7 years.\n- 3. Radioactive decay is the most significant energy source to affect the planet radius, and the latent heat from solidification is the second. In the long term, the planet radius variation as a result of uncertainty in these core energy sources is at most 15%.\n- 4. After solidification of the CEB the heat flux from the core decreases further. We calculate that the variation in radius due to uncertainties regarding the post-solidification phase is no larger than 6%.\n- 5. Overall, for typical model parameters, the contribution from the thermal state of the core to the planet radius is rather limited (a few percent at most). Therefore, the inferred envelope mass from mass-radius relation is mostly proportional to the envelope (H/He) mass fraction.  \n- 6. The atmospheric opacity significantly prolongs the magma ocean phase, and amplifies the core effects on the envelope evolution. Irradiation (for d > 0.3AU), on the other hand, has only a minor effect on the core evolution.  \nWe thank Dave Stevenson, Wim van Westrenen, and Arie van den Berg for useful comments and discussions. A.V. thanks Attay Kovetz for helpful ideas for code modifications. A.V. acknowledges support by the Amsterdam Academic Alliance (AAA) Fellowship, and by the Netherlands Origins Center. C.W.O. is supported by the Netherlands Organization for Scientific Research (NWO; VIDI project 639.042.422)', 'A1. Rock melting curve': "For the calculation of rock melting curve we use solidus and liquidus pressure-temperature relation as followed: for pressures below 10 GPa we take melting curves for peridotite, based on de Smet et al. (1999), and consistent with melting profiles from Hirschmann (2000) (see also Noack et al. 2017):  \nT solid = 1409 . 15 K +134 . 2 p -6 . 581 p 2 +0 . 1054 p 3  \nT liquid = 2035 . 15 K +57 . 46 p -3 . 487 p 2 +0 . 0769 p 3  \nFor pressures above 10 GPa we use melting curves for perovskite, based on Stamenkovi'c et al. (2011) for the solidus. We use a fixed difference between solidus and liquidus for lack of experimental data:  \nT solid = 1835 K +36 . 918 p -0 . 065444 p 2 +0 . 000076686 p 3 -3 . 09272 × 10 -8 p 4  \nT liquid = 1980 K +36 . 918 p -0 . 065444 p 2 +0 . 000076686 p 3 -3 . 09272 × 10 -8 p 4  \nWe define our melting curve at 40% between the solidus and liquidus. We assume that up to 40% melt the rock behaves like a solid, and above that as a magma ocean. Thus, the critical melting temperature that we use as a transition between a rock behaving as a solid and a rock behaving as a liquid is:  \nT melt = 1659 . 55 K +103 . 504 p -5 . 3434 p 2 +0 . 094 p 3  \nfor p < 10 GPa, and  \nT melt = 1893 K +36 . 918 p -0 . 065444 p 2 +0 . 000076686 p 3 -3 . 09272 × 10 -8 p 4  \nfor p > 10 GPa. In figure 7 we show the resulting rock melting temperature as a function of the CEB pressure for the range of sub-Neptune planets.", 'A2. Envelope mass loss': "As we show in Section 3, the initial core energy content doesn't affect the long term radius evolution. However, this energy can bloat the envelope until it isn't bound to the planet any more (Ginzburg et al. 2016, 2017). Here, we test this idea by adding high fractions of the formation core binding energy to the early phase core, and following the radius of the planet in comparison to its Hill radius. The expansion of the envelope is an outcome of the self-consistent core and envelope thermal evolution, accounting for the core heat flux and cooling, material properties (e.g., tabular EOSs and opacity) and their time dependency. We find that part of the envelope is being lost (R p > R Hill ) when high fractions of formation energy is initially stored in the core. As the envelope expands from the core luminosity, the inner part of the envelope remains gravitationally bound to the core, and later on it cools more efficiently (than a thick envelope) and contracts. What fraction of the envelope is removed depends mainly on the initial core energy content, atmospheric opacity profile, and the distance from the star (Hill sphere). Under the conditions of our model, we find that envelope loss by the core energy starts when more than 30% of the accretion energy remains in the core after its formation. For core energy content of 40%, for example, our standard model of 4.5 M ⊕ core cannot retain more than 0.1 M ⊕ envelope. If less than 20% of the core binding energy is left in the core after its formation, all the planetary envelopes within the range of this work survive.", 'REFERENCES': "Abe, Y. & Matsui, T. 1986, J. Geophys. Res., 91, E291  \nAnders, E. & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197  \nBaraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315  \nBatalha, N. M., Rowe, J. F., Bryson, S. T., Barclay, T., Burke, C. J., Caldwell, D. A., Christiansen, J. L., Mullally, F., Thompson, S. E., Brown, T. M., Dupree, A. K., Fabrycky, D. C., Ford, E. B., Fortney, J. J., Gilliland, R. L., Isaacson, H., Latham, D. W., Marcy, G. W., Quinn, S. N., Ragozzine, D., Shporer, A., Borucki, W. J., Ciardi, D. R., Gautier, III, T. N., Haas, M. R., Jenkins, J. M., Koch, D. G., Lissauer, J. J., Rapin, W., Basri, G. S., Boss, A. P., Buchhave, L. A., Carter, J. A., Charbonneau, D., Christensen-Dalsgaard, J., Clarke, B. D., Cochran, W. D., Demory, B.-O., Desert, J.-M., Devore, E., Doyle, L. R., Esquerdo, G. A., Everett, M., Fressin, F.,  \nGeary, J. C., Girouard, F. R., Gould, A., Hall, J. R., Holman, M. J., Howard, A. W., Howell, S. B., Ibrahim, K. A., Kinemuchi, K., Kjeldsen, H., Klaus, T. C., Li, J., Lucas, P. W., Meibom, S., Morris, R. L., Prˇsa, A., Quintana, E., Sanderfer, D. T., Sasselov, D., Seader, S. E., Smith, J. C., Steffen, J. H., Still, M., Stumpe, M. C., Tarter, J. C., Tenenbaum, P., Torres, G., Twicken, J. D., Uddin, K., Van Cleve, J., Walkowicz, L., & Welsh, W. F. 2013, ApJS, 204, 24  \nBean, J. L., D'esert, J.-M., Kabath, P., Stalder, B., Seager, S., Miller-Ricci Kempton, E., Berta, Z. K., Homeier, D., Walsh, S., & Seifahrt, A. 2011, ApJ, 743, 92  \nChachan, Y. & Stevenson, D. J. 2018, ApJ, 854, 21  \nChatterjee, S. & Chen, H. 2018, ApJ, 852, 58  \nChen, H. & Rogers, L. A. 2016, ApJ, 831, 180  \nCoughlin, J. L., Mullally, F., Thompson, S. E., Rowe, J. F., Burke, C. J., Latham, D. W., Batalha, N. M., Ofir, A., Quarles, B. L., Henze, C. E., Wolfgang, A., Caldwell, D. A., Bryson, S. T., Shporer, A., Catanzarite, J., Akeson, R., Barclay, T., Borucki, W. J., Boyajian, T. S., Campbell, J. R., Christiansen, J. L., Girouard, F. R., Haas, M. R., Howell, S. B., Huber, D., Jenkins, J. M., Li, J., Patil-Sabale, A., Quintana, E. V., Ramirez, S., Seader, S., Smith, J. C., Tenenbaum, P., Twicken, J. D., & Zamudio, K. A. 2016, ApJS, 224, 12  \nde Smet, J. H., van den Berg, A. P., & Vlaar, N. J. 1999, Lithos, 48, 153  \nD'esert, J.-M., Bean, J., Miller-Ricci Kempton, E., Berta, Z. K., Charbonneau, D., Irwin, J., Fortney, J., Burke, C. J., & Nutzman, P. 2011, ApJ, 731, L40  \nEfroimsky, M. 2012, ApJ, 746, 150  \nElkins-Tanton, L. T. 2008, Earth and Planetary Science Letters, 271, 181  \n-. 2012, Annual Review of Earth and Planetary Sciences, 40, 113  \nFortney, J. J., Mordasini, C., Nettelmann, N., Kempton, E. M.-R., Greene, T. P., & Zahnle, K. 2013, ApJ, 775, 80  \nFulton, B. J., Petigura, E. A., Howard, A. W., Isaacson, H., Marcy, G. W., Cargile, P. A., Hebb, L., Weiss, L. M., Johnson, J. A., Morton, T. D., Sinukoff, E., Crossfield, I. J. M., & Hirsch, L. A. 2017, AJ, 154, 109  \nGinzburg, S., Schlichting, H. E., & Sari, R. 2016, ApJ, 825, 29  \n-. 2017, ArXiv e-prints  \nGuenther, D. B., Demarque, P., Kim, Y.-C., & Pinsonneault, M. H. 1992, ApJ, 387, 372  \nGuillot, T. 2010, A&A, 520, A27  \nGuillot, T., Chabrier, G., Gautier, D., & Morel, P. 1995, ApJ, 450, 463  \nHirschmann, M. M. 2000, Geochemistry, Geophysics, Geosystems, 1, 1042  \n-. 2012, Earth and Planetary Science Letters, 341, 48  \nHirschmann, M. M., Withers, A. C., Ardia, P., & Foley, N. T. 2012, Earth and Planetary Science Letters, 345, 38  \nHoward, A. W., Marcy, G. W., Bryson, S. T., Jenkins, J. M., Rowe, J. F., Batalha, N. M., Borucki, W. J., Koch, D. G., Dunham, E. W., Gautier, III, T. N., Van Cleve, J., Cochran, W. D., Latham, D. W., Lissauer, J. J., Torres, G., Brown, T. M., Gilliland, R. L., Buchhave, L. A., Caldwell, D. A., Christensen-Dalsgaard, J., Ciardi, D., Fressin, F., Haas, M. R., Howell, S. B., Kjeldsen, H., Seager, S., Rogers, L., Sasselov, D. D., Steffen, J. H., Basri, G. S., Charbonneau, D., Christiansen, J., Clarke, B., Dupree, A., Fabrycky, D. C., Fischer, D. A., Ford, E. B., Fortney, J. J., Tarter, J., Girouard, F. R., Holman, M. J., Johnson, J. A., Klaus, T. C., Machalek, P., Moorhead, A. V., Morehead, R. C., Ragozzine, D., Tenenbaum, P., Twicken, J. D., Quinn, S. N., Isaacson, H., Shporer, A., Lucas, P. W., Walkowicz, L. M., Welsh, W. F., Boss, A., Devore, E., Gould, A., Smith, J. C., Morris, R. L., Prsa, A., Morton, T. D., Still, M., Thompson, S. E., Mullally, F., Endl, M., & MacQueen, P. J. 2012, ApJS, 201, 15  \nHowe, A. R. & Burrows, A. 2015, ApJ, 808, 150  \nIkoma, M. & Hori, Y. 2012, ApJ, 753, 66  \nJin, S. & Mordasini, C. 2018, ApJ, 853, 163  \nJontof-Hutter, D., Ford, E. B., Rowe, J. F., Lissauer, J. J., Fabrycky, D. C., Van Laerhoven, C., Agol, E., Deck, K. M., Holczer, T., & Mazeh, T. 2016, ApJ, 820, 39  \nKarato, S.-I. & Wu, P. 1993, Science, 260, 771  \nKennedy, G. M. & Kenyon, S. J. 2008, ApJ, 682, 1264  \nKovetz, A., Prialnik, D., & Shara, M. M. 1988, ApJ, 325, 828  \nKovetz, A., Yaron, O., & Prialnik, D. 2009, MNRAS, 395, 1857  \nLopez, E. D. & Fortney, J. J. 2014, ApJ, 792, 1  \nLopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59  \nMordasini, C., Alibert, Y., Georgy, C., Dittkrist, K.-M., Klahr, H., & Henning, T. 2012, A&A, 547, A112  \nMorley, C. V., Fortney, J. J., Kempton, E. M.-R., Marley, M. S., Visscher, C., & Zahnle, K. 2013, ApJ, 775, 33  \nMorschhauser, A., Grott, M., & Breuer, D. 2011, Icarus, 212, 541  \nMulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 814, 130  \nMurray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23  \nNettelmann, N., Fortney, J. J., Kramm, U., & Redmer, R. 2011, ApJ, 733, 2  \nNoack, L., Honing, D., Rivoldini, A., Heistracher, C., Zimov, N., Journaux, B., Lammer, H., Van Hoolst, T., & Bredehoft, J. H. 2016, Icarus, 277, 215  \nNoack, L., Rivoldini, A., & Van Hoolst, T. 2017, Physics of the Earth and Planetary Interiors, 269, 40  \nOrmel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1  \nOwen, J. E. & Wu, Y. 2013, ApJ, 775, 105  \n- -. 2016, ApJ, 817, 107\n- -. 2017, ApJ, 847, 29  \nPoirier, J.-P. 2000, Introduction to the Physics of the Earth's Interior, 326  \nPollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., & Greenzweig, Y. 1996, Icarus, 124, 62  \nRogers, L. A. & Seager, S. 2010, ApJ, 712, 974  \nSaumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713  \nSharp, C. M. & Burrows, A. 2007, ApJS, 168, 140  \nSolomon, S. C. 1979, Physics of the Earth and Planetary Interiors, 19, 168  \nStamenkovi'c, V., Breuer, D., & Spohn, T. 2011, Icarus, 216, 572  \n- Stamenkovi'c, V., Noack, L., Breuer, D., & Spohn, T. 2012, ApJ, 748, 41 Stevenson, D. J. Fluid dynamics of core formation., ed. H. E. Newsom & J. H. Jones, 231-249 Stevenson, D. J., Spohn, T., & Schubert, G. 1983, Icarus, 54, 466 Tackley, P. J., Ammann, M., Brodholt, J. P., Dobson, D. P., & Valencia, D. 2013, Icarus, 225, 50 Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64 Turcotte, D. L. & Oxburgh, E. R. 1967, Journal of Fluid Mechanics, 28, 29 Unterborn, C. T., Johnson, J. A., & Panero, W. R. 2015, ApJ, 806, 139 Valencia, D., O'Connell, R. J., & Sasselov, D. D. 2007, ApJ, 670, L45 van den Berg, A. P., Yuen, D. A., Beebe, G. L., & Christiansen, M. D. 2010, Physics of the Earth and Planetary Interiors, 178, 136 van den Berg, A. P., Yuen, D. A., Umemoto, K., Jacobs, M. H. G., & Wentzcovitch, R. M. 2019, Icarus, 317, 412 Vazan, A., Helled, R., Kovetz, A., & Podolak, M. 2015, ApJ, 803, 32 Vazan, A., Kovetz, A., Podolak, M., & Helled, R. 2013, MNRAS, 434, 3283 Vazan, A., Ormel, C. W., & Dominik, C. 2018, A&A, 610, L1  \nWolfgang, A. & Lopez, E. 2015, ApJ, 806, 183  \nFig. 2.- Radius (left) and CEB temperature (right) of 4.5 M ⊕ cores with 0.5 M ⊕ (10%) hydrogen-helium envelope. The different curves are for different assumptions for core energy sources: minimal radioactive heating and low formation energy ( E acc =0 . 05 E binding ) (blue), 4 times higher accretion energy ( E acc =0 . 2 E binding ) and including differentiation energy (reddashed), Earth ratio of radioactive heating (yellow-dashed), including latent heat (blue), and maximal radioactive heating together with all previous energy sources (green-dashed). The horizontal line indicates the core surface solidification.  \n<!-- image -->  \nFig. 3.- CEB temperature of 4.5 M ⊕ cores with different envelope masses. The percentages denote the mass fractions of the envelope with respect to the total planet mass. The horizontal lines indicate the (pressure-dependent) CEB solidification temperature of each planet (see appendix A1). The evolution is modelled by efficient magma ocean cooling, which overestimates the cooling after solidification (thin curves).  \n<!-- image -->  \nFig. 4.- Radius (left) and CEB temperature (right) for a 4.5 M ⊕ core with an envelope. Top panels: 10% envelope mass (0.5 M ⊕ ); bottom panels: 0.1% envelope mass (0.005 M ⊕ ). The solid curve is for efficient core cooling (magma ocean). The evolution with a conductive CEB layer (dashed) and with a conductive core (dashed-dotted) are shown after the CEB reaches the solidification temperature (horizontal dashed).  \n<!-- image -->  \nFig. 5.- Radius (left) and CEB temperature (right) for 5 M ⊕ planets with 0.5 M ⊕ (10%) envelope. The different curve styles are for different envelope conditions: standard case of solar metallicity opacity located at 1 AU (solid blue), planet with enhanced atmospheric opacity of 30 × solar (dashed-dotted black), and planet with higher irradiation as located at 0.3 AU (dashed red). The different curves of each color represent the different core cooling scenarios after core solidification (phase 3).  \n<!-- image -->  \nFig. 6.- Radius (left) and CEB temperature (right) for 10 M ⊕ (dashed black) and 5 M ⊕ (solid blue) planets with the same envelope mass of 0.5 M ⊕ . The different curves of each planet represent the different core cooling scenarios after solidification (phase 3). The cores feature a similar thermal evolution until solidification.  \n<!-- image -->  \nTable 3: Parameter study  \nFig. 7.- Rock melting curve we use in this work to determine the rocky core solidification. We define our melting curve (solid black) at 40% between the solidus (red-dashed) and liquidus (blue-dashed). See appendix A1 for details.  \n<!-- image -->"}