© The Authors 2026

1 Introduction

Sub-Neptunes and Neptunes are among the most common planetary types in the current sample of over 6000 detected exoplanets. A key objective of exoplanetary science is to understand their internal structure and origin by measuring their mass, radius, age, and atmospheric composition. However, determining their composition is challenging because multiple solutions can fit the same observational data. To reduce this degeneracy in bulk composition and internal structure, it is important to connect results from planet formation theories with planetary evolution simulations, and then link these to observations representing the planet’s current-state.

Exoplanet interiors are often modeled as separate layers of distinct composition. Each layer is assumed to be homogeneous in composition and fully convective, allowing the temperature profile to be approximated using the adiabatic gradient. However, several arguments challenge this layered structure. For instance, the water distribution within planets is still under debate. Under high pressure water can preferentially dissolve into the core (Luo et al. 2024), while also being produced at a magma hydrogen interface (Horn et al. 2025). Additionally, advanced formation models indicate that composition gradients naturally arise from the accretion of solids and hydrogen–helium gas (e.g., Ormel et al. 2021; Valletta & Helled 2022). Other processes taking place during the planetary evolution such as rain-out and demixing can also lead to the build up of boundary layers (e.g., Piaulet-Ghorayeb et al. 2025; Tejada Arevalo et al. 2026). Indeed, it was proposed that water rains out of the hydrogen-dominated envelope inside Uranus. This creates a sharp transition between a water-rich and water-poor layer (Cano Amoros et al. 2024; Howard et al. 2025). Such composition gradients can inhibit large-scale convection(Stevenson & Salpeter 1977; Guillot 1995; Markham et al. 2022). Recently, randomly generated density profiles of Uranus and Neptune that match the available gravity data show large regions without convection (Morf & Helled 2025). This supports the claim that planets are not fully convective. Therefore, the temperature gradient is expected to deviate from the adiabatic gradient.

In non-convective regions, the energy transport depends on thermal conduction and radiation. In the outer envelope, where the density is small, most of the energy is carried by thermal photons that diffuse throughout the planet. This process is usually modeled using tables for the Rosseland mean opacity (Ferguson et al. 2005; Freedman et al. 2014; Marigo et al. 2024). At higher densities and temperatures, photon diffusion is inefficient and thermal conductivity dominates. The thermal conductivity is often modeled by accounting for free electrons. However, the atomic nuclei can also contribute to the thermal conductivity, where the contribution depends on the exact planetary composition (Ross et al. 1984; Stamenković et al. 2011; French 2019). For example, it was shown that neglecting the nuclei contribution in water significantly underestimates the thermal conductivity (French 2019). We summarize the contribution from the nuclei or lattice vibrations as “vibrational conductivity”.

In Eberlein & Helled (2025) (erratum Eberlein & Helled 2026) (hereafter Paper I), we simulated the evolution of sub-Neptunes and Neptunes assuming that they have composition gradient that is stable against convection. For these types of planets, the composition gradient is expected to be relatively close to the atmosphere, i.e., above most of their internal energy. Hence, the treatment of the energy transport is crucial for their evolution. We compared four commonly used approaches to model the thermal conductivity. We found significant deviations in the thermal evolution of sub-Neptunes and Neptunes, which led to a radius difference of up to 20% depending on the chosen model. Furthermore, we showed that a layer with low conductivity high up in the planet makes the initial entropy state in the deep interior more important. A hotter start can inflate the radius over a long period of time. However, the results presented in Paper I neglect the possibility of convective mixing and, therefore, the change of the internal structure with time. The objective of this study is to include this effect and then test the importance of the thermal conductivity in a self-consistent model when mixing is considered.

Our paper is organized as follows. In Section 2, we describe the key aspects of the model. In Section 3, we present the results, and we address double-diffusive convection in Section 4. We discuss our results in Section 5, and we present our summary and conclusion in Section 6.

2 Methods

We followed the general procedure presented in Paper I. To model the evolution of sub-Neptunes and Neptunes we solved the stellar evolution equations (e.g., Kippenhahn et al. 2013) using the code Modules for Experiments in Stellar Astrophysics (MESA) (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). For the hydrogen–helium equation of state (EoS) we used the tables first presented in Müller et al. (2020b,a) in the updated version (Müller & Helled 2021, 2024), which includes nonideal interactions (Chabrier & Debras 2021). For the heavy-element EoS we assumed a fixed 50/50 water-to-rock ratio represented by H₂O and SiO₂ tables (More et al. 1988; Vazan et al. 2013). The tables were implemented in MESA based on a modified version of the extension custom EoS (Knierim & Helled 2024; Helled et al. 2025).

For the atmosphere boundary condition we used the irradiated gray model (Guillot & Havel 2011) with a ratio between the visible and thermal opacity of κ_V/κ_th = 0.03 based on the fit provided by Poser & Redmer (2024). We assumed an equilibrium temperature of T_eq = 400 K. A detailed discussion of the model and its simplifications including our choice of atmosphere, discrepancy between the assumed material for EoS and thermal conductivity, uncertainty in the initial conditions, choice of EoS, and demixing can be found in Paper I.

2.1 Thermal conductivity and opacity

In non-convective regions, the thermal transport depends on the total thermal conductivity k_tot: $$k_{\mathrm{t}\mathrm{o}\mathrm{t}}=k_{\mathrm{r}\mathrm{a}\mathrm{d}}+k_{\mathrm{v}\mathrm{i}\mathrm{b}}+k_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{c}},$$(1)

which is the sum of the conductivity contribution from photons, k_rad; vibrations in a dense fluid, k_vib (lattice or nuclei); and electrons, k_elec. The radiative opacity is related to the radiative conductivity by $$K_{\mathrm{r}\mathrm{a}\mathrm{d}}=\frac{16\sigma }{3}\frac{T^3}{\rho k_{\mathrm{r}\mathrm{a}\mathrm{d}}},$$(2)

with σ being the Stefan–Boltzmann constant, T the temperature, and ρ the density. For the vibrational and electron conductivity, we used three different models. We summarize the three models in Table 1 and the valid data region of each contribution in Figure 1. For the radiative opacity, k_rad, we used tables calculated with the AESOPUS2.1 web interface (Marigo & Aringer 2009; Marigo et al. 2024; Eberlein & Helled 2025). In Appendix A we describe how we treated high-density values that fall outside the tabulated area. For the vibrational conductivity, we used an empirical fit to the ab initio calculated conductivities of water (Eq. (7) in French 2019). For the electron conductivity, we used three different models: the results for partially ionized water (Eq. (30) in French & Redmer 2017); the fully ionized electron conductivity implemented in MESA (Cassisi et al. 2007, privately communicated by A.Y. Potekhin); and a constant electron conductivity with a value that is expected for Earth’s core–mantel boundary (e.g., Stevenson et al. 1983; Lobanov et al. 2021). In Paper I, we had a fourth conductivity model, which included the vibrational conductivity model based on Stamenković et al. (2011) for MgSiO₃. We concluded that this model is less relevant because the conductivity contribution from water is much higher and therefore more dominant. Additionally, the reference conductivity of MgSiO₃ is given at ρ = 3.87 g/cm³ and T = 2000 K, which is at least ~1000 K below the temperatures reached in most of our models at similar densities.

Fig. 1 Boundaries of different conductivity sources in density-temperature space. The red area (krad AESOPUS) shows the tabulated region from the AESOPUS2.1 tables (Marigo et al. 2024; Eberlein & Helled 2025), and the dark blue (kelec FR17) and light blue (kvib F19) area shows the region with ab initio data points for partially ionized water (French & Redmer 2017; French 2019). The green (kelec C07) area shows the tabulated electron conductivity as implemented in MESA (Cassisi et al. 2007). Hatched regions with squared lines (kvib dominated) or dots (kelec dominated) indicate where kvib or kelec contribute more than 50% to the total conductivity assuming the Cond-1 model with Z = 20 for the radiative conductivity. The density temperature regimes around A, B, C, and D are of particular interest and are further discussed in Section 5.

2.2 Convection

We considered convective mixing and changed the composition in convective regions accordingly. A region is unstable against large-scale convection if the Ledoux criterion is fulfilled: $${\mathrm{\nabla }}_{\mathrm{T}}\ge {\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}+\frac{\phi }{\delta }{\mathrm{\nabla }}_{\mu },$$(3)

with the temperature gradient ${\mathrm{\nabla }}_{\mathrm{T}}=\frac{\mathrm{d}\ln T}{\mathrm{d}\ln P}$, the adiabatic temperature gradient ∇_ad, the mean-molecular weight gradient ${\mathrm{\nabla }}_{\mu }=\left(\frac{\mathrm{d}\ln \mu }{\mathrm{d}\ln P}\right)$, and the thermodynamic derivatives $\delta =-{\left(\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }\ln T}\right)}_{P,\mu }$ and $\phi ={\left(\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }\ln \mu }\right)}_{P,T}$ (Ledoux 1947; Kippenhahn et al. 2013). Given a local luminosity, l, and pressure, P, the logarithmic temperature gradient by thermal conduction ${\mathrm{\nabla }}_{\text{th}}={\left(\frac{\mathrm{d}\ln T}{\mathrm{d}\ln P}\right)}_{\text{th}}$ is related to the conductivity by $${\mathrm{\nabla }}_{\mathrm{t}\mathrm{h}}=\frac{1}{4\pi }\frac{lP}{\rho T{k}_{\mathrm{t}\mathrm{o}\mathrm{t}}}.$$(4)

In a stable region we set ∇_T = ∇_th. As a result, the possibility of double-diffusive convection (semi-convection) was not considered (see detailed discussion on double-diffusive convection in Section 4). We assumed a mixing length of α_MLT = 0.1H_P, where H_P is the pressure scale height. Simulating convective mixing in 1D hydrostatic codes is often challenging. We used the extension gentle_mixing to MESA (Knierim & Helled 2024; Helled et al. 2025) to limit the maximal change in the composition profile between two time steps. The extension controls the time-step and convective-diffusion parameters to avoid sudden jumps in the composition that lead to convergence issues. We set the maximum squared difference of two consecutive composition profiles to $h_i^2={\int }_0^M{(X_i^{\mathrm{\prime }}-X_i)}^2\mathrm{d}m/M=0.001$, where X_i and $X_i^{\mathrm{\prime }}$ are mass fractions of the specie i before and after a time step, respectively. Furthermore, we limited the maximum time step to Δt = 10⁵ years during the first t = 0.1 billion years.

2.3 Initial model

We considered three different planetary masses of M_p = 5, 10, 15 M_⊕. The heavy-element mass fraction was set to Z_env = 0.20 in the envelope, and to Z_int = 1.0 in the deep interior. Both regions were connected by a transition region with a smooth composition gradient. We considered three different composition gradients: (i) a wide transition that extends over a region with a mass of Δq = 0.1M_p, (ii) a medium transition of Δq = 0.01M_p, and (iii) a narrow transition of Δq = 0.001M_p. In all the models, the total heavy-element mass fraction in the planet is Z_bulk = 0.95. The wide and narrow composition gradients are the same as in Paper I¹. We investigate the dependence of the evolution on the gradient width in further detail in Appendix B.

We used three different initial specific entropies, which we refer to as “cold”, “warm”, and “hot”, respectively. Contrary to Paper I, we used an increasing entropy gradient toward the surface to aid numerical stability in the convective mixing phase. This led to the configuration of a stable planetary interior, with convective zones growing predominately inwards from the surface rather than emerging within the composition gradient. Before setting the composition profile we set the initial entropy s_i(q) as a function of the normalized mass coordinate q = m/M_p using the relation $$s_{\mathrm{i}}(q)=s_{\text{center},\mathrm{i}}+q\mathrm{\Delta }s,$$(5)

with the initial central entropy s_center,i and the slope set to Δs = 0.1 k_B/m_u. The parameter s_center,i was chosen such that s_i (q = 0.8) = 0.5, 0.6, 0.7 k_B/m_u. Having used this approach, the initial entropy slightly below the composition profile (at q = 0.8) is thus the same as in Paper I. The initial models were created using the relax options of MESA in multiple steps (see Paper I for further details). Figure 2 shows the initial entropy and temperature profile for the wide composition gradient models.

Fig. 2 Initial profiles for the wide composition gradient, showing specific entropy (top) and temperature (bottom) as a function of normalized mass. The heavy-element mass fraction is overlaid in black in both panels, with its axis on the right. Blue (cold), orange (warm), and yellow (hot) correspond to different primordial entropies. The dotted (5 M⊕), solid (10 M⊕), and dashed (15 M⊕) lines indicate different planetary masses.

Fig. 3 Thermal conductivity as a function of temperature and density (top row) and as a function of density along the initial and final planet profile (bottom row). Top: total conductivity heat map as a function of density and temperature for the three different conductivity models (from left to right): Cond-1, Cond-2, and Cond-3. The color represents the total conductivity value. Black lines show an example temperature-density profile for a planet with Mp = 10 M⊕ under the hot start scenario, calculated with the respective conductivity model: AESOPUS (Marigo et al. 2024; Eberlein & Helled 2025), FR17 (French & Redmer 2017), F19 (French 2019), and C07(Cassisi et al. 2007). The upper line (Initial profile) corresponds to the initial state of the evolution, while the lower line (Final profile) represents the state after 10 Gyr. The dotted line indicates a convective region within the planet. Hatched regions with diagonal lines (radiation dominated), squared lines (vibration dominated), and dots (electron dominated) indicate where a specific conductivity mechanism contributes more than 50% to the total conductivity. Colored lines with labels mark the boundaries of the conductivity models similarly to Figure 1. For this plot, we assumed Z = 0.20 for the radiative conductivity and Z = 1.0 for the electron conductivity of C07. We note that the distribution of heavy elements changes within the planetary interior. Bottom: conductivity as a function of density for the initial (gray) and final (black) profiles shown in the top panel. The thick solid line shows the total conductivity along the density–temperature profile of the planet, while the thin solid (radiative), dashed (vibrational), and dotted (electron) lines show the individual conductivity contributions.

3 Results

3.1 Comparison of conductivity models

Figure 3 shows the total conductivity using each conductivity model. To illustrate which density temperature regions are important for the evolution of planets we overlay example profiles at the beginning of the simulation and after ten billion years. We note that these models are for illustration purposes and vary for different input parameters and assumptions. For example, much colder planets or faster cooling planets will enter regions with lower temperatures. The colored lines indicate the valid regions of the models, outside these regions the models are extended.

For Cond-1, the vibrational conductivity divides the region that is dominated by radiation and electron conduction. Given its extensive coverage in the density-temperature space, it cannot be neglected. Therefore, a planet that contains significant amounts of water should have a large region where non-convective energy transport is dominated by the vibrational conductivity. However, the fits of the vibrational conductivity were extended beyond the recommended range by the authors (see ρ ≲ 0.2 g/cm³ and T ≳ 4000K, and ρ ≳ 2 g/cm³ and T ≳ 1000K). In these extrapolated regions no data points validate the fit, yet the predicted values exceed those of the radiative conductivity of the AESOPUS2.1 tables. Electron conductivity dominates only at high densities and high temperatures. A significant part of the planet is outside the region that is covered by any of the tabulated data. As the planet cools, this region becomes larger.

For Cond-2, the electron conductivity dominates a larger part of the temperature-density space. In particular, the vibrational conductivity of Cond-1 is insignificant compared to the electron conductivity of Cond-2. The electron conductivity of Cassisi et al. (2007) is tabulated for T ≥ 10³ K. Hence, in the case of Cond-2, the part of the planet around ρ ~ 10⁻³ g/cm³ reaches temperature-density values where neither the electron conductivity nor the radiative conductivity is valid.

In the case of Cond-3, the electron conductivity dominates at higher densities, where radiation transport is ineffective. Again, with this approach, a large part of the planet is modeled outside the original AESOPUS2.1 table, and conductivities must be extrapolated.

3.2 Convective mixing

The key difference from Paper I is the consideration of convective mixing. As the planet cools, the convective region moves from outside in and erodes the composition gradient. This increases the heavy-element mass fraction in the envelope and can create convective staircases (e.g., Vazan et al. 2018; Müller et al. 2020b; Knierim & Helled 2024; Tejada Arevalo et al. 2025). These staircases should not be confused with the process of double-diffusive convection (Garaud 2018), which is not included in our simulations. The upper panel of Figure 4 shows the composition profile at different times. Convective zones appear in the region that was previously a smooth composition gradient. Small non-convective regions separate the convective zones with jumps in the heavy-element mass fraction. We find that this process is very sensitive to the used model assumptions. Nevertheless, we can identify certain trends. If the planet starts hotter, fewer and larger steps appear during the planetary evolution. If the planet starts very cold, some of the composition gradient remains with a large non-convective region.

Furthermore, the thermal conductivity influences when a region becomes unstable against convection. To illustrate this, we show the final composition profile using the different conductivity models in the lower panel of Figure 4. Because Cond-2 has a higher conductivity, the temperature gradient in non-convective regions is shallower (see Equation (4)). It follows that the shallow temperature gradient in non-convective regions leads to higher stability against convection (see Equation (3)). Therefore, the composition profile does not develop a convective staircase throughout the entire composition gradient. The lower conductivities of Cond-1 and Cond-3 lead to steeper temperature gradients that fulfill Equation (3) such that staircases appear.

Fig. 4 Heavy-element mass fraction versus radius at different times for the Cond-1 model (top) and for the three different models at t = 5 Gyr (bottom). Solid lines indicate non-convective regions, while dashed lines indicate convective regions. This plot corresponds to the simulations of a planet with Mp = 10 M⊕, a wide composition gradient, and a warm start.

Fig. 5 Radius over time with and without mixing. The simulations with mixing are shown in dark colors, while the simulations without mixing are shown in light colors. The solid lines represent the models with a wide composition gradient, while the dotted lines represent the model with a narrow composition gradient. The blue color indicates a cold start, while the yellow color indicates a hot start. The simulations for this plot assume Mp = 10 M⊕, a wide gradient, and Cond-1.

3.3 Radius evolution

Next, we compare the radius evolution with and without mixing. Figure 5 shows the radius evolution for different initial entropies and a wide and a narrow composition gradient. For a cold start the entropy is not high enough to erode the composition gradient. In this case the composition profile remains mostly intact, and therefore the difference between the simulations with and without mixing are low. The narrow gradient is too steep to be eroded by convection, and a convective staircase does not appear. Hence, a thermal boundary layer remains between the envelope and the deep interior such that the radius evolves similarly to the case without mixing. For the hot start with a wide gradient the entropy is high enough to erode some of the composition gradient and to create a composition staircase. This significantly increases the thermal transport such that the interior can effectively cool down. The radius contracts over the entire evolution.

In Figure 6 we show the radius evolution for a planet with M_p = 10 M_⊕ including convective mixing (see Appendix C for M_p = 5 M_⊕ and M_p = 15 M_⊕). As seen before if the planet starts cold the radius evolution is similar as in Paper I. Trapped heat creates smaller radii for Cond-1 and Cond-3 compared to Cond-2. The deviation from the previous paper arises at higher initial entropies. In case of the warm scenario, the Cond-2 model has a short phase where the radius increases at t ≈ 0.5 Gyr. This increase happens when the core becomes convective and the deep interior can rapidly lose energy. We find this behavior in most of the warm and hot initial models. In the case of the warm Cond-3 model, the low conductivity gives rise to convective instabilities. The resulting composition profile contains few steps with large convective zones (see Figure D.1). This results in an efficient energy transport through the staircase such that the radius evolution is determined by the cooling through the atmosphere. The hot models start much more inflated but cool down over the lifetime of the planet and eventually evolve to similar radii. Within the first ~100 Myr, the mixing process is complete and all hot models develop convective staircases (see Figure D.1). The small non-convective regions between these staircases are not large enough to effectively trap heat. The radius evolution is therefore mostly governed by the cooling through the atmosphere. This explains why the different conductivity models result in a similar radius evolution.

Figure 7 shows the planetary radius at t = 5 Gyr for different masses and conductivities, both with and without convective mixing. The results demonstrate that for high initial entropies, the thermal conductivity is less relevant for the late radius evolution (after a few gigayears). This is because convective mixing can transform the composition gradient into staircases, where heat trapping is less efficient. Interestingly, for a M_p = 15 M_⊕ planet, the hot Cond-2 model is stable against convection, and therefore the planet preserves most of its primordial composition profile (see Figure D.1). The lower panel of Figure 5 shows the relative radius difference (see caption for details). For low initial entropies, most of the composition gradient is sustained, leading to a significant radius difference between the different conductivity models. Overall, the figure indicates that the minimal primordial entropy required to destabilize most of the composition gradient depends on the planetary mass.

Fig. 6 Radius evolution for a planet with Mp = 10 M⊕ and a wide gradient. The blue, orange, and yellow lines represent the cold, warm, and hot initial entropies, respectively. Solid lines correspond to the Cond-1 model, dashed lines to the Cond-2 model, and dotted lines to the Cond-3 model.

4 Double-diffusive convection

In this work, we used the Ledoux criterion to distinguish between stable and convective regions. In stable regions, we assumed that energy transport occurs only by thermal diffusion and neglected material diffusion. In reality, regions that have been identified as stable could develop double-diffusive convection (DDC) under some circumstances. This can occur in regions where the thermal gradient is destabilizing, while the compositional gradient is stabilizing. The regimes that could develop DDC can be identified from the inverse density parameter $R_{\rho }^{-1}$ defined as (e.g., Leconte & Chabrier 2012): $$R_{\rho }^{-1}=\frac{\phi }{\delta }\frac{{\mathrm{\nabla }}_{\mu }}{{\mathrm{\nabla }}_{\mathrm{T}}-{\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}}.$$(6)

The transition between the purely diffusive energy transport regime and the DDC regime (layered or oscillatory) can be estimated using the critical inverse density parameter $R_{\rho ,\text{crit}}^{-1}$, $$R_{\rho ,\text{crit}}^{-1}=\frac{Pr+1}{Pr+\tau },$$(7)

where Pr is the Prandtl number and τ is the ratio between solute to thermal diffusivity. Convection occurs when $R_{\rho }^{-1}$ ≤ 1 (equal to Ledoux unstable), DDC occurs when $1<R_{\rho }^{-1}\le R_{\rho ,\text{crit}}^{-1}$, and stable when $R_{\rho ,\text{crit}}^{-1}<R_{\rho }^{-1}$ (see Rosenblum et al. (2011), Mirouh et al. (2012), and Leconte & Chabrier (2012) for a detailed discussion). The value of $R_{\rho ,\text{crit}}^{-1}$ is very uncertain because it depends on the material properties that determine the Prandtl number and diffusivities.

In the context of this work, the main difference in the radius evolution originates from whether a stable layer can be sustained over timescales of several gigayears. We find that hot models develop convective stair cases, while colder models retain a large stable region, and therefore the possibility of DDC may be more relevant for these cold cases. We thus explored how DDC could affect our results for the cold Cond-1 model with M_p = 10 M_⊕ and the wide composition gradient. Figure 8 presents the composition profile at two different ages of t = 1 Gyr (dark orange) and t = 10 Gyr (light orange) indicating the regions where DDC could occur. We identified the lower and upper bounds of the region that can potentially become DDC unstable. For simplicity, we adopted a constant value of $R_{\rho ,\text{crit}}^{-1}$ = 2.5 as suggested for Uranus using the same conductivity models for k_vib and k_elec (French & Nettelmann 2019). However, we note that the value of $R_{\rho ,\text{crit}}^{-1}$ is rather uncertain and would have a different value for different compositions and that its value is expected to vary with temperature and density. For dense water, values range from 1–7 (French & Nettelmann 2019). We also note that the time when DDC could start within the planetary depends on the value of $R_{\rho ,\text{crit}}^{-1}$. Higher (lower) values than considered here would lead to DDC occurring earlier (later). However, we find that varying $R_{\rho ,\text{crit}}^{-1}$ by a factor of 2 and of 0.5 barely changes the size of the zone that could become DDC unstable and therefore will not affect our conclusion (see Appendix E). In addition, we note that the self-consistent implementation of DDC in planetary evolution models is nontrivial and is still being investigated (Leconte & Chabrier 2012; Wood et al. 2013; Kurokawa & Inutsuka 2015; Fuentes et al. 2022; Anders et al. 2022; Tulekeyev et al. 2024; Dude & Hansen 2025; Fuentes 2025).

In DDC regions, both the thermal and compositional fluxes are higher relative to a fully stable configuration (Stevenson & Salpeter 1977; Rosenblum et al. 2011; Mirouh et al. 2012; Garaud 2018). In our particular example, we expect the thermal flux to be limited by the stable region. As a result, the radius evolution should be similar with and without DDC. However, we note that layered DDC (occurs at $R_{\rho }^{-1}$ close to 1) can lead to large scale convection. The exact details of DDC and its effect on planetary evolution are complex and the uncertainties originating from the initial composition and entropy profile are likely more important. Again, we note that DDC is unlikely to change the results for the warm and hot models since they are mostly convective. The fact that we considered different conductivities highlights the importance of understanding heat (and material) transport when composition gradients exist, and we therefore encourage more studies that constrain $R_{\rho ,\text{crit}}^{-1}$ for different compositions. We also hope that future work develops a more comprehensive implementation of DDC at planetary conditions.

Fig. 7 Planetary radius at t = 5 Gyr (top) and relative difference between the conductivity models (bottom). The left (Mp = 5 M⊕), middle (Mp = 10 M⊕), and right (Mp = 15 M⊕) columns show different planet masses. The circular (Cond-1), triangular (Cond-2), and square (Cond-3) markers represent the conductivity model. The light shades in the upper panel show results without convective mixing. We connected the different entropies with the same models for visibility reasons. The relative difference in the radius was calculated using ΔR/R = (R′ − RCond-1)/RCond-1, where RCond-1 is the radius for Cond-1 and R′ the radius for Cond-2 or Cond-3.

Fig. 8 Heavy-element mass fraction as a function of normalized mass for the cold Mp = 10 M⊕ planet with the wide composition gradient at t = 1 Gyr (dark orange) and t = 10 Gyr (light orange). The white (convective), gray (stable,) and orange (double diffusive) background colors indicate the stability regimes assuming $R_{\rho ,\text{crit}}^{-1}$ = 2.5.

5 Discussion

5.1 Thermal conductivity

With the Cond-1 model, we tested the case where water dominates the non-convective thermal transport. In this case, there is a region within the planet where the thermal conductivity is governed by the vibration of the nuclei rather than by radiation or electron conduction. We identified a region with ρ ≤ 0.2 g/cm³ (Figure 1, region A) where the fit of French (2019) for the vibrational conductivity becomes the dominant contribution to the thermal conductivity. Since these densities lie outside the fitted range, the extrapolated conductivity in this region is highly uncertain. It remains unknown which values of the vibrational conductivity are realistic and whether the radiative conductivity is smaller. The electron conductivity contribution from water always dominates the deep interior (Figure 1, region B). Because the thermal conductivity of pure water is likely a lower boundary of the real conductivity inside sub-Neptunes and Neptunes (Scheibe et al. 2021), it remains unclear how the conductivity would change if a mixture (e.g., rock+water) is considered.

The electron conductivity from Cassisi et al. (2007) covers most of the required temperature-density regime, but its assumption of fully ionized material likely overestimates the conductivity. For example, at Earth’s core–mantle boundary, the expected thermal conductivity ranges from k_tot = 4–11 W/m/K, with vibrational contribution dominating at k_vib = 3–10 W/m/K (Lobanov et al. 2021). This is in contrast with the Cond-2 model, where electron conductivity consistently exceeds vibrational contribution.

Pure water conductivity is a simplification and a lower limit (Scheibe et al. 2021), while fully ionized electron conductivity likely overestimates the true value. Moreover, there is evidence suggesting that vibrational (or nuclear) contributions are non-negligible. We conclude that thermal conductivity under planetary interior conditions requires further study, in particular for mixtures expected in sub-Neptunes and Neptunes.

5.2 Rosseland mean opacity tables

In Section 3.1, we identify the temperature-density regimes that are important for the planetary evolution. We identified a large portion within the planet where thermal transport is dominated by photons, but the radiative opacities do not extent to high enough densities (Figure 1, region C). Although the opacities from Freedman et al. (2014) as implemented in MESA provide higher densities (log $\tilde{R}$[g/cm³] ≤ 9), they do not properly cover temperatures above log T [K] ≥ 3.6. Additionally, Figure A.1 shows that the high-density regions have been extrapolated. The temperature-density profile of a planet roughly follows curves with constant density parameter of $\tilde{R}$. Therefore, it is a good indicator of where opacity calculations are required. All the planetary models presented in this study have values within 3 ≤ log $\tilde{R}$ [g/cm³] ≤ 8. We note that $\tilde{R}$ can have much lower values in the planet’s atmosphere. Higher values are expected for much colder and denser planets. To properly account for the radiative transport of energy inside sub-Neptunes and Neptunes, Rosseland-mean-opacity tables are required in the density parameter range at least up to log $\tilde{R}$[g/cm³] ≤ 8 for various heavy-element mass fractions.

From Figure 3, we further constrained the temperature-regime where opacity calculations are required. At higher densities of ρ ≳ 0.1 g/cm³, the temperature regime in which radiative energy transport dominates becomes smaller. In our setup with AESOPUS2.1 and the conductivity models, the highest required temperature is ~5000 K. However, at lower densities of ρ ≲ 0.1 g/cm³ (Figure 1, region D), the regime in which radiative energy transport dominates extends to temperatures above 5000 K.

Furthermore, we note that the interior structure of sub-Neptunes and Neptunes is still relatively uncertain. It is possible that such planets have more complex interiors that include composition gradients and layers of different heavy-element mass fraction (Valletta & Helled 2022; Cano Amoros et al. 2024; Morf & Helled 2025). Therefore, it is desirable to have opacity tables for various heavy-element mass fractions reaching values up to Z = 1.

5.3 Importance for exoplanets

Connecting the radius, mass, and age (time) to the planetary internal structure is very important for connecting planet formation with current-day observations. Our study clearly shows that the radius evolution of sub-Neptunes and Neptunes strongly depends on the assumed conductivity and the initial entropy profile. In particular, for sub-Neptunes and Neptunes, the conductivity plays a key role. In the case of intermediate-mass planets, the primordial composition gradient is expected to exist above a significant fraction of the mass compared to gas giant planets (Helled & Stevenson 2017; Ormel et al. 2021; Valletta & Helled 2022). Therefore, the gradient that acts as a thermal boundary layer is located above a significant amount of the internal energy budget. As a result, differences in the thermal transport (determined by the conductivity) can significantly affect the cooling. In addition, the density of the hydrogen–helium atmosphere is very sensitive to temperature. Hence, the hydrogen–helium envelope strongly correlates the thermal flux from the deep interior with the radius. Therefore, we expect that the evolution and internal structures of intermediate-mass planets are particularly sensitive to thermal conductivities.

6 Summary and conclusions

We studied the impact of the thermal conductivity on the evolution of sub-Neptunes and Neptunes. We improved on previous results from Paper I by including the effect of convective mixing. The main findings of this paper can be summarized as follows:

• The available data of thermal conductivity and radiative opacities are still insufficient for sub-Neptunes and Neptunes. Publicly available Rosseland mean opacities do not cover high enough temperatures, high enough densities, or high metallicities. The thermal conductivity of mixtures (water–rock, water–hydrogen–helium, water–methane) is uncertain;

• The thermal conductivity affects convective mixing and the final composition profile. High conductivity can inhibit convection;

• The inferred effect on the radius is more complex than in Paper I. We identified two cases: First, for low entropies convective mixing does not create a convective staircase over the entire composition gradient. A low conductivity creates a thermal boundary layer (see also Scheibe et al. (2021) and Paper I). Second, for high entropies a convective staircase replaces the composition gradient. In this case, the energy transport into the envelope is enhanced. Depending on the shape of the final composition profile, the radius can converge to similar values regardless of the conductivity model;

• The composition gradient can be too step to be eroded by convective mixing. In this case a thermal boundary layer remains, and the thermal evolution of the planet strongly depends on the assumed thermal conductivity.

Our results clearly indicate that further work is needed to better model the thermal transport inside sub-Neptunes and Neptunes. Improved calculations (and experiments) of the radiative opacities at higher densities and temperature for high metallicities are required. The determination of the thermal conductivity of various mixtures is also desirable. Finally, we note that in order to connect observations of “evolved planets” using evolution simulations, further constrains on the initial entropy and compositions are needed.

Acknowledgements

This work was supported by the Swiss National Science Foundation (SNSF) through a grant provided as a part of project number 215634: https://data.snf.ch/grants/grant/215634. We thank the referee for their helpful comments. We also thank Simon Müller and Henrik Knierim for many valuable discussions and technical support. Software: gentle_mixing (Knierim & Helled 2024; Helled et al. 2025), MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023), PyMesaReader, AESOPUS2.1 (Marigo & Aringer 2009; Marigo et al. 2024), Jupyter Notebook (Kluyver et al. 2016; Granger & Pérez 2021), NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration 2013, 2018; Astropy Collaboration 2022).

References

Appendix A Extension of the radiative opacity tables

In Paper I we created a set of Rosseland mean opacity tables for planetary evolution simulations, using the AESOPUS2.1 code (Marigo & Aringer 2009; Marigo et al. 2024). The tables span a wide range of temperature with log(T/K) = [2, 4.5] and density parameter $\tilde{R}$ = ρ/(10⁻⁶T/K)³ with log($\tilde{R}$/(g/cm³)) = [−8, 6]. Planets usually reach higher values of $\tilde{R}$ during their evolution. Here we compare two possible extensions of the tables to values of log($\tilde{R}$/(g/cm³)) = 9. The first method assumes a constant opacity for values above the table limit such that κ(T, $\tilde{R}$ ≥ 10⁶ g/cm³) = κ(T, $\tilde{R}$ = 10⁶ g/cm³) (MESA default). The second method extrapolates linearly from the last two table entries in a logarithmic space with $$\log \kappa (\log \tilde{R}\ge 6)=\log \kappa (6)+\frac{\log \kappa (6)-\log \kappa (5.8)}{0.2}(\log \tilde{R}-6),$$(A.1)

where the dependence of κ(log $\tilde{R}$) on log T is omitted for readability reasons and $\tilde{R}$ is measured in g/cm³. A linear extrapolation of κ instead of log κ would result in negative opacities. In Figure A.1 we show the scaling of the opacity with respect to the density. To compare our method of choice to extent the tables to higher densities we plot both methods in the region above log($\tilde{R}$/(g/cm³)) = 6. From these plots we conclude that the extrapolation seems to be the better choice. Keeping the opacities constant at high densities seems to significantly underestimate the opacities at high density. Especially at low temperatures, where the Freedman opacities show an increase in the opacity over multiple orders of magnitude. The Freedman opacities extend to log($\tilde{R}$/(g/cm³)) = 9 yet the tables have been extended by keeping the opacity constant at high densities. The density at which the opacity is constant is between log ρ[g/cm³] = (−2.75, −1.5) depending on the temperature.

Fig. A.1 Radiative opacity as a function of the density parameter $\tilde{R}$ (with density shown on the upper x-axis) from the AESOPUS2.1 tables. For comparison, the Freedman opacity curves as implemented in MESA are shown in light shades. The three panels, from left to right, correspond to temperatures of log T [K]= 2, 3, and 3.6. The blue, red, and yellow lines represent metallicities of Z = 0.80, 0.20, and 0.02, respectively. Solid lines indicate regions with tabulated data points (log($\tilde{R}$/(g/cm3)) ≤ 6), dotted lines show linear extrapolation on logarithmic values (see eq. A.1), and dashed lines denote constant values outside the data region.

We note that the higher density opacities are less relevant where the thermal transport is not dominated by photons. This can either be the case when electron or vibrational thermal conductivity is higher or the region is convective.

Fig. A.2 Radius evolution as a function of time for a Mp = 10 M⊕ planet, comparing capped and extrapolated radiative opacities. The blue lines represent a model with a composition gradient, while the orange lines show a model with flat composition profile. Both models assume a bulk heavy-element mass fraction of Z = 95. Lighter shades indicate the evolution using the extrapolated opacity method and darker shades indicate the evolution using the capped opacity method for out of table values.

In Figure A.2 we show the radius evolution for two different composition profiles. One model uses a flat composition profile, while the other uses the wide composition gradient. For both models we assume a bulk heavy-element mass fraction of Z = 0.95 and mass of M_p = 10M_⊕. Note, the initial central entropy is different for the two composition profiles. After 10 billion years the difference in radius is a few percent because of the different extension method of the radiative opacity tables.

Appendix B Comparison of different composition gradient widths

Fig. B.1 Initial profiles for the medium and narrow composition gradients, showing specific entropy (left) and temperature (right) as a function of normalized mass. The heavy-element mass fraction is overlaid in black in both panels, with its axis on the right. The colors blue (cold), orange (warm), and yellow (hot) correspond to different primordial entropies. The darker (Δq = 0.01) and lighter (Δq = 0.001) shades indicate different transition widths.

Fig. B.2 Radius evolution for a planet with Mp = 10 M⊕ with the medium (left) and narrow (right) composition gradient. The blue, orange, and yellow lines represent the cold, warm, and hot initial entropies, respectively. Solid lines correspond to the Cond-1 model, dashed lines to the Cond-2 model, and the dotted lines to the Cond-3 model.

Figure B.1 shows the initial entropy, heavy-element mass fraction, and temperature profiles for planets with a medium (Δq = 0.01) and narrow (Δq = 0.001) composition gradients. The narrow composition gradient is similar to the narrow composition profile presented in Paper I. Figure B.2 shows the corresponding radius evolution. We find that a thinner transition width creates a steeper composition gradient, which increases stability against convection (see eq. 3). For the medium composition gradient (Δq = 0.01), we find that no staircase appears for a cold start. For the hot start, all the assumed conductivities lead to the formation staircases and eventually, the planets cool down to have similar radii. For the narrow composition gradient (Δq = 0.001), the composition gradient is too steep to be eroded by convective mixing regardless of the initial entropy or thermal conductivity model. Therefore, the thermal evolution is the same as in Paper I. For Cond-2 the conductivity is sufficiently high to effectively transport the thermal energy from the deep interior to the outer envelop. In the cases of Cond-1 and Cond-3, the thermal energy is released over a longer period of time.

Appendix C Radius evolution for different masses

Figure C.1 shows the radius evolution for a planet with M_p = 5 M_⊕ and M_p = 15 M_⊕.

Fig. C.1 Radius evolution for planets with Mp = 5 M⊕ (left) and Mp = 15 M⊕ (right). The blue, orange, and yellow lines represent the cold, warm, and hot initial entropies, respectively. Solid lines correspond to the Cond-1 model, dashed lines to the Cond-2 model, and the dotted lines to the Cond-3 model.

Appendix D Final Composition Profiles

Fig. D.1 Heavy-element mass fraction as a function of radius at t = 10 Gyr for the planet with Mp = 5 M⊕ (first row), Mp = 10 M⊕ (second row), and Mp = 15 M⊕ (third row). The three columns show the three conductivity models Cond-1 (left), Cond-2 (middle), and Cond-3 (right). Different colors indicate the different initial entropies cold (blue), warm (orange), and hot (yellow). Solid lines represent non-convective regions and dotted lines represent convective regions.

Figure D.1 shows the final composition profile for the wide composition gradient.

Appendix E Double diffusive regions

Fig. E.1 Same as Figure 8 but with $R_{\rho ,\text{crit}}^{-1}$ = 1.25 (left), and 5 (right) for t = 1 Gyr (dark orange) and t = 10 Gyr (light orange). Note that our models do not consider the onset of DDC.

Figure E.1: same as Figure 8 but assuming different values of $R_{\rho ,\text{crit}}^{-1}$.

All Tables

All Figures

Fig. 1 Boundaries of different conductivity sources in density-temperature space. The red area (krad AESOPUS) shows the tabulated region from the AESOPUS2.1 tables (Marigo et al. 2024; Eberlein & Helled 2025), and the dark blue (kelec FR17) and light blue (kvib F19) area shows the region with ab initio data points for partially ionized water (French & Redmer 2017; French 2019). The green (kelec C07) area shows the tabulated electron conductivity as implemented in MESA (Cassisi et al. 2007). Hatched regions with squared lines (kvib dominated) or dots (kelec dominated) indicate where kvib or kelec contribute more than 50% to the total conductivity assuming the Cond-1 model with Z = 20 for the radiative conductivity. The density temperature regimes around A, B, C, and D are of particular interest and are further discussed in Section 5.

In the text

Fig. 2 Initial profiles for the wide composition gradient, showing specific entropy (top) and temperature (bottom) as a function of normalized mass. The heavy-element mass fraction is overlaid in black in both panels, with its axis on the right. Blue (cold), orange (warm), and yellow (hot) correspond to different primordial entropies. The dotted (5 M⊕), solid (10 M⊕), and dashed (15 M⊕) lines indicate different planetary masses.

In the text

Fig. 3 Thermal conductivity as a function of temperature and density (top row) and as a function of density along the initial and final planet profile (bottom row). Top: total conductivity heat map as a function of density and temperature for the three different conductivity models (from left to right): Cond-1, Cond-2, and Cond-3. The color represents the total conductivity value. Black lines show an example temperature-density profile for a planet with Mp = 10 M⊕ under the hot start scenario, calculated with the respective conductivity model: AESOPUS (Marigo et al. 2024; Eberlein & Helled 2025), FR17 (French & Redmer 2017), F19 (French 2019), and C07(Cassisi et al. 2007). The upper line (Initial profile) corresponds to the initial state of the evolution, while the lower line (Final profile) represents the state after 10 Gyr. The dotted line indicates a convective region within the planet. Hatched regions with diagonal lines (radiation dominated), squared lines (vibration dominated), and dots (electron dominated) indicate where a specific conductivity mechanism contributes more than 50% to the total conductivity. Colored lines with labels mark the boundaries of the conductivity models similarly to Figure 1. For this plot, we assumed Z = 0.20 for the radiative conductivity and Z = 1.0 for the electron conductivity of C07. We note that the distribution of heavy elements changes within the planetary interior. Bottom: conductivity as a function of density for the initial (gray) and final (black) profiles shown in the top panel. The thick solid line shows the total conductivity along the density–temperature profile of the planet, while the thin solid (radiative), dashed (vibrational), and dotted (electron) lines show the individual conductivity contributions.

In the text

	Fig. 4 Heavy-element mass fraction versus radius at different times for the Cond-1 model (top) and for the three different models at t = 5 Gyr (bottom). Solid lines indicate non-convective regions, while dashed lines indicate convective regions. This plot corresponds to the simulations of a planet with Mp = 10 M⊕, a wide composition gradient, and a warm start.
In the text

Fig. 5 Radius over time with and without mixing. The simulations with mixing are shown in dark colors, while the simulations without mixing are shown in light colors. The solid lines represent the models with a wide composition gradient, while the dotted lines represent the model with a narrow composition gradient. The blue color indicates a cold start, while the yellow color indicates a hot start. The simulations for this plot assume Mp = 10 M⊕, a wide gradient, and Cond-1.

In the text

	Fig. 6 Radius evolution for a planet with Mp = 10 M⊕ and a wide gradient. The blue, orange, and yellow lines represent the cold, warm, and hot initial entropies, respectively. Solid lines correspond to the Cond-1 model, dashed lines to the Cond-2 model, and dotted lines to the Cond-3 model.
In the text

Fig. 7 Planetary radius at t = 5 Gyr (top) and relative difference between the conductivity models (bottom). The left (Mp = 5 M⊕), middle (Mp = 10 M⊕), and right (Mp = 15 M⊕) columns show different planet masses. The circular (Cond-1), triangular (Cond-2), and square (Cond-3) markers represent the conductivity model. The light shades in the upper panel show results without convective mixing. We connected the different entropies with the same models for visibility reasons. The relative difference in the radius was calculated using ΔR/R = (R′ − RCond-1)/RCond-1, where RCond-1 is the radius for Cond-1 and R′ the radius for Cond-2 or Cond-3.

In the text

	Fig. 8 Heavy-element mass fraction as a function of normalized mass for the cold Mp = 10 M⊕ planet with the wide composition gradient at t = 1 Gyr (dark orange) and t = 10 Gyr (light orange). The white (convective), gray (stable,) and orange (double diffusive) background colors indicate the stability regimes assuming $R_{\rho ,\text{crit}}^{-1}$ = 2.5.
In the text

Fig. A.1 Radiative opacity as a function of the density parameter $\tilde{R}$ (with density shown on the upper x-axis) from the AESOPUS2.1 tables. For comparison, the Freedman opacity curves as implemented in MESA are shown in light shades. The three panels, from left to right, correspond to temperatures of log T [K]= 2, 3, and 3.6. The blue, red, and yellow lines represent metallicities of Z = 0.80, 0.20, and 0.02, respectively. Solid lines indicate regions with tabulated data points (log($\tilde{R}$/(g/cm3)) ≤ 6), dotted lines show linear extrapolation on logarithmic values (see eq. A.1), and dashed lines denote constant values outside the data region.

In the text

Fig. A.2 Radius evolution as a function of time for a Mp = 10 M⊕ planet, comparing capped and extrapolated radiative opacities. The blue lines represent a model with a composition gradient, while the orange lines show a model with flat composition profile. Both models assume a bulk heavy-element mass fraction of Z = 95. Lighter shades indicate the evolution using the extrapolated opacity method and darker shades indicate the evolution using the capped opacity method for out of table values.

In the text

Fig. B.1 Initial profiles for the medium and narrow composition gradients, showing specific entropy (left) and temperature (right) as a function of normalized mass. The heavy-element mass fraction is overlaid in black in both panels, with its axis on the right. The colors blue (cold), orange (warm), and yellow (hot) correspond to different primordial entropies. The darker (Δq = 0.01) and lighter (Δq = 0.001) shades indicate different transition widths.

In the text

	Fig. B.2 Radius evolution for a planet with Mp = 10 M⊕ with the medium (left) and narrow (right) composition gradient. The blue, orange, and yellow lines represent the cold, warm, and hot initial entropies, respectively. Solid lines correspond to the Cond-1 model, dashed lines to the Cond-2 model, and the dotted lines to the Cond-3 model.
In the text

	Fig. C.1 Radius evolution for planets with Mp = 5 M⊕ (left) and Mp = 15 M⊕ (right). The blue, orange, and yellow lines represent the cold, warm, and hot initial entropies, respectively. Solid lines correspond to the Cond-1 model, dashed lines to the Cond-2 model, and the dotted lines to the Cond-3 model.
In the text

Fig. D.1 Heavy-element mass fraction as a function of radius at t = 10 Gyr for the planet with Mp = 5 M⊕ (first row), Mp = 10 M⊕ (second row), and Mp = 15 M⊕ (third row). The three columns show the three conductivity models Cond-1 (left), Cond-2 (middle), and Cond-3 (right). Different colors indicate the different initial entropies cold (blue), warm (orange), and hot (yellow). Solid lines represent non-convective regions and dotted lines represent convective regions.

In the text

	Fig. E.1 Same as Figure 8 but with $R_{\rho ,\text{crit}}^{-1}$ = 1.25 (left), and 5 (right) for t = 1 Gyr (dark orange) and t = 10 Gyr (light orange). Note that our models do not consider the onset of DDC.
In the text