© The Authors 2025

1 Introduction

Intermediate-mass planets make up the majority of the detected exoplanet population. Their interiors are assumed to be composed of mostly rocks and ices surrounded by a hydrogen–helium (H–He) envelope. Unlike super-Earth planets, their radii are inflated by the presence of light elements (e.g. Rogers et al. 2011; Lozovsky et al. 2018). The mass fraction of H-He strongly influences the radius of these planets regardless of the internal energy or the external stellar flux. Apart from the existence of an H-He atmosphere and a deep interior that is dominated by heavy elements, their detailed structure is highly uncertain. Constrained by only radius and mass, the composition profile of the detected exoplanet population is highly degenerate. Even for Uranus and Neptune, where gravitational data provide additional constraints on the planetary density distribution, the structure and bulk composition is still uncertain (e.g. Nettelmann et al. 2013; Neuenschwander & Helled 2022; Morf et al. 2024).

Since the internal planetary structure can change as a planet evolves, static models may be further constrained by modeling planetary evolution. Often, planets are assumed to be fully convective, i.e., adiabatic, which imposes strong assumptions on both their evolution and structure. However, it recently became evident that assuming a homogeneous and adiabatic structure is an oversimplification. Composition gradients, which inhibit convection, can form due to the continuous accretion of solids and gas (e.g. Valletta & Helled 2022) or by accreting different materials when migrating over ice lines in the protoplanetary disk (e.g. Schneider & Bitsch 2021; Eberlein et al. 2024). Indeed, it has been suggested that both Uranus and Neptune have composition gradients in their deep interiors that can survive for billions of years (Vazan & Helled 2020; Tejada Arevalo 2025). Composition gradients are not always primordial and can be caused by other processes, such as phase separation and the rainout of heavy elements. This further complicates the picture and can create (additional) composition gradients during the evolution of a planet (e.g., Stevenson 1980; Lorenzen et al. 2009; Redmer et al. 2009; Püstow et al. 2016; Schöttler & Redmer 2018; Chang et al. 2024; Cano Amoros et al. 2024). In addition, depending on the opacity, a deep radiative layer can also form within the planet (e.g., Howard et al. 2023; Müller & Helled 2024).

In non-convective regions, the heat transport is governed by multiple thermal transport mechanisms – radiative, diffusive, and conductive. Usually, the diffusive transport where hot particles diffuse into regions with colder temperatures is insignificant compared to radiation and conduction. In the outer parts of the planet, it is mostly the photons that contribute to thermal transport. Material emits and reabsorbs thermal radiation, and given a temperature gradient, this leads to a net diffusion of thermal energy that depends on the Rosseland-mean opacity (Freedman et al. 2014; Marigo & Aringer 2009; Marigo et al. 2024). In denser regions, where the temperature is still sufficiently low so that the material is not ionized, conductivity is governed by the vibration of molecules, similar to phonons in solids. Ross et al. (1984) reviewed the phonon model for metals applied to anharmonic non-metallic crystals. This model was used for the thermal conductivity of the Earth’s core-mantle boundary (Stamenković et al. 2011). More recently, French (2019) used ab initio calculations of water to create an empirical fit of the thermal conductivity contribution from the nuclei. In the deep hot interior of the planet, free electrons enhance the thermal conductivity. Cassisi et al. (2007) elaborated on the chemical picture of the electron conductivity by considering the density and scattering processes of electrons. French & Redmer (2017) complemented the conductivities of the water nucleus by calculating the electron contribution.

In this paper, we investigate the influence of thermal conductivity on the evolution of Neptunes and sub-Neptunes. We compare different cases for the conductivity to assess how it affects the radius evolution and the inferred internal structure. Our paper is structured as follows. In Section 2, we summarize the physical processes and the methods. In Section 3 we show the results of our study and discuss them in Section 4. Finally, our summary and conclusions are presented in Section 5.

2 Methods

We modeled the planetary evolution by solving the stellar structure equations (e.g. Kippenhahn et al. 2013) using the Henyey method (Henyey et al. 1965) with the Modules for Experiments in Stellar Astrophysics (MESA) code (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). We modified the code to make it applicable to medium-sized planets. We used the equation of state (EoS) presented by Müller et al. (2020b,a); Müller & Helled (2021, 2024) that combines the EoS with non-ideal interaction (Chabrier & Debras 2021) for H–He with the heavy-element EoS for H₂O and SiO₂ (More et al. 1988; Vazan et al. 2013) in a 50/50 rock-water mixture. The tables were implemented using a modified version of the MESA extension ‘custom EoS’ (Knierim & Helled 2024; Helled et al. 2025). The modifications enabled the solver to utilize linear interpolations within the EoS tables rather than the standard bicubic spline interpolation scheme employed by MESA.

Furthermore, we changed the radiative opacities and thermal conductivities as discussed below. We started the evolution with a given initial thermal state, composition profile, and mass as discussed in Section 2.2 and evolved the planet for 10 Gyr. The atmospheric cooling was simulated using the irradiated gray implementation in MESA based on Guillot & Havel (2011). We adopted an equilibrium temperature of T_eq = 400 K and a mean visible-to-thermal opacity ratio of κ_V/κ_th = 0.03 based on the fit provided by Poser & Redmer (2024). Although our chosen equilibrium temperature lies outside the original fit range of 500–4000 K, the authors estimate the fit uncertainty to be approximately ±0.1, which exceeds the difference in the opacity ratio between 400 K and 500 K. We used the Ledoux criterion to determine if a region is convective or non-convective. Because of numerically challenging simulations, we did not change the composition over time.

2.1 Thermal conductivity and opacity

A temperature gradient, ∇T, within a planet creates a thermal flux, F_thermal. The flux within the planetary interior can be modeled as a diffusive process such that it follows Fourier’s law: $$F_{\text{thermal}}=-k\nabla T,$$(1)

given the thermal conductivity, k. The total conductivity, k_tot, which is calculated as $$k_{\text{tot}}=k_{\text{rad}}+k_{\text{vib}}+k_{\text{elec}},$$(2)

is the sum of all processes, where κ_rad is the conductivity contribution from photons, while k_vib and k_elec are the conductivity contributions from vibrations in a dense fluid and from electrons, respectively.

The radiative conductivity can be converted to radiative opacity, κ_rad, which measures how much light is absorbed in a medium using $${\kappa }_{\text{rad}}=\frac{16\sigma }{3}\frac{T^3}{\rho k_{\text{rad}}},$$(3)

where σ is the Stefan-Boltzmann constant, T the temperature, and ρ the density. In the following, we use the term radiative conductivity to refer to thermal transport by photons. The standard low temperature opacity tables in MESA that cover high metallicities are an extended version of the Freedman opacities (Freedman et al. 2008, 2014). They have been tabulated for a fixed solar H-He ratio and different metallicity fractions up to Z = 1, relevant for planets. However, the reliable range stops at log(T/K) = 3.8 and significantly underestimates the opacity. In this region, MESA uses high temperature opacities instead, which do not cover high Z ≥ 0.2 materials. To improve on that, we created tables similar to the Freedman tables using the AESOPUS2.1 (Marigo & Aringer 2009; Marigo et al. 2024) web interface. The recently updated method allowed us to the create tables in the range for temperatures up to log(T/K) = 4.5 with arbitrary values for the metallicity. The AESOPUS2.1 opacities were compared to the Freedman tables and to another set of opacities from Malygin et al. (2014). In the regime relevant for planets (high log R), the AESOPUS2.1 opacities are generally higher than the opacities from the other two studies, especially below temperatures of log(T/K) ≲ 3.3 for Malygin et al. (2014) and log(T/K) ≲ 3.5 for the Freedman tables. The differences might originate from a more complete list of molecular lines over the entire temperature region (Marigo et al. 2024). As a result, the AESOPUS2.1 opacities seem to be more reliable. We note that the latest release of MESA includes a set of opacity tables made with the AESOPUS2.1 code. However, these tables only include low metallicities and are therefore inappropriate for modeling a large range of planetary types. (For details on how we created the tables and used them in MESA, see the Appendix A.)

For the vibrational conductivity, we used two models and compared them with each other. From the theory of phonons in metals, Ross et al. (1984) reviewed the application to anhar-monic non-metallic crystals. We used the formulation from Stamenković et al. (2011) and took the reference conductivity for MgSiO₃ k_{vib, ref}(ρ_ref = 3.87 g/cm³, T_ref = 2000 K) = 0.71 W/m/K (see Table 1 in Stamenković et al. 2011). More details can be found in Appendix B. The second model we used is based on the empirical fit to ab initio calculated conductivities of water (Eq. (7) in French 2019). They provide an analytical function that is suitable for densities in the range of ρ ∈ (0.2, 10) g/cm³ and temperatures of T ∈ (600, 50 000) Κ.

The electron contribution to the thermal conductivity in MESA is modeled by default with improved tables from Cassisi et al. (2007). This conductivity assumes fully ionized matter, which is not suitable for planets. We compared this conductivity with the fit to the ab initio calculations for the electron conductivity of partially ionized water (Eq. (30) in French & Redmer 2017). Their fit is valid in the same temperature density region as the work in French (2019).

Because MESA uses the opacity formulation of the thermal transport, we combined the contributions from vibrations and electrons as a sum of the conductivities and converted the result into an opacity. This opacity was then combined with the radiative opacity using the reciprocal sum to get the total opacity. To avoid truncation errors in the reciprocal sum due to very low conductivities, we imposed a lower bound of k ≥ 10⁻⁶ W/m/K on all the conductivity values. This threshold has little effect on the planetary evolution since it is at least five orders of magnitude smaller than the lowest conductivity obtained in any of our models. To verify this, we performed a simulation with a reduced lower limit of k ≥ 10⁻⁸ W/m/K and found a deviation of only ~10⁻⁹ R_⊕ in radius after 10 Gyr.

2.2 Model setups

Our models correspond to planets shortly after their formation process has been completed. The initial model was constructed given a mass, composition profile, and primordial entropy. We simulated the evolution for three different planetary masses of M_p = 5, 10, and 15 Μ_⊕. To account for the uncertainty in the primordial planetary entropy, we considered three different central specific entropies of s_{center, i} = 0.5, 0.6, and 0.7 k_B/m_H, which we refer to as “cold,” “warm,” and “hot,” respectively. We created models with a wide and narrow composition gradient. Except for what we report in Section 3.2 and if not otherwise stated, we used the wide transition profile. The envelope metallicity was fixed to Z_env = 0.20 and the deep interior has Z_int = 1.0. The transition region has a width of Δq = 0.1 in terms of the normalized mass coordinate q for the wide gradient and Δq = 0.01 for the narrow gradient. The position was fixed to achieve a bulk composition of Z_bulk = 0.95 with a hydrogen-helium ratio of Y/X = 0.34 based on the solar photosphere (Asplund et al. 2009). Figure 1 shows the initial specific entropy, composition, and temperature profile for the wide composition gradient (see Figure C.1 for the narrow composition gradient). Details on how we used MESA’s relax options to create our initial model are presented in Appendix C.

To compare the different approaches to the conductivity, we created multiple model cases: Cond-1 uses the vibrational and electronic conductivity for water under high pressure (French & Redmer 2017; French 2019), Cond-2 uses the default electron conductivity of MESA for fully ionized material based on Cassisi et al. (2007), Cond-3 assumes a constant electron plus vibration conductivity based on the conductivity at the Earth’s core-mantle boundary of k_elec = 4 W/m/K (e.g. Stevenson et al. 1983; Lobanov et al. 2021), and Cond-4 uses the phonon model for the vibrational conductivity (Stamenković et al. 2011) and the electron conductivity for partially ionized water (French & Redmer 2017). All the cases use the AESOUPS tables for the radiative conductivity. We summarize the models in Table 1.

Overall, we recommend using Cond-1 because it is the most updated model. However, the conductivity of Cond-1 is currently only available for H₂O, and we hope that future work investigates the thermal conductivity of various materials under a large range of pressures and temperatures.

Fig. 1 Initial profiles for the wide composition gradient, showing specific entropy, temperature, and composition as functions of normalized mass. The colors blue, orange, and yellow correspond to different primordial entropies. The dotted, solid, and dashed lines correspond to planets with a mass of 5 Μ⊕, 10 Μ⊕, and 15 Μ⊕, respectively. Top: Specific entropy vs. normalized mass of the initial model for the wide composition profile. The metallicity of the composition profile is shown by a black line, with its y-axis on the right-hand side. Bottom: Temperature vs. normalized mass of the initial model for the wide composition profile. The black line corresponds to the metallicity.

3 Results

3.1 The importance of thermal conductivity

We first compared the different conductivity cases: Cond-1, Cond-2, and Cond-3. We compared Cond-4 with Cond-1 separately because they only differ in a small region within the planet. We ran the evolution for 10 Gyr and compared the three cases. Figure 2 shows the conductivity as a function of the planetary radius at t = 5 Gyr. We chose this age to better reflect the typical ages of observed exoplanets, in contrast to the final simulation time of t = 10 Gyr. The qualitative differences between the conductivities remained similar throughout the evolution.

In the outer low-density parts of the planet, the radiative conductivity dominates. In the inner parts, where material is at least partially ionized, most of the thermal conduction is carried by electrons. In the intermediate part, the temperature is not high enough for the material to be ionized, but the material is too dense for photons to carry the energy. In this part, the vibration of nuclei and molecules take over. Although it is the most effective means of non-convective energy transport in this temperature-density regime, it is still inefficient. We note that only Cond-1 includes a model for the vibrational part, but assuming a fully ionized electron conductivity, as we did in Cond-2, makes the vibrational part insignificant. In Appendix Ε we show how the uncertainty in Cond-1 affects the inferred radius of a Μ_p=10 Μ_⊕ planet.

In Figure 2, the difference in the outer part of the planet (from ~3 R_⊕) is caused by the different temperature-density regime of the radiative-opacity tables. The inner part, which is mostly dominated by the electron conductivity, is similar in shape for Cond-1 and Cond-2 but different by roughly one order of magnitude. In Cond-1, the planet has a region where the vibration of the molecules conducts most of the thermal energy. Yet, this mode is inefficient and leads to low conductivity.

Next, we analyze the temperature profiles at t = 5 Gyr in Figure 3 to explain the general energy transport and its effect on the planetary radius. In this study, the planetary radius is set to be the radius at the 1 bar pressure level. We also present the conductivity values throughout the interiors. The plots are intended to qualitatively illustrate which parts of the planet cool efficiently and which parts can contract or remain inflated. In Cond-2 the conductivity remains high throughout the planet. Therefore, energy can be transported out of the deep interior efficiently. The large convective zone carries this energy toward the outer envelope and heats it from beneath. While the interior cools, the cooling of the outer envelope is controlled by the efficiency of the radiative cooling through the atmosphere, resulting in a generally hotter outer envelope. Because the outer envelope consists mostly of Η-He, which have a high thermal expansion coefficient, the higher temperature envelope remains puffy throughout the entire evolution. On the other hand, in Cond-1 and Cond-3 the energy is efficiently trapped beneath non-convective regions with low conductivity. Both cases have hot interiors with temperatures above 10^4.4K beyond a radius of r ~ 2 R_⊕, which is a much larger region compared to Cond-2. Only a small amount of energy reaches the outer envelope, and thus the outer envelope can cool down efficiently and contract, resulting in smaller radii for Cond-1 and Cond-3 in comparison to Cond-2.

Next, we study the evolution of the planets. Figure 4 shows the radius against time for three initial energy states and a planet mass of M_p = 10 Μ_⊕ (see Figure D.1 in Appendix D for the 5 M_⊕ and 15 Μ_⊕ mass planets). We find the same phases in the evolution of a planet with a thermal boundary layer as Scheibe et al. (2021). The envelope first cools down and contracts, depending on the flux emitted throughout the atmosphere. This phase occurs on a timescale of 0.1–1 Gyr. Next, the evolution is governed by the flux from the deep interior to the outer envelope throughout conductive layers. The balance between the energy that is released throughout the atmosphere and the energy that is received from the interior determines the size of the planet (i.e., how inflated it is). Then the contraction is mainly determined by the contraction of the deep interior. However, the high-Z material has a low thermal expansion coefficient such that cooling does not lead to a large contraction and the planet evolves slowly. In the last phase, the energy reservoir of the deep interior is depleted, and contraction accelerates again.

We highlight that Cond-1 and Cond-3 are very similar during their early evolution and deviate more as time progresses. Initially, the contraction is very fast until it slows down after t ≈ 2 Gyr. The initial contraction is caused by the cooling of the outer envelope and stalls as soon as the internal heat prevents further cooling of the envelope. The planetary radius then depends on the heat flux from the deep interior to the outer envelope. Cond-2 is very different from the other two cases and has a larger radius during most of the evolution. Apart from the outer region of the planet, the conductivity in Cond-2 is higher than in the other cases. The higher conductivity leaks more energy from the deep interior to the outer envelope. This keeps the envelope hotter and therefore more extended compared to the cases with lower conductivities. As a result, the envelope is inflated throughout the entire planetary evolution in this case.

In order to understand which regions can cool and contract and which remain inflated, we show the temperature profiles at various times in Figure 5. In Cond-1, the temperature in the deep interior near the center decreases over time. Due to the slow energy transport into the envelope, the envelope contracts and increases the pressure on the layers beneath. This increasing pressure and the low thermal conductivity through the region with the composition gradient cause a slight increase in temperature around the region of r = 3 R_⊕. The higher conductivity of Cond-2 leads to faster cooling of the deep interior. The profiles of the outer envelope change little after t = 0.5 Gyr. The heat from the interior keeps the envelope more extended earlier than in the other cases.

In Cond-3, the planet cools even slower than Cond-1 in the deep interior. Since the heat remains trapped in the deep interior, the outer envelope can cool down and contract most efficiently in Cond-3. As soon as the energy of the deep interior is depleted in Cond-1 and Cond-2, the planets can contract more, while Cond-3 still fuels its slightly heated envelope. Hence, we expect that if the simulation time were much longer, Cond-3 would lead to an inflated radius over a longer period of time than Cond-1 and Cond-2.

Fig. 2 Total conductivity vs. radius at t = 5 Gyr for the warm Mp= 10 Μ⊕ planet. The blue, green, and red curves represent Cond-1, Cond-2, and Cond-3, respectively. The largest contribution to the conductivity is shown by the line styles, i.e., solid, dotted, and dashed, for radiation (opacity), vibration, and electron dominated, respectively.

Fig. 3 Temperature and conductivity as a function of the radius at t = 5 Gyr for our three cases of the “warm” planet with Mp = 10 Μ⊕. In each panel, the temperature is shown in the left half of the cone, colored in blue, red, and yellow from cold to hot. The order of magnitude of the conductivity is shown in the right half of each cone, colored in different shades of brown. The areas with circular arrows indicate convective regions within the planet. The 1 bar radius is also indicated.

Fig. 4 Radius vs. time for a planet with Mp = 10 Μ⊕. The solid, dashed, and dotted lines correspond to Cond-1, Cond-2, and Cond-3, respectively. The different colors correspond to the different initial energy states.

3.2 A narrow composition transition

Next, we investigate the effect of the conductivity value for a model with a narrow composition transition between the envelope and the interior metallicity. The structure is similar to a differentiated interior with an envelope of Η-He and heavy elements on top of a H-He-free interior. Figure 6 shows the radius evolution for the three conductivity cases for the warm 10 M_⊕ planet. As the evolution occurs faster than for the wide composition transition models, we show a logarithmic time axis for t ≤ 1 Gyr to emphasize the behavior at early times. The narrow composition region leads to a small region that is dominated by thermal transport in a non-convective region, and therefore, energy is transported much more efficiently.

Compared to the wide composition transition models we considered, Cond-1 and Cond-3 cool faster, and contraction continues throughout the entire evolution simulation. Cond-3 is slightly enlarged compared to Cond-1 for most of the evolution, and they converge to similar radii at 10 Gyr.

The high conductivity of Cond-2 allows for very efficient cooling. The initial cooling phase, where only the envelope contracts, occurs on a timescale of 0.1 Gyr. Before t = 1 Gyr, most of the internal energy is depleted such that the flux into the envelope decreases and the planet can contract again. Interestingly, the cold Cond-2 model (dashed blue) expands significantly due to the energy flux from the deep interior. We find that all the initial energy states for Cond-2 converge to a similar radius within 10 Gyr. After that, the contraction is governed by the radiative cooling of the atmosphere.

To better demonstrate the cooling behavior of the interior, Figure 7 shows the warm Cond-1 and Cond-2, and the cold Cond-2 temperature-radius profile for various times. Since Cond-3 is similar to Cond-1, it is not shown.

In the case of Cond-1, the sharp composition change leads to a clear layered thermal profile. Two large convective layers are separated by a non-convective layer that is stabilized by the composition gradient. The outer envelope is again non-convective because of efficient energy transport by radiation. After the initial contraction of the outer envelope, the slow-cooling deep interior acts as an energy reservoir that keeps the envelope inflated for several gigayears.

In Cond-2, the thermal energy transport is sufficiently high to decrease the internal energy reservoir by an order of magnitude. After most of the energy is emitted, the outer envelope can cool down and contract much more than in Cond-1 and Cond-3. As mentioned above, the radius of the cold Cond-2 model even increases in the beginning. The early profile at t = 0.01 Gyr (black lines in Figure 7) has a more contracted envelope than the slightly more evolved profile at t = 0.5 Gyr. The central temperature drops by almost 5000 Κ within these two time steps, releasing more energy into the outer parts than the atmosphere can emit. At t ~ 5 Gyr the temperature gradient across the composition transition approaches zero. The deep interior becomes non-convective, as the temperature gradient in the deep interior is lower than the adiabatic gradient. At this point, the evolution is governed by the thermal flux that escapes via the planetary atmosphere.

Fig. 5 Planetary temperature profile at various times. The three different conductivity cases are shown in the panels from left to right for the warm Mp = 10 Μ⊕ model. The color indicates the planetary age, while the style of the line indicates whether the zone is convective (solid) or non-convective (dotted).

Fig. 6 Planetary radius vs. time for a planet with Mp = 10 Μ⊕. The model is similar to the one presented in Figure 4 but has a much narrower composition transition (Δq = 0.01). The solid, dashed, and dotted lines correspond to Cond-1, Cond-2, and Cond-3, respectively. The different colors show the different initial energy states. We note that the x-axis is logarithmic between 0.01 and 1 Gyr and linear between 1 and 10 Gyr.

3.3 Vibrational conductivity model comparison

Next, we compare the models Cond-1 and Cond-4, which differ only in the treatment of the vibrational conductivity. Figure 8 shows the total conductivity and the vibrational component at t = 5 Gyr for the warm M_p = 10 M_⊕ planet. In Figure F.1 we show the radius over time for Cond-1 and Cond-4. We find only small variations in all the tested masses and initial energy states, where the largest difference is found to be 1.6% for the hot M_p = 10 M_⊕ model. The key difference lies in the smaller conductivity in Cond-4. We note that we used the same electronic conductivity in Cond-1 and Cond-4. A different model for the electronic conductivity can increase the region in the planet where the vibrational conductivity dominates and therefore increase the deviations between Cond-1 and Cond-4. We also note that the scaling law assumed in Cond-4 might not be suitable for the density range in the planet, where the vibrational conductivity dominates. This is because in this region, the densities are lower by about a factor of ten in comparison to the reference conductivity value of ρ_ref = 3.87 g cm⁻³. More information on the vibrational conductivity at various conditions is therefore important.

3.4 Implications for exoplanet characterization

We found large differences in the inferred radii depending on the model assumptions. This has important implications for the characterization of exoplanets in the relevant mass regime. Below, we analyze the results for the models with the wide composition gradient. Figure 9 shows the planetary radii at t = 5 Gyr for the first three conductivity cases; the cold, warm, and hot initial entropy states; and the three masses (see Table 2 for the exact values of the radii). The bottom panel shows the relative difference in the inferred radius when considering different conductivity models and initial entropies. In particular, we compare the warm Cond-1 with the warm Cond-2 and warm Cond-3 models. For the comparison in the inferred primordial entropy, we compare the warm Cond-1 with the cold and hot Cond-1 models. We note that all models have the same composition profile and the same bulk composition, yet the evolution still leads to very different radii. Therefore, for a given mass, the radius uncertainty is found to be of the order of ~20%, depending on the conductivity, and of the order of ~25%, depending on the primordial entropy, when considering the spread between the minimum and maximum values in Figure 9. These differences are much larger than the typical measured uncertainty of the planetary radius, which is of the order of ~5% and could go even lower with upcoming missions (e.g., PLATO). As a result, it is clear that theory plays a key role in the characterization of Neptunes and sub-Neptunes and that in order to take full advantage of the data, improvements in theory are required. In particular, making such improvements is critical to determining the conductivities of various elements at the relevant pressure-temperature regime and to better constraining the post-formation entropy (and composition gradient) of such planets. We hope that future studies will focus on these topics.

Fig. 7 Temperature vs. radius at different times for the narrow composition transition. We show the warm Cond-1 and Cond-2 in the left and middle panel and the cold Cond-2 in the right panel for the M = 10 Μ⊕ model. The color indicates the planetary age, and the solid (dashed) line indicates whether the zone is convective (non-convective).

Fig. 8 Conductivity vs. radius at t = 5 Gyr for the warm Mp = 10 M⊕ planet. The blue and brown lines correspond to Cond-1 and Cond-4. respectively. The total conductivity is shown by the solid line, and the thick line shows the contribution from the vibrational conductivity. Outside the shown radius window, the conductivity is dominated by the electron and radiation contribution.

4 Discussion

While our study presents a step forward toward the understanding of Neptunes and sub-Neptunes, some simplifying assumptions have to be made, which we discuss below. First, we note that the planetary cooling strongly depends on the atmospheric model. In this work we used a standard semi-gray irradiated atmosphere model (Guillot & Havel 2011). Our choice of atmospheric model represents a compromise between accuracy and complexity. The semi-gray approximation has been shown to deviate by ≲10% from full numerical solutions to the radiative transfer problem across a wide range of gravities and effective temperatures (Parmentier et al. 2015). However, a different model would change the inferred thermal energy flux through the atmosphere, leading to a different temperature at the bottom of the atmosphere. If this flux is significantly different compared to the one we calculate here, the radius evolution would also be affected. Since the contraction of the outer envelope depends on the energy balance between the flux coming from the deep interior and the flux that is emitted through the atmosphere, our results are still valid in terms of the relative differences between different model assumptions.

Second, we have presented low-temperature opacity tables using the AESOPUS2.1 code (Marigo et al. 2024). The new tables cover a range of 2 ≤ log(T/K) ≤ 4.5 for different metallicities and a constant Y/X = 0.34 fraction. They are an improvement of the commonly used Freedman opacities (Freedman et al. 2008) included in MESA since they cover a higher temperature range. Their disadvantage compared to the Freedman opacities is the limited range in the parameter log R = log ρ – 3 log T + 18 (all units in cgs). The Freedman tables cover a range of up to log R ≤ 9, while the AESOPUS2.1 code only supports values up to log R ≤ 6. The planet models in this study can reach values of log R ~ 7–8, and we hope to further improve on this in future work. In cases where log R is out of the range of the tables, MESA uses the last available grid point. This means that for a given temperature and density with a higher value of log R > 6, the opacity is used for the given temperature and a lower density to match the last available log R = 6.

Third, for simplicity, we kept the composition profile fixed over time and only considered the effect of the energy transport from convection. In reality, convective mixing erodes composition gradients (e.g. Knierim & Helled 2024; Tejada Arevalo 2025). Including the effect of convective mixing would affect the inferred evolution, especially for the composition profile with a wide transition region. Such a composition gradient could mix from outside inward, and the composition may start building a step-like profile with narrow non-convective regions separating convective regions. As already stated in the introduction, Uranus is suspected to have a stable non-convective boundary layer not too far inside its interior. Furthermore, if a planet forms cold (low entropy), mixing is likely to be less efficient, and stable composition gradients can be expected in many planets. Therefore, as long as there is a stable non-convective region, it is clear that the choice of the conductivity is significant and should be carefully considered. We hope to further improve on this choice in a future study that self-consistently treats convection and mixing.

Fourth, due to a lack of data, the assumed material of the EoS is not always the same as the material in our conductivity models. We took four choices for the conductivity models that all assume a different material. This material is different from the water-rock mixture we assumed for the EoS and does not scale with the metallicity. This inconsistency is hard to improve, as the data for the conductivities are sparse. We hope that further investigations of the thermal conductivity of different materials will provide the required information, allowing for more consistent planetary models.

Fifth, although we have investigated multiple cases with different masses, the initial temperature profile in the models is nearly adiabatic. However, formation models suggest that composition gradients form during the planetary formation process (Valletta & Helled 2020, 2022). Therefore, Neptunes and sub-Neptunes could have temperature profiles that differ from the adiabatic ones after their formation. The composition profile we used here is based on the idea of having a smooth transition between the central metallicity and an envelope metallicity similar to what is expected from Uranus interior models. Because mixing likely starts from the outside in, we assumed a flat composition profile for the envelope, which might have been an eroded composition gradient. We suggest that future studies include large ranges of bulk metallicities and internal structures.

Sixth, the planetary evolution depends on the choice of the EoS. In Appendix G, we compare the radius evolution obtained with different H–He EoSs and find that simulations using the SCvH EoS (Saumon et al. 1995) yield radii up to ≲2% larger after 10 Gyr compared to those with the CMS EoS. Variations in the EoS of heavy elements would also affect the results. For example, Howard et al. (2025) compared two water EoSs for static adiabatic models for Uranus, which leads to a temperature difference of ~1000 K at a few hundred kbar. Given the large uncertainty in the planetary composition and the unknown ratio between different elements (e.g., ice and rock), the uncertainty of the EOS is expected to be small. At the same time, we note that a systematic study of the uncertainties arising from both the EoS and the assumed materials is desirable.

Seventh, for simplicity, we did not consider demixing and settling, which are processes that could take place and affect the evolution of sub-Neptunes and Neptunes. For example, the phase separation of water (Bailey & Stevenson 2021; Cano Amoros et al. 2024; Howard et al. 2025) or carbon-bearing species (He et al. 2022; Cheng et al. 2023; Militzer 2024) could delay the planetary cooling and lead to radius inflation. We also note that progress in characterizing material properties of various mixtures in planetary conditions is required. We hope to address these topics in future research.

Finally, we note that for planets with a significant amount of water in the deep interior, we used the easy-to-use fits for the conductivity from French & Redmer (2017) for the electronic contribution and those of French (2019) for the vibrational contribution. Scheibe et al. (2021) employed the same models for thermal conductivity and conclude that the presence of additional materials would enhance the conductivity relative to pure water. Thus, using the conductivity of pure water provides a lower bound. For example, the higher electrical conductivity of ammonia compared to water could lead to higher thermal conductivity (Ravasio et al. 2021). It is therefore clear that in order to perform reliable planetary evolution simulations, one should use the appropriate conductivity data. We suggest that it is highly desirable to investigate the properties of different H–C–N–O mixtures at planetary conditions. Our work shows that, similar to the EoS, the conductivity of material plays an important role in planetary evolution and therefore also in planetary characterization.

Fig. 9 Planetary radius at t = 5 Gyr and relative difference between the models and initial entropies. Top: planetary radius at t = 5 Gyr for the models with the wide composition gradient. Different markers indicate the different conductivity cases, while the different colors correspond to the different initial entropies. The three planet masses (Mp = 5, 10, and 15 Μ⊕) are shown from left to right. Bottom: relative difference of the inferred radius ΔR/R = (R′ – R)/R. Green corresponds to the deviation caused by the different assumed conductivities. We compare the warm Cond-1 with the warm Cond-2 and Cond-3. Blue corresponds to the deviation caused by different initial entropy states. We compare the warm Cond-1 with the cold and hot Cond-1.

5 Summary and conclusions

Our Galaxy contains many Neptunes and sub-Neptunes, but their characterization remains a challenge. Often, these planets are modeled assuming a differentiated structure and adiabatic interior. In this work, we have investigated how the planetary thermal evolution, and therefore also the inferred internal structure, is affected when the planets consist of non-convective layers.

We find that the treatment of the stable layer, the thermal conductivity in particular, significantly affects the results. We have shown that the commonly adopted simplifications, such as electron conductivity based on fully ionized hydrogen or a constant conductivity profile, lead to very different radius predictions. In addition, we find that the vibrational contribution to the conductivity must be included. In particular, the treatment of the conductivity is critical if a significant fraction of the planet’s total energy is trapped beneath a non-convective layer. The initial energy state also has a key effect on the planetary radius and its evolution. For a given mass and composition, we find that the inferred radius can differ by up to ~20% depending on the model for the conductivity and by up to ~25% depending on the primordial energy state. Our results show that knowledge of the conductivity (of various mixtures) at planetary conditions and further constraints on the primordial state of Neptunes and sub-Neptunes are crucial for the characterization of this planetary type.

Acknowledgements

We thank the referee and editor for their helpful comments. We also thank Simon Müller and Henrik Knierim for many valuable discussions and technical support. The authors acknowledge support from SNSF under grant 200020_215634. Software: MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023), AESOPUS2.1 (Marigo & Aringer 2009; Marigo et al. 2024), Jupyter Notebook (Kluyver et al. 2016; Granger & Pérez 2021), MesaReader, NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration 2013, 2018, 2022).

Appendix A A low-temperature opacity table with AESOPUS2.1

Fig. A.1 Radiative opacity as a function of temperature for constant log R. The left, middle and right panels correspond to metallicities of Ζ = 0.02, 0.20, and 0.95, respectively. The dashed line shows the radiative opacity from Freedman et al. (2014) as included in MESA and the solid line shows the custom-made AESOPUS2.1 tables using the web interface (Marigo et al. 2024). The three different colors show different values for log R. We note that the AESOPUS2.1 tables are only available up to log R = 6, resulting in the same AESOPUS2.1 opacities for log R ≥ 6 for a given temperature (orange and green solid-line).

We use the web interface of the AESOPUS2.1 code (Marigo & Aringer 2009; Marigo et al. 2024) to create new low temperature radiative opacity tables that cover low and high metallicities. We use a fixed helium to hydrogen mass ratio of Y/X = 0.34 corresponding to the solar photosphere (Asplund et al. 2009). Following the Freedman tables included in MESA, we create 7 tables for Ζ ∈ {0.01, 0.02, 0.04, 0.1, 0.2, 0.63, 0.95}. Due to numerical limitations, we set the maximum heavy-element mass fraction to Ζ = 0.95. For the elemental abundance ratios, we select the reference solar composition of Lodders (2003).

We use the full available temperature range of 2 ≤ log Τ ≤ 4.5 in steps of Δ log T = 0.02 between 2 ≤ log Τ ≤ 3.5 and Δ log Τ = 0.05 between 3.5 ≤ log Τ ≤ 4.5. The density range is given by −8 ≤ log R ≤ 6 in steps of Δ log R = 0.2. The Web interface does not allow one to calculate the full grid in one query. Therefore, we split the temperature range into three regions and manually merge the tables. The tables provided by Freedman et al. (2014) are available up to log R = 9, while the AESOPUS2.1 tables are limited to log R = 6 at a given temperature. When a requested log R falls outside the available range, MESA’s opacity module keeps the temperature fixed and substitutes a lower density so that the inferred opacity corresponds to the highest available log R value. Figure A.1 shows the opacity as a function of the temperature for constant log R and compare our new table with the Freedman opacity. Marigo et al. (2024) suggested that the differences between their results and those of Freedman et al. (2014) could be due to pressure broadening and a lower number of molecules considered for molecular transitions.

Appendix B Vibrational conductivity

Ross et al. (1984) reviewed the phonon model for a completely anharmonic crystal and derived the temperature and density scaling of the thermal conductivity of liquids under high pressure. We follow the formulation of Stamenković et al. (2011), such that the partial derivatives of the vibrational conductivity k_vib are given by $${\left(\frac{\partial \ln k_{\text{vib}}}{\partial \ln T}\right)}_{\rho }=-1,$$(B.1) $$\psi :={\left(\frac{\partial \ln k_{\text{vib}}}{\partial \ln \rho }\right)}_T=3\gamma -2{\left(\frac{\partial \ln \gamma }{\partial \ln \rho }\right)}_T-\frac{1}{3},$$(B.2)

with γ = (∂ ln v/∂ In ρ)_T being the Grüneisen parameter for the average vibrational frequency v. The parameter γ can be approximated using the thermodynamic Grüneisen parameter $${\gamma }_{\text{TH}}=\frac{\alpha K_T}{\rho c_V},$$(B.3)

where α = – (∂ ln ρ/∂ ln T)_P is the thermal expansivity, K_T = – (∂ ln ρ/∂ ln P)_T the isothermal compressibility, and c_V the specific heat capacity at constant volume. Given a reference conductivity k_{vib, ref}(T_ref, ρ_ref) and an EOS to calculate γ_ΤΗ, one can construct the vibrational conductivity from $$k_{\text{vib}}=k_{\text{vib},ref}{\left(\frac{\rho }{{\rho }_{\text{ref}}}\right)}^{\psi }\left(\frac{T_{\text{ref}}}{T}\right),$$(B.4)

with ψ defined by eq. (B.2) using the approximation γ ≈ γ_ΤΗ.

Appendix C Initial model

In this section, we describe how we use MESA to produce the initial models. The basic outline is to first create an adiabatic model and utilize MESA’s relax options to create the desired initial state. The procedure ensures that the initial model has a given mass M_{p, i}, a specific composition profile Z_i(q) as a function of the normalized mass coordinate q, and an initial energy budget defined by some central specific entropy s_{central, i}- If not stated otherwise, we let MESA run for 1000 years with maximal time steps of 100 years, between each step, to allow the model to settle in for better convergence in the long run.

Step 1: We create an adiabatic model using the option create_initial_model with M_p = 20 Μ_⊕, R_p = 9.4R_⊕, Z_bulk = 0.8, Y/X = 0.34. Note, that this model has a flat composition profile. We force the model to follow the adiabat given by the EoS.

Step 2: Next, we relax the composition to Z_bulk = 0.94 using relax_initial_Z to avoid convergence issues for planets with masses of M_p ≲ 8 Μ_⊕. Ideally, one would set the metallicity in this step to the final central metallicity to ensure a consistent entropy, which is set in Step 4. However, this led to issues, which is why we set Z_bulk to a lower value. This is accounted for by setting a slightly higher entropy in Step 4. The model is forced to follow the adiabat.

Step 3: Now we use the relax option relax_initial_mass to shrink the current model down to its desired mass M_{p, i}. This step can be difficult if M_{p, i} is low. In this case, increasing the metallicity as described in Step 2 can help. Or before changing the mass one can set a different entropy profile as described in Step 4. Again, we force the model to follow the adiabat.

Step 4: We set the initial energy budget, which we parametrize by the central specific entropy s_{center, i}. Using the option relax_initial_entropy MESA applies an extra specific heat Q_l in each layer l according to $$Q_l=\frac{1-s_l}{s_{\text{i}}}\frac{E_l}{{\tau }_{\text{relax\hspace{0.17em}entropy}}},$$(C.1)

with s_l being the current specific entropy, Ε_l the current specific energy and τ_{relax entropy} = 1 × 10⁻³ yr an relax timescale. Usually, the adiabat, i.e., a constant entropy profile, agrees with the profile that is predicted by the mixing length theory. But especially for the outer parts the entropy profile can significantly deviate from an adiabat. This can lead to convergence issues where the solver reacts with rapidly changing profiles. We avoid this issue by forcing the model to follow the adiabatic gradient given by the EoS. The specific entropy depends on the metallicity, and our final composition profile features a central metallicity of Z_{int, i} = 1, which differs from the value set in step Step 2. As a result, in this step we assign a higher entropy such that the initial model after relaxing the composition has the desired entropy.

We construct three initial energy states — hot, warm, and cold — corresponding to central entropies of s_{center, i} = 0.7, 0.6 and 0.5 k_B/m_H respectively. Although models with different masses and narrow composition gradients may show slight variations in s_{center, i} they agree to within rounding to the first decimal place.

Fig. C.1 Initial profiles for the narrow composition gradient, showing specific entropy, temperature, and composition as functions of normalized mass. The colors blue, orange, and yellow correspond to different primordial entropies. The dotted, solid, and dashed lines correspond to planets with a mass of 5 Μ⊕, 10 Μ⊕, and 15 Μ⊕, respectively. Top: Specific entropy vs. normalized mass of the initial model for the narrow composition profile. The blue (cold), orange (warm), and yellow (hot) colors represent different initial entropies. The dotted, solid, and dashed lines correspond to planets with a mass of 5 Μ⊕, 10 Μ⊕, and 15 M⊕, respectively. The metallicity of the composition profile is represented by a black line, with its y-axis on the right-hand side. Bottom: Temperature vs. normalized mass of the initial model for the narrow composition profile. The black line shows the metallicity.

Step 5: Next, we set the composition profile using the option relax_initial_composition. MESA relaxes the composition profile of the planet by reading in a file that contains the normalized mass coordinates for each layer with the desired mass fraction of each species. We use a reduced isotope network that only includes hi, h4, and ο 16 representing X, Y, and Z, respectively. This step sometimes leads to convergence issues. It is advised to keep the parameter num_steps_to_relax_composition not too high and the relax timescale timescale_for_relax_composition short. We find values with 50 and −1 to be well-behaved. Because the entropy is a function of temperature, density, and composition, the final profile will not have the flat entropy profile we set in Step 4. The initial specific entropy profiles and composition profiles (after the completion of the four steps described above) are shown in the upper panels of Figure 1 for the wide gradient. Figure C.1 shows the entropy and temperature profiles for the narrow gradient. In this step, the model is not forced to follow the adiabatic gradient and the outer boundary condition is given by the default photosphere without irradiation.

We consider two cases for the composition profile, one with a wide transition region (transition over Δq = 0.1) and one with a narrow transition region (transition over Δq = 0.01). Both profiles are fixed to the same bulk metallicity of Z_bulk = 0.95 with an envelope metallicity of Z_env = 0.2.

Step 6: Using the options irradiation_flux we heat the outer parts of the planet over 100 000 years sufficiently high such that changing to a different atmosphere model does not lead to convergence issues. We find a flux of 8 × 10⁶ erg/s/cm² and a column depth of 100 g/cm² a good choice of parameters for our desired equilibrium temperature.

Step 7: Finally we change the atmosphere setting atm_option to ’irradiated_grey’. The initial model is then evolved for another 10 000 years before saving it as a . mod file.

This process is repeated for the three initial energy states and the two composition profiles. The initial temperature profiles are shown in the lower panels of Figure 1 for the wide gradient, and in Figure C.1 for the narrow gradient. For all initial models, we use the thermal conductivity model corresponding to Cond-1.

Appendix D Radius evolution for the small and large planet

For completeness, in this section we show the radius evolution for the three initial entropy cases discussed in Section 3.1, focusing on planets with 5M_⊕ and 15 M_⊕, as illustrated in D.1.

Fig. D.1 Radius vs. time for the planets with Mp = 5 Μ⊕ (top panel) and Mp = 15 Μ⊕ (bottom panel). The solid, dashed, and dotted lines correspond to Cond-1, Cond-2, and Cond-3, respectively. The different colors represent the different initial energy states.

Appendix E Effects of Cond-1 uncertainties

We investigate how the uncertainty in the conductivity inferred with model Cond-1 affects the planetary evolution. French & Redmer (2017) state that the systematic error is expected to be up to 5% and suggest that the statistical error can be estimated based on the scattering of the data points around the fit (see Fig. 4 in French and Redmer 2017). Therefore, we take an uncertainty for the vibrational and electronic conductivity of 30% in Cond-1 and scale the conductivities by a factor of k′ = k · (1 ± 0.3), where k is the unsealed conductivity of the electronic or vibrational part. Figure E.1 shows the effect of the uncertainty in conductivity on the radius evolution for the 10 M_⊕ planet. We find small deviations, with the highest being −2% for the hot 15 Μ_⊕ planet. Generally, the deviations are larger for the hotter models because more energy is trapped within the planet.

Fig. E.1 Planetary radius vs. time for a planet with Mp = 10 Μ⊕. The shaded area shows the uncertainty originating from the statistical error of the conductivity model in Cond-1. The different colors correspond to the different initial energy states. We annotate the final radius and the uncertainty associated with the chosen conductivity model.

Appendix F Radius evolution comparison for Cond-1 and Cond-4

We compare the radius evolution between Cond-1 and Cond-4 for the M_p = 10 Μ_⊕ planet in F.1.

Fig. F.1 Planetary radius vs. time for a planet with Mp = 10 Μ⊕. The solid and dashdot line indicates Cond-1 and Cond-4, respectively. The different colors correspond to the different initial energy states.

Appendix G The effect of the Η-He EoS

In this study we used the CMS EoS with non-ideal interactions for Η-He (see Section 2 for details). Here, we present the planetary evolution when using the SCvH EoS for Η-He (Saumon et al. 1995) instead of CMS. Figure G.1 shows the radius evolution for the models with a warm start and a wide composition gradient using both EoS. We find small deviations of ≲ 2% for the planetary radius after 10 Gyr. Comparing the water EoS using MESA is currently not possible due to the rather complex implementation procedure (e.g. Helled et al. 2025).

Fig. G.1 Planetary radius vs. time for the models with a wide composition gradient and warm initial conditions. The solid and dashed line correspond to the CMS and SCvH EoS for Η-He, respectively. The colors, blue, green and red show different planetary masses of Mp =5, 10, and 15Μ⊕, respectively. For each simulations we show the radius after t = 10 Gyr on the right hand side.

As noted above, detailed investigations of different EoS s and mixtures that systematically quantify the role of the used EoS and assumed composition are required. This is in particularly timely now as exoplanetary data is improving and theoretical uncertainties can be significant (Müller et al. 2020a).

References

All Tables

All Figures

Fig. 1 Initial profiles for the wide composition gradient, showing specific entropy, temperature, and composition as functions of normalized mass. The colors blue, orange, and yellow correspond to different primordial entropies. The dotted, solid, and dashed lines correspond to planets with a mass of 5 Μ⊕, 10 Μ⊕, and 15 Μ⊕, respectively. Top: Specific entropy vs. normalized mass of the initial model for the wide composition profile. The metallicity of the composition profile is shown by a black line, with its y-axis on the right-hand side. Bottom: Temperature vs. normalized mass of the initial model for the wide composition profile. The black line corresponds to the metallicity.

In the text

	Fig. 2 Total conductivity vs. radius at t = 5 Gyr for the warm Mp= 10 Μ⊕ planet. The blue, green, and red curves represent Cond-1, Cond-2, and Cond-3, respectively. The largest contribution to the conductivity is shown by the line styles, i.e., solid, dotted, and dashed, for radiation (opacity), vibration, and electron dominated, respectively.
In the text

Fig. 3 Temperature and conductivity as a function of the radius at t = 5 Gyr for our three cases of the “warm” planet with Mp = 10 Μ⊕. In each panel, the temperature is shown in the left half of the cone, colored in blue, red, and yellow from cold to hot. The order of magnitude of the conductivity is shown in the right half of each cone, colored in different shades of brown. The areas with circular arrows indicate convective regions within the planet. The 1 bar radius is also indicated.

In the text

	Fig. 4 Radius vs. time for a planet with Mp = 10 Μ⊕. The solid, dashed, and dotted lines correspond to Cond-1, Cond-2, and Cond-3, respectively. The different colors correspond to the different initial energy states.
In the text

	Fig. 5 Planetary temperature profile at various times. The three different conductivity cases are shown in the panels from left to right for the warm Mp = 10 Μ⊕ model. The color indicates the planetary age, while the style of the line indicates whether the zone is convective (solid) or non-convective (dotted).
In the text

Fig. 6 Planetary radius vs. time for a planet with Mp = 10 Μ⊕. The model is similar to the one presented in Figure 4 but has a much narrower composition transition (Δq = 0.01). The solid, dashed, and dotted lines correspond to Cond-1, Cond-2, and Cond-3, respectively. The different colors show the different initial energy states. We note that the x-axis is logarithmic between 0.01 and 1 Gyr and linear between 1 and 10 Gyr.

In the text

	Fig. 7 Temperature vs. radius at different times for the narrow composition transition. We show the warm Cond-1 and Cond-2 in the left and middle panel and the cold Cond-2 in the right panel for the M = 10 Μ⊕ model. The color indicates the planetary age, and the solid (dashed) line indicates whether the zone is convective (non-convective).
In the text

	Fig. 8 Conductivity vs. radius at t = 5 Gyr for the warm Mp = 10 M⊕ planet. The blue and brown lines correspond to Cond-1 and Cond-4. respectively. The total conductivity is shown by the solid line, and the thick line shows the contribution from the vibrational conductivity. Outside the shown radius window, the conductivity is dominated by the electron and radiation contribution.
In the text

Fig. 9 Planetary radius at t = 5 Gyr and relative difference between the models and initial entropies. Top: planetary radius at t = 5 Gyr for the models with the wide composition gradient. Different markers indicate the different conductivity cases, while the different colors correspond to the different initial entropies. The three planet masses (Mp = 5, 10, and 15 Μ⊕) are shown from left to right. Bottom: relative difference of the inferred radius ΔR/R = (R′ – R)/R. Green corresponds to the deviation caused by the different assumed conductivities. We compare the warm Cond-1 with the warm Cond-2 and Cond-3. Blue corresponds to the deviation caused by different initial entropy states. We compare the warm Cond-1 with the cold and hot Cond-1.

In the text

Fig. A.1 Radiative opacity as a function of temperature for constant log R. The left, middle and right panels correspond to metallicities of Ζ = 0.02, 0.20, and 0.95, respectively. The dashed line shows the radiative opacity from Freedman et al. (2014) as included in MESA and the solid line shows the custom-made AESOPUS2.1 tables using the web interface (Marigo et al. 2024). The three different colors show different values for log R. We note that the AESOPUS2.1 tables are only available up to log R = 6, resulting in the same AESOPUS2.1 opacities for log R ≥ 6 for a given temperature (orange and green solid-line).

In the text

Fig. C.1 Initial profiles for the narrow composition gradient, showing specific entropy, temperature, and composition as functions of normalized mass. The colors blue, orange, and yellow correspond to different primordial entropies. The dotted, solid, and dashed lines correspond to planets with a mass of 5 Μ⊕, 10 Μ⊕, and 15 Μ⊕, respectively. Top: Specific entropy vs. normalized mass of the initial model for the narrow composition profile. The blue (cold), orange (warm), and yellow (hot) colors represent different initial entropies. The dotted, solid, and dashed lines correspond to planets with a mass of 5 Μ⊕, 10 Μ⊕, and 15 M⊕, respectively. The metallicity of the composition profile is represented by a black line, with its y-axis on the right-hand side. Bottom: Temperature vs. normalized mass of the initial model for the narrow composition profile. The black line shows the metallicity.

In the text

	Fig. D.1 Radius vs. time for the planets with Mp = 5 Μ⊕ (top panel) and Mp = 15 Μ⊕ (bottom panel). The solid, dashed, and dotted lines correspond to Cond-1, Cond-2, and Cond-3, respectively. The different colors represent the different initial energy states.
In the text

	Fig. E.1 Planetary radius vs. time for a planet with Mp = 10 Μ⊕. The shaded area shows the uncertainty originating from the statistical error of the conductivity model in Cond-1. The different colors correspond to the different initial energy states. We annotate the final radius and the uncertainty associated with the chosen conductivity model.
In the text

	Fig. F.1 Planetary radius vs. time for a planet with Mp = 10 Μ⊕. The solid and dashdot line indicates Cond-1 and Cond-4, respectively. The different colors correspond to the different initial energy states.
In the text

	Fig. G.1 Planetary radius vs. time for the models with a wide composition gradient and warm initial conditions. The solid and dashed line correspond to the CMS and SCvH EoS for Η-He, respectively. The colors, blue, green and red show different planetary masses of Mp =5, 10, and 15Μ⊕, respectively. For each simulations we show the radius after t = 10 Gyr on the right hand side.
In the text