© The Authors 2025

1. Introduction

Unraveling the origin of gas giants is a fundamental goal of planetary science. Naturally, the advent of accurate atmospheric measurements from space missions such as JWST (Gardner et al. 2006) and ARIEL (Tinetti et al. 2018) raises the question of how these data can advance this effort. Planets form in disks with diverse chemical and physical conditions. It has been suggested that the planetary formation history leaves an imprint on their present-day composition (e.g., Öberg et al. 2011; Schneider & Bitsch 2021; Turrini et al. 2021; Knierim et al. 2022), although deciphering this imprint is not always straightforward (e.g., Mollière et al. 2022). Translating atmospheric abundance measurements into constraints on the bulk composition – and, by extension, formation pathways – remains a major challenge.

Traditionally, atmospheric abundances have been interpreted as the planetary bulk composition (i.e., assuming that the planet is homogeneously mixed) or as the envelope’s composition above a heavy-element core (e.g., Sing et al. 2024). Recent models have begun refining this framework; for instance, by integrating atmospheric and interior retrievals (Acuña et al. 2024; Wilkinson et al. 2024). However, the underlying assumption of a compositionally homogeneous planet or a distinct core-envelope structure remains largely unchanged. In contrast, internal structure models of Jupiter and Saturn that fit the gravity data from the Juno (e.g., Bolton et al. 2017) and Cassini (e.g., Spilker 2019) missions suggest that both gas giants deviate from this simplified picture. Instead, they have inhomogeneous interiors and fuzzy cores (see recent review by Helled & Stevenson 2024, and references therein). Building on these insights, Bloot et al. (2023) developed an interior structure retrieval code with more intricate composition profiles. Although an important step forward, this approach assumes a static composition, neglecting compositional changes through mixing. To properly interpret atmospheric measurements of gas giant exoplanets, it is essential to model the mixing (of a given primordial composition profile) and determine the composition profiles inside giant planets for various conditions.

In Knierim & Helled (2024) (hereafter KH24), we investigated the erosion of primordial composition gradients and showed that the survival of planetary cores depends on their shapes and primordial entropies. Moreover, we found that steep composition gradients with high metallicities and low primordial entropies (i.e., initially cold) can persist for gigayears, leading to atmospheric abundances that differ significantly from the planetary bulk composition. In contrast, low-metallicity profiles with high primordial entropies (i.e., initially hot) erode rapidly (∼10⁷ yr), leading to the opposite trend of well-mixed planets, whereby the atmospheric composition is the same as the bulk composition. Focusing on self-consistently reproducing present-day Jupiter and Saturn, Tejada Arevalo et al. (2025) and Sur et al. (2025) independently confirmed these findings using their evolution code APPLE (Sur et al. 2024). Our previous study laid the theoretical foundation for understanding mixing processes in giant planets. In this work, we examine the observational implications, in particular the effect on the atmospheric abundances and the planetary radius. Moreover, we show how combining these two measurements can constrain the primordial planetary conditions, providing important insights for formation models. This Letter is structured as follows. Section 2 outlines our methods. In Sect. 3, we investigate how convective mixing influences the atmospheric metallicity and radius evolution. Our discussion and conclusions are presented in Sects. 4 and 5, respectively.

2. Methods

The methods of this study are similar to the ones introduced in KH24. In short, we simulated the evolution of the planets using a modified version of the “Modules for Experiments in Stellar Astrophysics” code (MESA; Paxton et al. 2010, 2013, 2015, 2018, 2019; Jermyn et al. 2023). Our modified version uses equations of state (EoSs) that are appropriate under planetary conditions (Müller et al. 2020) and an improved treatment of convection and convective boundaries (KH24). For hydrogen and helium, we use the CD EoS (Chabrier & Debras 2021), while all elements heavier than these are represented by a 50–50 mixture of water (H₂O) and rock (SiO₂; Müller et al. 2020. All models utilize the semi-gray atmosphere model from Guillot (2010), assuming an equilibrium temperature of T_eq = 500 K (see Appendix A for details). Unlike in KH24, we assume that semi-convective regions are stable, thereby damping the erosion of the core compared to the previous study (see Sect. 4 for details). Helium rain is not included in the models shown here due to its limited effect on the planetary evolution. A further discussion on helium rain is presented in Appendix B.

We created the initial models for our numerical experiments with a proto-solar composition (Lodders 2021), before using the relax_initial_entropy and relax_initial_composition algorithms to create the desired entropy and composition profiles. We note that the (specific) entropies we refer to throughout this study are the initial values for a protosolar composition before relaxing the composition profile (see Appendix C for details). We varied these initial entropies from 8 to 11 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$, ranging from cold to hot start scenarios. The composition profiles we relaxed are the “core,” “extended,” “Helled 2023,” and “large core” models from KH24. These profiles range from shallow, extended profiles to more traditional core-envelope structures (see also Appendix C). Although these profiles are not meant to accurately represent any specific formation scenario, they qualitatively cover the expected behavior from formation studies (e.g., Helled 2023, and references therein). The planetary masses range from 0.3 M_J to 3 M_J, where M_J is Jupiter’s mass. We did not scale the composition profiles to the respective planetary mass. Instead, we simply added an envelope of protosolar composition on top of the composition profile. As a result, more massive planets have lower metallicities by construction. This setup is based on the assumption that most of the heavy elements in the planet are accreted during the early phases (phase-1, phase-2), followed by gas accretion that essentially determines the final planetary mass Helled (2023). We evolved all models for 10³ yr before initiating mixing, and then up to 10 Gyr.

3. Results

3.1. Atmospheric abundance evolution

The left panel of Fig. 1 shows the atmospheric heavy-element mass fraction versus time for different primordial entropies of the Helled 2023, extended, and core models for a 1 M_J planet. In all the simulations, the atmospheric heavy-element mass fraction increases with time. Higher primordial entropies lead to higher atmospheric heavy-element mass fraction values and faster mixing. Relative to their bulk metallicities of Z_bulk = 0.25 for the Helled 2023 model and Z_bulk = 0.18 for the extended model, the shallower composition profile of the extended model mixes more thoroughly. Indeed, for a primordial entropy of 11 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$, the extended model becomes fully mixed within 10⁷ yr, whereas the Helled 2023 model retains a core. The core model has the steepest and highest-Z composition profile and is therefore the most stable against convection. As a result, while the outer part of the gradient erodes, the inner core remains intact, leading to a plateau in the atmospheric metallicity plot.

Fig. 1. Left: Atmospheric heavy-element mass fraction vs. time for different primordial entropies and composition profiles (see legend in right plot). All models assume a 1 MJ planet. The data with errors are hypothetical atmospheric abundance measurements (see Sect. 3.3). Right: Radius vs. time for the same models as in the left plot. The age and radius of NGTS-30b and Kepler-167e are measured (see Sect. 3.3).

The evolution of Z_atm is dictated by three factors: the planetary cooling rate, the entropy profile, and the composition profile. Initially, the planet has an outer convective region and an inner radiative region (i.e., it is not adiabatic). As the planet cools, its outer convective region expands inward, thereby eroding the compositionally inhomogeneous interior. This expansion is opposed by any negative Z gradient and any positive entropy gradient. Therefore, mixing is inefficient in planets with low-entropy (“cold”) interiors and steep composition gradients. Furthermore, because hotter planets cool faster, they also mix faster (for details on convective mixing, see KH24). This interplay between the primordial entropy and the primordial composition profile leads to an increase in the planet’s atmospheric metallicity over time. Importantly, the primordial planetary structure dictates the evolution of its atmospheric abundances (see Appendix D).

Therefore, if gas giants, or a subpopulation of them, form in a similar fashion (i.e., with similar composition profiles and entropies), their origin would be reflected in their atmospheric composition over time. On the other hand, the absence of such a pattern would suggest that gas giants originate with widely varying primordial structures. Thus, if a predominant formation pathway exists, atmospheric abundance changes caused by convective mixing could provide a novel constraint on the planetary origin, particularly their primordial structure.

3.2. Mixing influences the observed radius

The right panel of Fig. 1 shows the planetary radius versus time for the simulations considered in Sect. 3.1. Initially, hotter planets exhibit larger radii, but as they cool over billions of years, their sizes decrease, converging with the ones of their cooler counterparts. Differences between primordial composition gradients, however, lead to more pronounced variations in radius. While this result is well known, our models include an often-overlooked aspect: convective mixing. Mixing significantly alters the radius evolution in two ways. First, it redistributes material, which itself leads to different densities and material functions (e.g., opacities) throughout the planet. Second, it consumes energy, thereby reducing the planet’s thermal energy content, and ultimately leading to smaller radii (see Appendix E for a detailed comparison). These processes can change the radius beyond the measured accuracy, and therefore affect the inferred planetary composition (see Fig. E.1). Therefore, convective mixing plays a crucial role in the determination of the radius evolution. Overall, our results suggest that the primordial planetary state is not always erased with time (as in adiabatic models), and can leave an observational imprint.

3.3. Mock data

Since both atmospheric metallicity and radius depend on the primordial planetary structure, combining these measurements can help infer the initial state. In this section, we illustrate this synergy by analyzing mock exoplanet data alongside available observations. For the mock data, we assumed relative errors of 15% for the atmospheric heavy-element mass fraction, 50% for the age, and 3% for the radius, based on optimistic constraints from recent JWST studies (e.g., Sing et al. 2024; Welbanks et al. 2024). For the observational data, we selected the exoplanets NGTS-30b ((0.96 ± 0.06) M_J, (0.928 ± 0.032) R_J, (1.1 ± 0.4) Gyr; Battley et al. 2024) and Kepler-167e ($1.{01}_{-0.15}^{+0.16}{\mathrm{M}}_{\mathrm{J}}$, (0.91 ± 0.04) R_J, $7.1_{-4.6}^{+4.4}\mathrm{Gyr}$; Chachan et al. 2022). These planets have masses comparable to Jupiter, equilibrium temperatures below 1000 K (i.e., are not inflated, see Sect. 4), and radii that align with our simulations. For NGTS-30b and Kepler-167e, we considered two and three hypothetical atmospheric measurements, respectively, with uncertainties matching the ones of the mock data. This combined dataset, shown in Fig. 1, forms the basis for the analysis presented below.

Planet “Mock 1” serves as our first example. With an estimated age of (2 ± 1)×10⁷ yr and an atmospheric metallicity of Z_atm = 0.075 ± 0.011, this young object aligns with core models above 9 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ and the extended model at 8 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$. These two models, however, predict notably different radii: ∼1.25 R_J for the core models above 9 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ and ∼1.12 R_J for the 8 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ extended model. This significant difference suggests that accurate radius measurements, as is shown in Fig. 1, could effectively distinguish between these competing formation pathways. In contrast, the planet “Mock 2” presents a more ambiguous case. Its radius of (1.1 ± 0.03) R_J at an age of (5.0 ± 2.5)×10⁷ yr is consistent with most models, excluding only the 8 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ Helled 2023 and larger than 9 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ core configurations. A measured atmospheric metallicity of Z_atm = 0.22 ± 0.03, however, would exclude all extended and core models, leaving only the Helled model with primordial entropies of 10 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ or higher as possible options. Similarly, planet “Mock 3” demonstrates the power of combining multiple observables. With a radius of (1.01 ± 0.03) R_J and an age of (1.50 ± 0.75)  ×  10⁸ yr, this exoplanet aligns with all Helled 2023 models. Including an atmospheric heavy-element abundance of Z_atm = 0.15 ± 0.02 would rule out the 8, 10, and 11 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ Helled 2023 models, setting limits on the primordial entropy.

Turning to the observed data, the mass, radius, and age of NGTS-30b are consistent with all Helled 2023 models and 10 to 11 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ extended models. However, the atmospheric metallicity predicted by the 8 and 9 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ Helled 2023 models is Z_atm ≲ 0.15. An accurate measurement of Z_atm, such as the ones indicated in Fig. 1, could exclude some of these scenarios.

Finally, Kepler-167e is more challenging when it comes to discriminating between the different models. Its observed radius, mass, and age align with all the models considered in this study. Nonetheless, atmospheric abundance measurements could still discern between competing theories. Specifically, these measurements could distinguish the 11 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ Helled 2023 model, the 8 and 9 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ extended models, and the core models.

These examples clearly show the importance of combining the measured mass, radius, and atmospheric metallicity with evolution models and the key role of convective mixing in exoplanetary characterization. In addition, combining these measurements can be used to exclude certain formation pathways or place upper or lower bounds on the system’s primordial entropy. Younger planetary systems are particularly promising, as both their radii and atmospheric metallicities tend to exhibit the largest differences between various primordial structures. For example, different primordial composition gradients such as the 8 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ extended model and the 10 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ core model differ by more than 0.1 R_J, making them readily distinguishable with current observational capabilities. Detecting and characterizing young exoplanets remains challenging, but recent discoveries (e.g., Barber et al. 2024) show that these systems are observable, with the number of detected young exoplanets steadily increasing. As such planets evolve, their radii tend to converge within the observational uncertainties assumed in this study. Even in such cases, atmospheric metallicity measurements remain valuable, as they can be used to exclude certain primordial structures. The effectiveness of these constraints, however, depends on both the quality of the observational data and the region of parameter space being probed. Some areas of the radius-metallicity-age space are inherently more degenerate, making it more challenging to distinguish between formation scenarios (see Sect. 4). Also, to fully interpret any individual object, a broader range of primordial entropies and composition profiles would need to be explored.

3.4. Radius evolution for different planetary masses

Figure 2 shows the radius evolution of all models for different planetary masses. Beyond the trends already discussed in Sect. 3.2, we observe that planetary radii converge with increasing planetary mass. This behavior is primarily driven by the well-known mass-independent radius of hydrogen-helium (H-He) -dominated planets (Zapolsky & Salpeter 1969). Within our model framework, more massive planets possess larger H-He envelopes, resulting in lower bulk metallicities and radii that increasingly approach the ones expected from a pure H-He planet. A more subtle effect is the reduced mixing efficiency in more massive planets at a given primordial entropy (KH24). As a result, more massive planets develop similar metallicity profiles (in absolute mass coordinates), leading to weaker mixing-induced variations in radius. On the other hand, the lower the planetary mass, the more the planetary radii differ as their bulk metallicities and their mixing properties increasingly deviate from each other. Consequently, under the assumption that more massive planets have larger H-He envelopes, the interior structures of lower-mass planets are less degenerate than the ones of higher-mass planets, allowing for stronger constraints on their primordial conditions. We note, however, that this is a strong assumption that may not always hold. It is clear that further investigations in this direction are required.

Fig. 2. Radius evolution of different planetary masses for the Helled 2023 model (top left), extended model (top right), core model (bottom left), and large core model (bottom right). The large core model for 0.3 MJ is below the shown range. The exoplanet data were obtained from the PlanetS catalog (Otegi et al. 2020; Parc et al. 2024). Data points consistent with the evolutionary tracks are highlighted with thick lines, while ones that deviate are shown with thin lines.

In order to compare our models with observations, we included all exoplanets from the PlanetS catalog (Otegi et al. 2020; Parc et al. 2024) that fall within our planetary mass range and have equilibrium temperatures below 1000 K (see Sect. 4). Exoplanets consistent with a given model’s evolutionary track (in terms of mass, age, and radius) are highlighted with thick lines. We find that although different models reproduce subsets of the observed data to varying degrees, none of them can explain all exoplanets in the sample. This is expected, given the significant spread in observed radii even for a fixed planetary mass. For example, at 1 M_J, the radii of observed exoplanets span more than 0.5 R_J, implying that there is no single formation pathway that can explain the entire population. Instead, a diversity of primordial structures is required. The current dataset remains relatively sparse, with only seven exoplanets in the PlanetS catalog that match 1 M_J and T_eq < 1000 K. As observational surveys expand, a clearer picture may emerge, potentially revealing new trends in planetary structure and evolution.

4. Discussion

The results presented here are subject to two key and distinct sources of uncertainty: (1) simplifications in the physical modeling, and (2) ambiguities in the initial conditions. On the modeling side, uncertainties arise from a simplified treatment of convection (mixing length theory), the atmosphere model, and neglected effects such as rotation (Fuentes et al. 2023). In addition, the material properties, specifically the opacities and EoSs, remain poorly constrained (e.g., Cozza et al. 2025) and can influence a planet’s evolution (Müller & Helled 2024; Howard et al. 2025). Also, the planetary metallicity was represented by an ideal mixture of water and rock, while in reality various heavy elements can exist in the interiors of giant planets (and at different depths), influencing the mixing efficiency, and therefore the long-term evolution.

Here, we have investigated four distinct primordial composition profiles and four different primordial entropy profiles. Specifically the primordial entropy, however, remains unknown and could also depend on the planetary mass. Therefore, we have also considered cases in which more massive planets have the same primordial core but higher primordial entropies than their less massive counterparts (see Appendix F). While our models only represent a small subset of all the physically plausible primordial structures, they illustrate the additional information gained by considering the atmospheric composition. To illustrate this, Fig. 3 presents a range of distinct primordial structures that are all consistent with mock observations.

Fig. 3. Primordial heavy-element mass fraction vs. mass. The panel in the top right corner shows the atmospheric metallicity vs. time for these curves, created with the analytic model from KH24 (for details, see Appendix G). The orange curves are consistent with Zatm = 0.12 ± 0.02 at (3.2 ± 1.6)  ×  108 yr (bold mock data point). The red curves are in addition consistent with with Zatm = 0.1 ± 0.02 at (3.2 ± 1.6)  ×  106 yr and Zatm = 0.15 ± 0.02 at (3.2 ± 1.6)  ×  109 yr (thin mock data points).

In addition, the ability to constrain primordial structures remains limited by the lack of observed non-inflated gas giants. The highest-quality observations typically come from close-in inflated gas giants. As long as the inflation mechanism and its magnitude remain unknown, it is not possible to infer their composition and internal structure accurately.

Furthermore, observations retrieve specific molecular abundances in the upper atmosphere (millibar to bar) while our inferred atmospheric metallicity, Z_atm, represents the heavy-element mass fraction in the outer convective envelope (1 to 10³ bar). In addition, we represent the heavy elements with an idealized mixture of water and rock. As a result, directly translating observed molecular abundances into an envelope metallicity is not straightforward and requires further investigation.

Finally, a better understanding of planet formation could provide more stringent constraints on the planetary initial compositional and thermal structure, thereby excluding some evolutionary models that, while consistent with present-day observations, are unlikely to arise from realistic formation processes.

5. Conclusions

We have explored the observational consequences of convective mixing. We have showed that accounting for convective mixing can break part of the degeneracy between hot and cold start scenarios, allowing one to constrain the primordial state of a planet with current-day observational capabilities. Our key conclusions are that:

1. The atmospheric metallicity increases with time. The shape of the Z_atm(t) curve is determined by the planet’s primordial entropy and composition profile.

2. Mixing alters the planetary radius by redistributing material.

3. Combining mass, radius, and atmospheric metallicity can break degeneracies, thereby constraining the primordial state of the planet, particularly in young planets.

4. Low-mass gas giants, for which the primordial composition profile represents a substantial fraction of the total mass, are prime candidates for having their primordial states retrieved. The higher the planetary mass (the higher the H-He mass fraction), the harder it is to discriminate among different primordial states’ formation histories.

Overall, our results indicate that the primordial planetary conditions are not always erased and can leave an observable imprint. Therefore, convective mixing in giant planets can provide an important link between the primordial planetary state and present-day observables, offering a new avenue to constrain their formation history. As we enter an era of unprecedented observational precision with JWST and soon ARIEL, we are closer than ever to delivering on their promise: unraveling the origin of giant exoplanets.

Acknowledgments

This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606.

References

Appendix A: Atmospheric boundary conditions

We adopt the default implementation of the semi-gray, globally averaged atmosphere model of Guillot (2010), as introduced in Paxton et al. (2013). Specifically, we use MESA’s atm_irradiated_opacity = ‘iterated’ option, in which the thermal opacity is iteratively adjusted to match the temperature and pressure at the base of the atmosphere. Below the atmosphere, the opacities are provided by MESA’s opacity module and depend on the local thermodynamic conditions, including composition. The visible opacity is fixed and taken from Guillot (2010), who validated this choice by comparing to more detailed atmospheric models. While a fully self-consistent treatment of atmospheric boundary conditions would be ideal, we emphasize that the mixing inside the planet is rather insensitive to the specific choice of atmospheric conditions (see KH24).

Appendix B: Helium rain

Helium rain refers to the process by which helium becomes immiscible in hydrogen under given pressure and temperature conditions, leading to phase separation between hydrogen and helium, and the formation of helium droplets that settle toward deeper layers (Stevenson & Salpeter 1977). This phase separation releases gravitational energy, which can delay the planetary cooling (Fortney & Hubbard 2003; Howard et al. 2024). Observational evidence suggests that helium rain occurs in both Jupiter and Saturn (see Helled & Stevenson 2024, and references therein), but the exact conditions for helium immiscibility remain uncertain. Several phase diagrams which differ in the extent and location of the immiscibility region have been proposed (e.g., Morales et al. 2009; Lorenzen et al. 2011; Schöttler & Redmer 2018). As a result, the impact of helium rain on the thermal and structural evolution of giant exoplanets remains unclear. To explore the effect of helium rain on our results, we conducted a series of simulations using the phase diagram of Schöttler & Redmer (2018), assuming instantaneous removal of immiscible helium. Figure B.1 compares the radius evolution of models with and without helium rain. We find that across the range of models considered, helium rain has only a small impact on the planetary evolution. This is because of two main reasons. First, the strong stellar irradiation delays the planetary cooling and therefore the onset of helium immiscibility. Second, convective mixing enriches the envelope with heavy elements, which reduces the helium concentration, and thus the amount available to separate once immiscibility begins. As a result, models with strongly enriched envelopes–such as Helled 2023 and extended, where Z_atm exceeds 10%–experience little to no helium rain. In contrast, models with more stable composition gradients, such as core and large core, retain lower envelope metallicities and do show a small increase in radius after several gigayears. Still, the effect remains small. In all cases, neither the increase in planetary radius nor the change in atmospheric metallicity caused by helium rain are currently observable. Our results imply that helium rain plays a limited role in the evolution of hot and warm Jupiters in comparison to cold giant planets.

Fig. B.1. Radius evolution of different planetary masses with and without helium rain for the core model (left) and the Helled 2023 model (right).

Appendix C: Initial model details

Figure C.1 (left) shows the heavy-element mass fraction as a function of mass for the four models investigated in this study. These models are identical to those used in KH24 (see Fig. 3 in KH24). Furthermore, Fig. C.1 (right) shows the initial entropy profiles for the numerical experiments presented in Fig. 1. These profiles are not isentropic but exhibit a nearly linear increase in specific entropy with decreasing metallicity, reflecting the imposed compositional gradient. To isolate the effects of thermal structure from those of composition, and to enable comparisons with previously published isentropic models, we define a reference entropy value based on an equivalent homogeneous planet with protosolar composition. Specifically, when we refer to an initial entropy of, for example, 8 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$, we mean that the planet would have this entropy if it were fully convective and of protosolar composition. In this sense, the quoted entropies can be interpreted as protosolar-equivalent entropies, providing a consistent benchmark across different interior structures.

Fig. C.1. Left: Initial metallicity profiles of the four models investigated in this study. Right: Entropy profiles of models Helled 2023, extended, and core investigated in this study. The entropy labels in the legend indicate the specific (homogeneous) entropy of the proto-solar composition prior to relaxing the composition gradient.

Appendix D: Analytic details

Let Z(m) be the primordial composition profile. As the outer convective region eats its way into the planet, the atmospheric metallicity increases to

$$\begin{array}{c}\hfill Z_{\mathrm{atm}}=\frac{1}{M-m_{\mathrm{RCB}}}{\displaystyle \underset{m_{\mathrm{RCB}}}{\overset{M}{\int }}}Z(m)\mathrm{d}m,\end{array}$$(D.1)

where M is the planetary mass and m_RCB is the mass coordinate of the boundary between the inner radiative and the outer convective region. As we showed in KH24, this boundary is given by

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{S}_{\mathrm{cool}}}{M-m_{\mathrm{RCB}}}=\mathrm{\Delta }s(m_{\mathrm{RCB}})+\mathrm{\Delta }{s}_{\mathrm{comp}}(m_{\mathrm{RCB}}),\end{array}$$(D.2)

where ΔS_cool is the (total) entropy loss, Δs is the (specific) entropy difference between the core and the envelope, and Δs_comp is the same difference but for the composition entropy. We define the composition entropy as the integral over the composition gradient multiplied by a response function that converts the gradient to an entropy. Using the X, Y, and Z notation for the mass fraction of hydrogen, helium, and everything else, respectively, as well as assuming a fixed ratio between X and Y (e.g., protosolar), the composition entropy simplifies to

$$\begin{array}{c}\hfill s_{\mathrm{comp}}={\displaystyle \underset{0}{\overset{M}{\int }}}{\left(\frac{\partial s}{\partial Z}\right)}_{P,\rho }\frac{\mathrm{d}Z}{\mathrm{d}m}\mathrm{d}m,\end{array}$$(D.3)

where P and ρ are the pressure and density, respectively. Thus, the function form of Δs_comp is directly influenced by the composition gradient dZ/dm. Over gigayears, planets cool to approximately the same entropy (i.e., most of the cooling happens early). Hence, in the absence of any additional mechanism that changes the primordial shape of Δs or reduced the cooling rate ΔS_cool, the evolution of the atmospheric metallicity is largely governed by the primordial entropy and the primordial composition gradient (for more details, see KH24).

Appendix E: Radius evolution details

Figure E.1 shows the radius evolution of a 1 M_J Helled 2023 model under three different scenarios: (i) with convective mixing (our default case), (ii) without mixing, and (iii) with the same bulk metallicity but with the heavy elements being homogeneously distributed. Planets with a homogeneous distribution of heavy elements converge to the same radius after a few gigayears as their initial conditions are erased. On the other hand, when mixing is suppressed, composition gradients inhibit convection, reducing the efficiency of the planetary cooling. This leads to a larger radius compared to the other models. Notably, for the 11 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ case, neglecting mixing results in an overestimate of the radius by about 0.08 R_J, which lies within current observational uncertainties. This highlights the importance of accounting for convective mixing when interpreting observational data.

Fig. E.1. Radius evolution of a 1 MJ Helled 2023 model for different entropies. The solid lines indicate models with convective mixing. Dashed lines show the same models without mixing, while dotted lines represent planets with the same bulk metallicity, but with the heavy elements being homogeneously distributed.

Appendix F: Uncertainty in primordial entropy and its mass dependence

In Sect. 3.4, we assumed a constant primordial entropy of 9 ${\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ across all planetary masses. In reality, however, the specific entropy is expected to vary with mass-particularly during runaway gas accretion, where shock heating establishes a positive entropy gradient in more massive envelopes (e.g., Cumming et al. 2018). However, the entropy gradient built up during runaway gas accretion may not affect the mixing of primordial composition gradients significantly, as mixing begins only once the outer envelope has cooled sufficiently (KH24).

The composition gradients in the deep interior correspond to the earlier formation phases and strongly depend on the accretion rates and the local formation conditions. It is yet to be determined whether the central entropy correlates with planetary mass. To explore the impact of a potential correlation between core entropy and planetary mass, we repeated the simulations of model “extended” at various masses, adjusting the entropy relative to the $9{{\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}}^{-1}$ baseline: by $-1{\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ for 0.3 M_J and 0.5 M_J, by $+1{\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ for 1.5 M_J and 2.0 M_J, and by $+2{\mathrm{k}}_{\mathrm{B}}{\mathrm{m}}_{\mathrm{u}}^{-1}$ for 3.0 M_J. Figure F.1 compares these simulations to the models used in the main text. We find that varying the entropy alters the mixing efficiency and, consequently, the atmospheric metallicity. This effect is most pronounced in the low-mass models, consistent with our findings in Sect. 3. In contrast, the radius evolution is rather insensitive to the entropy shift, with noticeable differences only for very cool, low-mass planets. These results demonstrate that while the exact outcome depends on the specific model assumptions, the overall trends remains. Therefore, it is clear that the consideration of convective mixing has an added value.

Fig. F.1. Left: Atmospheric metallicity over time for the “extended” composition model at different planetary masses. Solid curves correspond to the default entropy of 9 kB mu; dashed curves indicate models with adjusted entropies. Right: Radius evolution for the same set of simulations. See text for further details.

Appendix G: Details of the degeneracy plot

To compute all the curves shown in Fig. 3 in a reasonable time, we employed the analytic mixing model described in Appendix D. To obtain the Δs and Δs_comp, we assumed Δs(m) = α_s(s, Z)(Z(m)−Z(0)) and Δs_comp(m) = α_comp(s, Z)(Z(m)−Z(0)), where the response functions α_s and α_comp were obtained by interpolating MESA simulations. Similarly, we obtained the cooling rate ΔS_cool by interpolating MESA simulations. The metallicity curves were generated with a generalized logistic function,

$$\begin{array}{c}\hfill Z(m)=Z_{\mathrm{core}}-\frac{Z_{\mathrm{core}}-Z_{\mathrm{atm}}}{1+exp[\alpha (m-m_{\mathrm{mid}})]},\end{array}$$(G.1)

where Z_core and Z_atm are the heavy-element mass fractions in the limits of m → −∞ and m → ∞, respectively. The steepness of the transition between Z_core and Z_atm is quantified with α, the location of the transition by the midpoint m_mid. For sufficiently steep values of α, Z_core and Z_atm can be thought of as the heavy-element mass fraction at the core and in the atmosphere. To generate a large sample of initial models, we drew Z_core from a uniform distribution between 0 and 1, Z_atm from a uniform distribution between 0 and Z_core, m_mid from a uniform distribution between 0 and 1 M_J, and α from a log-uniform distribution between 1 and 200.

All Figures

Fig. 1. Left: Atmospheric heavy-element mass fraction vs. time for different primordial entropies and composition profiles (see legend in right plot). All models assume a 1 MJ planet. The data with errors are hypothetical atmospheric abundance measurements (see Sect. 3.3). Right: Radius vs. time for the same models as in the left plot. The age and radius of NGTS-30b and Kepler-167e are measured (see Sect. 3.3).

In the text

Fig. 2. Radius evolution of different planetary masses for the Helled 2023 model (top left), extended model (top right), core model (bottom left), and large core model (bottom right). The large core model for 0.3 MJ is below the shown range. The exoplanet data were obtained from the PlanetS catalog (Otegi et al. 2020; Parc et al. 2024). Data points consistent with the evolutionary tracks are highlighted with thick lines, while ones that deviate are shown with thin lines.

In the text

Fig. 3. Primordial heavy-element mass fraction vs. mass. The panel in the top right corner shows the atmospheric metallicity vs. time for these curves, created with the analytic model from KH24 (for details, see Appendix G). The orange curves are consistent with Zatm = 0.12 ± 0.02 at (3.2 ± 1.6)  ×  108 yr (bold mock data point). The red curves are in addition consistent with with Zatm = 0.1 ± 0.02 at (3.2 ± 1.6)  ×  106 yr and Zatm = 0.15 ± 0.02 at (3.2 ± 1.6)  ×  109 yr (thin mock data points).

In the text

	Fig. B.1. Radius evolution of different planetary masses with and without helium rain for the core model (left) and the Helled 2023 model (right).
In the text

	Fig. C.1. Left: Initial metallicity profiles of the four models investigated in this study. Right: Entropy profiles of models Helled 2023, extended, and core investigated in this study. The entropy labels in the legend indicate the specific (homogeneous) entropy of the proto-solar composition prior to relaxing the composition gradient.
In the text

	Fig. E.1. Radius evolution of a 1 MJ Helled 2023 model for different entropies. The solid lines indicate models with convective mixing. Dashed lines show the same models without mixing, while dotted lines represent planets with the same bulk metallicity, but with the heavy elements being homogeneously distributed.
In the text

	Fig. F.1. Left: Atmospheric metallicity over time for the “extended” composition model at different planetary masses. Solid curves correspond to the default entropy of 9 kB mu; dashed curves indicate models with adjusted entropies. Right: Radius evolution for the same set of simulations. See text for further details.
In the text