© The Authors 2025

1 Introduction

The evolution of giant planets plays a crucial role in our understanding of planetary formation and internal structure (Burrows et al. 2001; Hubbard et al. 2002; Baraffe et al. 2008; Helled et al. 2014). As we continue to discover exoplanets and deepen our insight into our Solar System’s giants, it becomes increasingly important to advance theoretical models and numerical simulations in order to interpret the rich data. The evolution of giant planets is typically modeled by solving the partial differential equations of mass, momentum, and energy conservation, as well as the transport of chemical elements and energy for example, Kippenhahn et al. (2013). Several stellar models have been adapted to model planetary evolution such as CEPAM (Guillot & Morel 1995) and Etoile (Kovetz et al. 2009; Vazan et al. 2013). Recently, a new planetary evolution code named APPLE that includes various relevant processes and parameters was presented (Sur et al. 2024). Unfortunately, these mentioned models are currently not public and unavailable to the scientific community.

On the other hand, the Modules for Experiments in Stellar Astrophysics code, i.e., MESA; Paxton et al. (2011, 2013, 2015, 2019); Jermyn et al. (2023) is open-source and widely utilized for stellar evolution simulations and simple planet models. However, enhancing its capability and relevant for giant planet modeling requires certain modifications. This paper presents several significant modifications to the MESA code (version 24.03.1). This allows the community to model the evolution of giant planets, accounting for the relevant physical processes and parameters. Our modifications include an equation of state (EoS) tailored for planetary materials, modifications of the radiative opacity, an improved treatment of convective mixing processes, and several methods of modeling helium rain or sedimentation of chemical species in general. By implementing these changes, we aim to produce more reliable models of giant planet evolution, thereby providing valuable insights into their formation and interiors.

This paper does not aim at presenting new results; instead, it provides the technical details of the modifications to the MESA code that make it a state-of-the-art planetary evolution model. The goal of this paper is to make this model accessible to the community, facilitating further research into the evolution and interiors of giant planets.

Our paper is organized as follows. In Section 2, we present the modifications to the EoS, while in Section 3 the opacity treatment. Sections 4 and 5 describe the treatments of convective mixing and helium rain, respectively. A short summary is presented in Section 7. The code is open source and is publicly available as a GitHub repository¹.

2 A planetary equation of state for MESA

A foundational aspect of modeling giant planets is the EoS, which describes how the properties of a material change under varying conditions. Since MESA was primarily developed to model the evolution of stars, the calculations rely on equations of state constructed mainly for stellar matter, which can lead to inaccuracies and limitations when applied to planets. One major issue is that for planetary conditions, the EoS tables provided with MESA are limited to heavy-element mass fractions of at most 4%. This limits the range of models that can be calculated, since the deep interiors and envelopes of giant planets can be much more metal-rich. To address this limitation, we incorporated a more suitable EoS that reflects the unique conditions within giant planets, including hydrogen-helium mixtures and other elements found in their interiors.

In this work, we present our most recent version of a planetary EoS developed to be used with MESA. Earlier versions of the EoS were already used to investigate, for example, the origin of the dilute core in Jupiter (Liu et al. 2019; Müller et al. 2020b; Meier et al. 2025), the planet’s enriched atmosphere (Müller & Helled 2024), mixing in giant planets (Knierim & Helled 2024, 2025), the characterization of giant exoplanets (e.g. Müller et al. 2020a; Müller & Helled 2021, 2023; Müller & Helled 2025), and planet formation models (Valletta & Helled 2020; Mol Lous et al. 2024; Shibata & Helled 2025).

Our EoS combines tables of single materials based on experimental data and theoretical calculations. Integrating this EoS into MESA involves (i) constructing new tables and (ii) modifying the existing EoS subroutines to ensure appropriate handling of the new tables for a large range of different compositions.

2.1 Construction

We adapted the widely used linear-mixing approximation (Amagat’s law) to construct an EoS for several components. The following relations can be readily derived by considering the Gibbs free energy of a mixture G(p, T, N) = ∑_iμ_iN_i as a thermodynamic potential, where p is the pressure, T the temperature, and μ_i are the chemical potentials. The absolute numbers of each component N_i are contained within the vector N. For an arbitrary mixture, the total density ρ of the mixture is given by $${\rho }^{-1}(p,T,\mathit{X})=\sum _i{X}_i{\rho }_i^{-1}(p,T),$$(1)

where X_i and ρ_i are the mass fraction and density of component i, the components of X are the individual mass fractions and ∑_i X_i = 1. It is further assumed that the mixture is in thermodynamic equilibrium, such that each component is at the same temperature and pressure. The total internal energy, u, and entropy, s, were calculated as the weighted sum of the individual quantities with an additional term for the entropy $$u(p,T,\mathit{X})=\sum _i{X}_i{u}_i(p,T)$$(2a) $$s(p,T,\mathit{X})=\sum _i{X}_i{s}_i(p,T)+s_{\mathrm{i}\mathrm{d}},$$(2b)

where u_i and s_i are the individual internal energies and entropies. Neglecting the free electrons, the ideal entropy of mixing s_id can be approximated as (Chabrier et al. 2019): $$s_{\mathrm{i}\mathrm{d}}=k_b(N\ln N-\sum _i{N}_i\ln N_i)=-\frac{k_b}{\overline{A}{m}_H}\sum _i{x}_i\ln x_i,$$(3)

where N = ∑_i N_i, x_i = N_i/N are the number fractions, k_b is the Boltzmann constant and m_H the atomic hydrogen mass. The mean atomic weight is $\bar{A}=\sum _i{x}_i{A}_i$, where A_i are the individual atomic mass numbers.

Following these rules, all necessary thermodynamic quantities to solve the evolution equations can be calculated, usually in more than one way. In Appendix A, we describe the calculation of all quantities required by the EoS module of MESA.

Due to its stellar-evolution heritage, MESA is set up to handle a three-component EoS given the hydrogen, helium and heavy-element mass fractions (X, Y, Z). This imposes a limitation that only one representative heavy-element component can be modeled in the interior. However, this representative does not need to be a pure substance, but can be any arbitrary mixture of several heavy elements. For our purposes, this was achieved by combining different individual heavy-element equations of state with the ideal mixing approximation.

MESA requires the logarithmic density and temperature as the thermodynamic basis, such that the EoS can be written as, for example, p = p(logρ, log T, X, Y, Z). The MESA tables use the linear combination log Q = logρ [g/cc] − 2 log T [K] + 12 in place of the density as the first independent variable. This simplifies the interpolation, since most stellar and planetary interiors tend to be close to a diagonal in the (log ρ, log T) space. Using (log ρ, log T) as the independent variables complicates the application of the ideal mixing approximation, since the individual densities ρ_i are unknown a priori. Calculating the EoS at runtime would require a root-finding procedure for each EoS call, which is numerically unfeasible. Therefore, MESA uses a set of precomputed tables for different compositions, and then interpolates between these. Our EoS adheres to this strategy; however, since it spans heavy-element fractions from zero to one, it necessitates the computation of numerous tables.

2.2 Implementation

We implemented the calculation of the EoS and the creation of the tables for MESA as a separate Python module called tinyeos². The module uses EoS tables for hydrogen, helium, and a heavy element to model arbitrary mixtures. The currently available tables are Saumon et al. (1995), Chabrier et al. (2019) and Chabrier & Debras (2021) for hydrogen and helium, and the QEOS (More et al. 1988) heavy-element tables for H₂O, SiO₂, Fe (A. Vazan, priv. comm.) and CO (Podolak et al. 2023), as well as the SESAME (Lyon 1978) and AQUA (Haldemann et al. 2020) tables for H₂O. The tables are evaluated using the bicubic spline scipy.interpolate.RectBivariateSpline, and derivatives are calculated directly from the splines.

The EoS is available with both (log ρ, log T) and (log p, log T) as the thermodynamic bases. If a table is only available for one pair of these independent variables, it is inverted by interpolating along isotherms with the piece-wise monotonic cubic spline scipy.interpolate.PchipInterpolator to create tables for the other pair. In addition to the density and temperature or pressure and temperature, the EoS uses the mass fractions of hydrogen and up to three heavy elements as inputs. While the total mixture is calculated during runtime, the three-component heavy-element part of the mixture is pre-calculated and stored as tables.

Using (log p, log T) as the thermodynamic basis significantly simplifies the mixture calculation, since each component is assumed to be at the same pressure and temperature. In our module, the default calculation of the (log ρ, log T) EoS (the required basis in MESA) performs a vectorized root-finding procedure with scipy.optimize.elementwise that uses the (log p, log T) EoS. This is numerically robust and, due to vectorisation, quite efficient³.

By default, our EoS builds the tables for MESA as follows: The EoS is calculated on rectangular grid in (log Q, log T), with log T_(min,max) [K] = 2, 6 and log Q(_min,max) = −6, 6 and a resolution of Δ log T [K] = 0.02 and Δ log Q = 0.05. The requirement that the tables be rectangular in (log Q, log T) leads to EoS calls with densities that are outside the range of the hydrogen, helium or heavy-element tables. It is therefore left to the responsible user to ensure that the planetary models are within the validity of the respective equations of state (see Figure 1).

The above procedure was repeated for all possible combinations of the elemental mass fractions with a spacing of Δ(X, Y, Z) = 0.1, and yielded 66 tables. To speed up computation, the building of the tables was parallelized using joblib.Parallel. The output of the EoS was checked for bad values, for example, infinities or nans. When that occurred, the bad values were replaced by a distance-weighted nearest-neighbour interpolation. In addition, since the calculation of spline derivatives could sometimes yield unphysical values, a few crucial quantities, for example, the adiabatic temperature gradient and pressure derivatives, were constrained to remain within certain bounds. For numerical robustness, the resulting tables were smoothed with the distance-weighted nearest-neighbour interpolation on the entire grid. The smoothing could be disabled or repeated; by default, it was done once.

These tables could then be used in MESA with our custom EoS subroutine. This subroutine contains a few crucial modifications, allowing for the full range of different compositions. Otherwise, it is largely a copy of the original MESA EoS subroutines and works in the same way. The custom EoS code and tables using the Chabrier & Debras (2021) EoS for hydrogen and helium, and a 50–50 H₂O SiO₂ mixture from the QEOS tables, is available on GitHub.

2.3 Example evolution calculations

In this section, we illustrate the regions of validity of our EoS in the temperature-density space and demonstrate the capabilities of the new EoS. We calculated two evolution models of homogeneously mixed M = 1 M_J planets with T_eq = 400 K and bulk metallicities of Z = 0.012 and Z = 0.20, and proto-solar hydrogen-helium ratios. A 50–50 water-rock mixture represented the heavy elements. We assumed hot-start initial conditions and used the irradiated-grey atmospheric model of Guillot (2010) as implemented in MESA.

The densities and temperatures for times between 1 Myr and 4.5 Gyr are shown in Figure 1 for both bulk compositions. The dashed black lines show the boundary of the EoS, and the entropy contours are shown in the background. Figure 2 shows the same models, but depicts the pressure and temperature profiles instead. Here, the contours show the difference in density Δ log ρ due to the different bulk compositions. The higher bulk metallicity leads to a significantly denser and hotter interior compared to a model with solar composition.

For modeling giant planets, the EoS presented here is a significant improvement over the current options available in MESA, as it allows for the modeling of planets with arbitrary bulk compositions. This is particularly important since giant planets are often metal-rich and do not share the same composition as their host stars (Thorngren et al. 2016; Teske et al. 2019; Müller & Helled 2023; Howard et al. 2025; Müller & Helled 2025).

3 Opacity

The thermal evolution of giant planets strongly depends on the atmospheric boundary conditions, which regulate the efficiency of energy loss to space. These boundary conditions are governed by the atmospheric temperature-pressure profile and opacities, which control radiative transfer and thus influence the cooling and contraction history of the planet (Guillot & Morel 1995; Fortney et al. 2008). Accurate treatment of opacities, particularly those of H₂, He, and trace species like H₂O, CH₄, and NH₃, is essential for modeling the radiative-convective boundary, which in turn sets the planet’s luminosity over time (Marley et al. 2007; Freedman et al. 2008). Small variations in opacity can lead to significant differences in thermal timescales, highlighting the sensitivity of evolutionary models to atmospheric assumptions.

For giant planets, opacities are also important in the deeper interior. Since giant planets may harbour composition gradients (Helled & Howard 2024), determining the modes of energy transportation requires knowledge of conductive and radiative opacities (Eberlein & Helled 2025). We note that in the deep-interior conditions relevant for giant planets, electron conduction is expected to be significantly more efficient than radiative transport by photons. In MESA, the conductive opacity is calculated from tabulated thermal conductivities. The tables are an extended version of the calculations from Cassisi et al. (2007), and were privately communicated by A.Y. Potekhin.

For the radiative contribution at high temperatures (beyond a few thousand K), the dominant opacity sources are electron scattering (Thomson scattering) and bound-free and free-free absorption by highly ionized atoms. MESA primarily uses OPAL opacities (Iglesias & Rogers 1996) for these regimes. These are extensive tables covering a wide range of temperatures, densities, and chemical compositions, computed using detailed atomic physics calculations.

In low temperatures (mostly relevant for the outer envelope and atmosphere), radiative opacities are hard to determine: At temperatures of ≲10⁴ K, the opacity becomes dominated by molecules (for example, H₂O, CO, TiO), the H⁻ ion, dust grains, or even clouds. In this regime, MESA relies on extensive pre-computed tables from specialized codes, such as those by Ferguson et al. (2005), Freedman et al. (2008, 2014), and AESOPUS (Marigo & Aringer 2009; Marigo et al. 2022, 2024).

While the low-temperature opacities of Ferguson et al. (2005) include grains, their contribution is absent from the other tables. Additionally, they do not account for potential cloud decks that could locally modify the opacity significantly. To alleviate this, we also implemented simple grain and cloud contributions that can be included in the opacity calculation and therefore mimic more complex atmospheric conditions (Poser et al. 2019).

We note that the following methods of modifying the opacities to include additional physics are clearly simplified. However, they provide an accessible way of studying how, and for which parameters, they affect the thermal evolution of giant planets. An alternative and more comprehensive approach would be to use tabulated atmospheric models to determine the boundary conditions for the evolution models. Currently, a few tables exist in MESA that can be used for stars and brown dwarfs. It would be clearly advantageous to have a variety of tabulated atmospheric models that can be used for irradiated giant planets, and we hope that they can be both constructed and implemented in the future.

Fig. 1 Time evolution of two homogeneous planets in the (log ρ, log T) space with Z = 0.012 (dash-dotted purple lines) and Z = 0.20 (solid red lines) for times between 1 Myr and 4.5 Gyr. More transparent lines correspond to earlier times. The contours show the entropy, and the dashed black lines show the boundaries of the EoS.

Fig. 2 Time evolution of two homogeneous planets in the (log p, log T) space with Z = 0.012 (dash-dotted purple lines) and Z = 0.20 (solid red lines) for times between 1 Myr and 4.5 Gyr. More transparent lines correspond to earlier times. The contours show the difference in density for the two compositions, and the dashed black lines show the boundaries of the EoS.

3.1 Grains

To model refractory grains that stay mixed in the background gas and do not settle into the clouds, we followed the methodology from Valencia et al. (2013). This simple model determines the conditions under which grains are present or absent, and provides a linear fit to the more complex calculations for a grain opacity from Alexander & Ferguson (1994). Namely, for $\log T<\log T_1^{\ast }=0.0245\log \overline{R}+3.096$, where $\log \overline{R}=\rho /{\left({10}^{-6}T\right)}^3$, the grain contribution is supposed to be maximal. This expression deviates from the one in Valencia et al. (2013). There is, in fact, a mistake in the published version, which results in an inconsistent fit to the Alexander & Ferguson (1994) grain opacity. Here, we provide the corrected version communicated to us privately by D. Valencia. For $\log T<\log T_1^{\ast }$, the opacity from the grains is $$\log {\kappa }_{\text{grains }}=0.430+1.3143(\log T-2.85),$$(4)

where T is in Kelvin, and κ_grains in cm²/g. The total opacity is simply the sum of the gas and grain opacity: κ = κ_gas + κ_grains. For $\log T>\log T_2^{\ast }=0.0245\log \overline{R}+3.221$ on the other hand, the grains are thought to have fully evaporated and κ_grains = 0. For $\log T_1^{\ast }<\log T<\log T_2^{\ast }$, the grain opacity is linearly interpolated between the two critical values ${\kappa }_{\text{grains }}(\log T_1^{\ast })$ and ${\kappa }_{\text{grains }}(\log T_2^{\ast })=0$.

The calculations from Alexander & Ferguson (1994) assumed that the grain-size distribution is that of the interstellar medium dominated by small, micro-sized grains. This results in a very high grain opacity, which is likely not very realistic for evolved giant planets (Movshovitz et al. 2010). To account for this, our implementation in MESA includes a user-defined scaling factor f_grains that is applied to Eq. (4). The MESA implementation of the grain opacity described here was already used in planet formation calculations (Valletta & Helled 2020; Mol Lous et al. 2024; Shibata & Helled 2025) and the characterization of giant exoplanets (Delamer et al. 2024).

3.2 Clouds

To account for the effect of clouds on the thermal evolution of giant planets, we implemented a simple model that treats cloud decks as an additional opacity source (Heng et al. 2012; Poser et al. 2019; Poser & Redmer 2024). The cloud opacity κ_clouds is then added to the local gas and grain opacity. The cloud opacity is parametrized as a Gaussian function: $${\kappa }_{\text{clouds }}={\kappa }_{\mathrm{c}}\cdot \exp {[-{\delta }_c(1-p/p_c)}^2],$$(5)

where δ_c and p_c are the cloud-deck thickness and location, and κ_clouds is the cloud-opacity normalisation. This yields a cloud opacity that is maximal at p_c and symmetrically spread around it with some width defined by δ_c. This formalism can easily be generalized to account for several cloud decks, in which case the total opacity is κ_clouds = ∑_i κ_{clouds, i}.

The three free parameters determining the cloud opacity are poorly constrained and depend on the molecules making up the cloud deck. Reasonable estimates for the three parameters are κ_clouds,0 ~ 0.05–1.0 cm²/g, δ_c ~ 10–100, and p_c ~ 10⁵−10⁷ Ba⁴ (Poser & Redmer 2024). For generality and simplicity, in our current implementation, the cloud-deck location is constant in time. However, we note that Poser & Redmer (2024) suggested an evolving cloud-deck location that matches pre-calculated condensation curves of the cloud species of interest (Visscher et al. 2006, 2010). As the planet cools, the condensation curves intersect the atmospheric pressure-temperature profile at different pressures, adjusting the position of the cloud deck.

3.3 Opacity windows

Recently, Müller & Helled (2024) investigated models of Jupiter in which the radiative opacity decreases significantly around T ~ 2000 K, caused by increased hydrogen transparency and potential depletion of alkali metals. This opacity window can lead to the formation of a deep radiative zone in the planetary envelope, typically between 1 and 10 kbar. The dip in the radiative opacity κ was modeled as a simple Gaussian reduction: $$\kappa ={\kappa }_0(1-\alpha {\mathrm{e}}^{-0.5\times ((\log T-\mu )/\sigma {)}^2}),$$(6)

where κ₀ is the unmodified opacity from Freedman et al. (2014), and log T is the logarithm of the local temperature in K. The parameters μ, σ and α determine the location, width and depth of the opacity reduction. The nominal values of μ = 3.3 and σ = 0.15, and α = 0.9 yield a qualitatively good match to the location and width of the opacity reduction from Guillot et al. (1994).

While detailed opacity calculations that can account for a depletion in alkali metals (Marigo et al. 2024; Siebenaler et al. 2025) exist, there are large uncertainties related to the depletion factor and individual opacities. The simple approach from Müller & Helled (2024) yields similar results, allowing for easy investigation of how different parameter values influence the evolution.

The radiative region suppresses convective mixing and thus decouples the composition of the outer atmospheric layers from the planet’s deep interior. They therefore argued that late-stage accretion of heavy elements can enrich the planetary atmosphere without significantly altering the deeper composition. It was shown that such a configuration is stable over gigayear timescales. This mechanism resolves the apparent inconsistency between observed atmospheric enrichment of Jupiter and predictions from homogeneous interior models, and has important implications for interpreting atmospheric compositions in both solar-system and extrasolar gas giants.

To demonstrate the modifications of the opacity, we calculated an evolution model of a 1 M_J planet with T_eq = 200 K using the Freedman et al. (2014) opacity as implemented in MESA. Then, we used the final model at about 4.5 Gyr to calculate the opacity profiles, including grains, clouds, or an opacity window.

The results are shown in Figure 3. The opacity changes significantly by these modifications. The increased opacity due to grains or clouds slows the cooling, yielding a different planetary radius and luminosity at a given time (Vazan et al. 2013; Poser & Redmer 2024). The opacity window does not only modify the cooling rate but also adds a deep radiative zone that could disconnect the deeper envelope from the atmosphere (Guillot et al. 1994; Howard et al. 2023; Müller & Helled 2024). All of these have implications for the characterization of planetary interiors and atmospheres, and therefore should be considered when modeling giant planets.

Fig. 3 Opacity profile of a 1 MJ planet with Teq = 200 K at about 4.5 Gyr. The pressure-temperature conditions were from the baseline model using the Freedman et al. opacity (blue line). The orange line includes grains with a scaling factor fgrains = 0.1. The effect of a single optically thick cloud deck at 10 bar with δc = 10 and κclouds,0 = 1 cm2/g is shown in the green line. The dashed purple line illustrates how an opacity window could be created due to a lack of alkali metals.

4 Convective mixing

Convective mixing also plays a vital role in giant planet evolution as it redistributes material throughout a planet, thereby directly altering its internal structure (Vazan et al. 2013; Müller et al. 2020b; Sur et al. 2024; Knierim & Helled 2025). This process requires energy (or releases it in the case of sedimentation) and modifies material properties (for example, the EoS or opacity), substantially changing the evolutionary trajectory. Accurately modeling giant planets therefore requires us to account for this fundamental process.

Since mixing plays an equally important role in stellar evolution, MESA already incorporates various convective mixing algorithms. In general, convective mixing is modeled as a diffusive process, with the diffusion coefficient determined by mixing length theory. Convective stability is evaluated using the Schwarzschild (or Ledoux) criterion. In its most basic form, this evaluation is done by computing the difference y = ∇_rad − ∇_ad (or ∇_L instead of ∇_ad in the case of the Ledoux criterion), and then locating the cell where y changes its sign. However, in the presence of discontinuous composition gradients, this “signchange” approach produces convective boundaries where ∇_rad ≠ ∇_ad on the convective side-failing to satisfy the Schwarzschild criterion (Gabriel et al. 2014). Hence, this approach can underestimate the extent of convective regions, yielding unreliable models. To address this issue, MESA first introduced the predictive_mixing algorithm, and later improved upon it with the convective_premixing algorithm. The general idea of the convective_premixing algorithm is to expand convective regions until ∇_rad = ∇_ad on the convective side of the boundary. It does this by tentatively treating adjacent radiative cells as convective, fully mixing this expanded region, and re-evaluating the Schwarzschild criterion. If the newly added cell remains radiative, the algorithm reverts to the original convective region; if it becomes convective, the expansion continues (for details, see Jermyn et al. 2023).

Although these algorithms perform well for stars, they tend to become numerically unstable when applied to giant planets. The primary reason is the more complex EoS. While most stars can be reasonably described by a fully ionized ideal gas, the planetary EoS is significantly more intricate (see Sect. 2). Under planetary conditions, small changes in composition can induce substantial variations in the overall structure. Moreover, the derivatives of the EoS often behave poorly, particularly near phase transitions. These issues result in ill-conditioned Jacobians, which can cause the Henyey solver to fail.

To address these instabilities, Knierim & Helled (2024) (hereafter KH24) introduced the gentle_mixing algorithm. This method monitors changes in the composition profile between steps. If the change exceeds a certain threshold, the algorithm reduces the mixing efficiency and shortens the timestep.

At its core, the gentle_mixing algorithm monitors changes in a cell’s composition and the model’s overall heterogeneity: $$h_i^2={\int }_0^M{(X_i^{\mathrm{\prime }}-X_i)}^2\frac{\mathrm{d}m}{M},$$(7)

where X_i and $X_i^{\mathrm{\prime }}$ denote the mass fractions of species i before and after a step, respectively. If either the abundance change or the heterogeneity h_i exceeds a user-defined threshold, MESA’s mix_factor and the time step are reduced, and a retry is initiated. The algorithm provides several options to customise this behavior (see GitHub repository).

Most notably, it includes an adaptive damping mechanism that tracks the number of retries and re-dos performed by the solver. If a retry was required in the previous step, the algorithm increases mix_factor gradually from its last (damped) value rather than resetting to the original value. If subsequent steps succeed, mix_factor is incrementally restored to its reference value. Conversely, if retries persist and mix_factor drops below a user-defined minimum, mixing is temporarily disabled to facilitate solver convergence. Similarly, the time step is not reduced indefinitely but only down to a specified lower limit. The gentle_mixing algorithm also allows for asymmetric h_i thresholds, distinguishing between high- and low-abundance regions. For example, high-Z cells, typically less stable under planetary conditions, can be treated with more stringent thresholds.

Many of these modifications were also applied to the convective_premixing algorithm. A key difference is that h_i can be evaluated during convective_premixing’s tentative region expansion, allowing the process to be aborted early if the heterogeneity becomes too large. To support these new capabilities, the do_eos_for_cell routine in star/private/micro.f90 was modified, and a new routine, solve_eos_given_PgasS, was added to eos/eos_support.f90. We note that the damping of the mixing is acceptable as long as the mixing timescale remains much shorter than the Kelvin–Helmholtz timescale—that is, as long as mixing continues to occur significantly faster than the planet’s thermal evolution. Since the convective mixing timescale in gas giants is typically six to seven orders of magnitude shorter than the Kelvin–Helmholtz timescale, we can afford to damp the mixing somewhat without affecting the overall evolution. Figure 4 illustrates this process. As a result, MESA mixes unstable regions over several smaller time steps instead of a single, large one. In addition, gentle_mixing also monitors changes during convective_premixing: If the newly expanded region alters the composition profile beyond the allowed threshold, gentle_mixing reverts the expansion and exits convective_premixing —again accompanied by a reduction in the time step.

Furthermore, we extended convective_premixing to better resolve convective boundaries under planetary conditions. In its default implementation, convective_premixing holds either pressure and temperature or density and temperature constant during the tentative expansion of a convective region. However, this can introduce significant errors, as temperature may vary strongly when mixing material in planetary interiors. To address this, we introduced two new modes that instead hold either pressure and density or pressure and total entropy in the convective region constant. While fixing any such pair of variables is, strictly speaking, unphysical, the pres-sure-entropy mode in particular provides substantially greater numerical stability than the default approach. In comparison to KH 24, gentle_mixing has undergone several improvements. Most notably, it now includes an adaptive scheme that attempts to mix as much as possible without failures. Figure 5 shows the evolution of two different initial composition profiles utilizing gentle_mixing. The extended profile (left panel) erodes significantly over 10 Gyr, leading to substantial enrichment of the envelope. On the other hand, the core-like profile (right panel) is more resilient, losing only its outer gradient during the same period. Despite the differing degrees of mixing, both models exhibit significant changes to their internal structures. Figure B.1 compares the resulting temperature profiles for simulations with and without mixing after 10 Gyr. If mixing is ignored, composition gradients artificially suppress convection, resulting in interiors that are considerably hotter than they should be. When mixing is included, only the innermost physically stable layers retain elevated temperatures; the outer interior approximately recovers the temperature structure of a homogeneous model with the same bulk composition. Thus, although the extent of convective mixing depends on the details of the entropy and composition profiles (see KH 24 for a full theoretical treatment), accurately modeling the mixing process is crucial to obtain realistic internal structures.

Fig. 4 Sketch of the gentle_mixing algorithm. Rather than transitioning directly from the initial to the final heavy-element profile, the algorithm inserts several intermediate steps. The inset in the top right corner shows how the damping of mixing efficiency is accompanied by a reduction in time step.

Fig. 5 Evolution of the composition profile for a 1 MJ planet and an entropy of 8 kB mu−1, starting from an extended composition profile (left) and a core-like structure (right).

5 Helium rain

While hydrogen and helium can be mixed, under some pressures and temperatures, helium becomes immiscible with hydrogen, leading to phase separation. This process is known as “helium rain” (Stevenson & Salpeter 1977). The settling of helium leads to energy release and can affect planetary evolution. Although the H–He phase diagram is still being investigated and the exact conditions for H–He demixing remain uncertain (Lorenzen et al. 2011; Schöttler & Redmer 2018; Brygoo et al. 2021), it is clear that helium rain has occurred in both Jupiter and Saturn (for example, Fortney & Hubbard 2004; Howard et al. 2024, and references therein). Therefore, we also incorporated helium rain and element sedimentation more generally in MESA (see Knierim & Helled 2025, for details). The implementation proceeds in three steps: (1) determine the immiscibility region using a user-supplied phase diagram, (2) remove the sedimenting species, and (3) repeat steps (1) and (2) until convergence is reached.

In step (1), the user-supplied phase diagram is read into a custom phase_diagram type. This type stores the lower and upper miscible number fractions, together with the corresponding temperature and pressure values. The user can specify the elements to be included and to apply a temperature shift. To determine whether a cell undergoes precipitation, we first linearly interpolate the lower and upper miscible fractions in pressure–temperature space. We then evaluate whether the local abundance ratio lies above, within, or below the miscibility gap. For ratios inside the gap, we calculate the depletion required for the composition to reach the lower miscibility boundary. The resulting abundance change is converted into a mass abundance difference and stored in the array dXi. We compute the changes in the other mass fractions assuming that their ratios remain constant throughout the sedimentation process. We note that in fact it is not trivial to determine the other mass fractions throughout the evolution (see Appendix C).

At this stage, dXi can be further modified by user-defined constraints. These include limiting the maximum abundance change per step, enabling many small adjustments rather than a single large rain-out event. Additional options allow one to close gaps within a precipitating region and enforce monotonicity of the composition gradient. The algorithm also accounts for core and ocean formation when the innermost cell undergoes precipitation. In this case, since there is no lower cell to transfer material to, dXi is adjusted to enrich the bottom cell up to the upper miscibility boundary. The process continues outward until all immiscible cells are resolved.

In step (2), the specified dXi values are subtracted from the sedimenting species, using one of three modes, listed here in order of increasing complexity. In instant_rain, the code evaluates the stability of the entire model and transfers the sedimenting element in a single step to the deepest miscible layer. The thermal structure is then updated globally. The bottom_up_rain mode proceeds cell by cell, starting from the deepest layer and moving upward, transferring material and updating the thermal structure incrementally. Finally, the advection_diffusion_rain mode solves the advection-diffusion equation, modeling element sedimentation as an advective process and convection as diffusion, with MESA’s D_mix used as the diffusion coefficient. Between MESA time steps, the solver iterates until a convergent solution is reached. This mode introduces two additional free parameters: (1) the sedimentation velocity of the element and (2) the time step used by the advection-diffusion solver.

In advection_diffusion_rain, the algorithm solves the standard advection-diffusion equation for an isotropic sphere: $$\frac{\mathrm{\partial }Y}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }m}\left[4\pi r^2\rho (-D_{\mathrm{mix}}\frac{\mathrm{\partial }Y}{\mathrm{\partial }r}+v_{\mathrm{s}\mathrm{e}\mathrm{d}}Y)\right]=0,$$(8)

where Y is the helium mass fraction, t, m and r the time, mass, and radial coordinates, D_mix is the convection diffusion coefficient from mixing-length theory, and v_sed is the sedimentation velocity. In the case of helium rain, the sedimentation velocity is only nonzero if helium is actively raining out according to the phase diagram. The equation is solved with the implicit backward-Euler scheme on a non-uniform grid. The method contains two additional free parameters: (1) the time step used by the solver and (2) the velocity of the sediment species. For all methods, the energy of the sedimentation is accounted for in the energy equation and MESA’s output quantities (for example, the total energy). Figure 6 compares the evolution of a homogeneous Saturn-mass planet across the different helium rain modes. The instant_rain and bottom_up_rain methods yield very similar radius evolutions and helium distributions. By contrast, the advection-diffusion solver produces a smoother core-envelope transition, shaped by the ratio of advection to diffusion. A higher sedimentation velocity results in a steeper gradient and transports more helium into the deep interior. In the case shown, we adopt a sedimentation velocity v_sed = 0.825 cm s⁻¹, though both higher and lower values are plausible given the current uncertainty in the microphysics of helium rain.

These modes provide users with significant flexibility. In light of the uncertainties surrounding element sedimentation, the implementation also supports options such as applying a temperature shift to the phase diagram to better match specific observational constraints.

The process of Helium rain is largely implemented in a new file, element_sedimentation.f90, located in the star/private directory. The core routine, do_element_sedimentation, is called during the evolve loop, immediately before MESA’s diffusion solver.

6 Comparison to other evolution models

Currently, MESA is the only open source Henyey code that is capable of simulating the full thermal and compositional evolution of giant planets. This makes direct comparisons with other studies rather challenging. Not only are alternative codes inaccessible to us, but they also often use different EoSs (and the numerical details of their implementation) and atmospheric models.

Under these constraints, we first tackled the seemingly simple case of a homogeneously mixed typical warm Jupiter with M = 1 M_J at T_eq = 500 K, and with zero or 5% bulk metallicity. However, as discussed below, even this comparison is not trivial.

We compared the following models:

• Two models calculated with MESA using the T − τ and irradiated grey (Guillot 2010) atmospheric models with the Freedman et al. (2014) gas opacities, the Chabrier & Debras (2021) hydrogen-helium EoS, and the QEOS H₂O for the heavy elements.

• One model calculated with the open-source code GASTLI (Acuña et al. 2024), which uses its own grid of atmospheric models, the hydrogen-helium equation of state from Chabrier et al. (2019) combined with the one from Howard & Guillot (2023), and the AQUA equation of state (Haldemann et al. 2020) for the heavy elements. We note that GASTLI is not a full Henyey code and is incapable of treating transport of energy or chemical elements. It calculates static models for different internal temperatures, and then derives the evolution from that.

• One model calculated with the closed-source code CEPAM (Guillot & Morel 1995) shared privately with us by S. Howard. The models used the non-grey atmospheric model from Parmentier et al. (2015), the hydrogen-helium equation of state from Chabrier et al. (2019) combined with the one from Howard & Guillot (2023), and the SESAME water equation of state (Lyon 1978) for the heavy elements. While CEPAM is a Henyey code, it currently lacks the ability to model the transport of chemical elements.

Figure 7 shows the radius evolution for both zero and 5% metallicity. We find that after a few Gyr, for Z = 0 all models agree on the planetary radius within about 2%, while for the Z = 0.05 case the results diverge more significantly, with the smallest and largest radii differing by about 5%. For the pure hydrogen-helium case, the EoSs used in the models are similar. Therefore, the observed variation is probably due to the different atmospheric models. For the Z = 0.05 case, there is an additional complication that all models use a different EoS for the heavy elements. Compared to many recently reported observational errors on the radii of giant planets, these differences between the models are rather large and would also translate into differences in the inferred bulk compositions. Given the expected accuracy of radius measurements (on the order of a few percent), theoretical uncertainties are important and should be considered (Müller et al. 2020a; Howard et al. 2025).

We also note that the evolution of the planetary radius is typically shown in units of Jupiter’s radius. However, not all models use the same radius for Jupiter: While planetary scientists commonly use the volumetric mean radius of Jupiter (6.9911 × 10⁹ cm), in exoplanetary sciences Jupiter’s equatorial radius (7.1492 × 10⁹ cm) is often used. This leads to a non-negligible difference that can affect results and conclusions drawn from models. For this comparison, for example, we had to renormalise the results from GASTLI, since it uses Jupiter’s equatorial radius.

We next present a rough comparison of the mixing and helium rain modules. For that, we reconstructed the Jupiter model presented by Sur et al. (2025) in MESA. Specifically, we adopted a pure-water EoS for the heavy elements and matched the initial temperature and composition profile. For helium rain, we used the instant_rain algorithm together with the phase diagram of Schöttler & Redmer (2018), shifted by 370 K. In addition, we compare our model to that presented in Howard et al. (2024), where Jupiter’s interior was modeled with a heavy-element core surrounded by a H–He envelope. The results are presented in Figure 8.

We find that all three models reproduce Jupiter’s present-day helium abundance. This agreement is expected, since each study applied a temperature shift to its phase diagram to achieve consistency with the measured helium abundance in Jupiter’s atmosphere. Although the detailed numerical values differ slightly, as expected given the varying assumptions, the qualitative trends are similar.

In addition, the degree of mixing in the dilute core region is similar in our model and the one of Sur et al. (2025). However, our simulation produces a smooth profile of the helium fraction in the deep interior while simultaneously maintaining a sharp boundary with the envelope. Also, our model does not include convective “staircases”, namely, convective regions of uniform composition separated by sharp composition gradients. In the absence of semi-convection, such staircases are not expected to form and are a numerical artifact (see Vazan et al. 2018, for discussion). Crucially, our model using gentle_mixing erodes much more of the outer dilute core than the model by Sur et al. (2025), which employs the sign-change approach to determine convective boundaries (see discussion in Sect. 4). In general, we find that the different codes produce similar interior profiles, but the exact shapes of the curves depend on the numerical details.

Overall, our results demonstrate that the different evolution models yield similar results, but also that a direct comparison is extremely difficult given all the details and subtleties that go into these models. It is clearly desirable to have detailed model comparisons and benchmarks, but this is beyond the scope of this paper. We note, however, that there is currently an effort underway within the consortium of the PLATO mission (Rauer et al. 2025) under work-package WP116100 to compare (and in the longer future, benchmark) planetary evolution models.

Fig. 6 Comparison of different helium rain implementations, showing the radius evolution (top) and the final helium mass fraction profile (bottom) of a homogeneous Saturn-mass model.

Fig. 7 Radius evolution for a 1 MJ planet with Teq = 500 K for zero (top) and 5% bulk metallicity (bottom). The lines correspond to different models: orange and red for MESA with theT − τ and irradiated grey atmospheric models, green for GASTLI, and blue for CEPAM.

Fig. 8 Distribution of helium (mass fraction) in present-day Jupiter. Shown are the results inferred by this work, Sur et al. (2025), and Howard et al. (2024).

7 Summary

We have introduced several modifications to the MESA code relevant for giant planets. These include a new equation of state, various treatments of the radiative opacity, and the processes of convective mixing and settling. As the MESA code is public, our modifications can be used by the community, paving the way for deeper insights into the formation and characteristics of giant planets in the Solar System and those orbiting around other stars.

Clearly, planetary evolution models can be further improved. Such improvements include more realistic atmospheric models, conductivities, and other thermodynamic properties. In addition, it would be valuable to benchmark the different evolution models in order to validate them and investigate how the simulation approach and numerical details affect the results.

We hope to address these topics in future research. Indeed, to take full advantage of the current and upcoming observational data, there is a need to enhance the capabilities for planetary evolution studies and push our understanding of these planetary objects to the next level.

Data availability

Our planetary MESA version can be found in the GitHub repository https://github.com/uzhplanets/mespa.

Acknowledgements

We thank the anonymous referee and the editor for valuable comments. We thank Andrew Cumming for guidance and valuable discussions, Sho Shibata and Mark Eberlein for contributions to the custom equation of state subroutines, and Allona Vazan for sharing the QEOS tables of H₂O, SiO₂ and Fe. We acknowledge support from the Swiss National Science Foundation (SNSF) grant 200020_215634 and the National Centre for Competence in Research ‘PlanetS’ supported by SNSF. Extensive use was also made of the Python packages NumPy (Harris et al. 2020), SciPy (Wes McKinney 2010; Virtanen et al. 2020), Matplotlib (Hunter 2007), Jupyter (Kluyver et al. 2016), fortranformat (Arnold 2025), joblib (Joblib Development Team 2020), mesa_reader (Wolf 2025), tinyeos (Müller 2025), mesatools (Müller 2024), py_mesa_helper (Knierim 2024b), and the Fortran code mesa_custom_eos (Knierim 2024a) based on the MESA EoS subroutines.

References

Appendix A Equation of state details

In Table A.1 we list the EoS quantities required by MESA. In Section 2, we explained how the density, pressure, entropy and internal energy are determined. Here, we also describe how the other quantities are calculated when using the recommended densitytemperature EoS, which employs the pressure-temperature tables and performs a root-finding procedure. In the following description, the quantities of the individual components are marked by the index i, and the partial derivatives involved are calculated directly from the bicubic splines – no finite differences are used.

With the aid of the inverse-function and triple-product rules, the pressure derivatives are calculated as $${\chi }_{\rho }={(\mathrm{\partial }\ln \rho /{\mathrm{\partial }\ln p|}_T)}^{-1},$$(A.1a) $${\chi }_T=-\frac{\mathrm{\partial }\ln \rho /{\mathrm{\partial }\ln T|}_p}{\mathrm{\partial }\ln \rho /{\mathrm{\partial }\ln p|}_T}.$$(A.1b)

The partial derivatives of the density with respect to pressure and temperature are $${\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }\ln p}|}_T=\rho \sum _i\frac{X_i}{{\rho }_i{\chi }_{\rho ,i}},$$(A.2a) $${\frac{\mathrm{\partial }\ln \rho }{\mathrm{\partial }\ln T}|}_p=-\rho \sum _i\frac{X_i{\chi }_{T,i}}{{\rho }_i{\chi }_{\rho ,i}}.$$(A.2b)

The adiabatic temperature gradient for the mixture is calculated using the triple-product rule as $${\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}=-\frac{\mathrm{\partial }\ln s/{\mathrm{\partial }\ln p|}_T}{\mathrm{\partial }\ln s/{\mathrm{\partial }\ln T|}_p},$$(A.3)

where the partial derivatives of the entropy with respect to pressure and temperature are $${\frac{\mathrm{\partial }\ln s}{\mathrm{\partial }\ln p}|}_T={s^{-1}\sum _i{X}_i{s}_i\frac{\mathrm{\partial }\ln s_i}{\mathrm{\partial }\ln p}|}_T,$$(A.4a) $${\frac{\mathrm{\partial }\ln s}{\mathrm{\partial }\ln T}|}_p={s^{-1}\sum _i{X}_i{s}_i\frac{\mathrm{\partial }\ln s}{\mathrm{\partial }\ln T}|}_p.$$(A.4b)

The other required quantities related to derivatives are calculated using the two pressure derivatives and the adiabatic temperature gradient (by employing the triple-product rule and Maxwell relations; see, for example, Hansen & Kawaler (1994)). The specific heat capacities at constant pressure and volume are $$c_p=\frac{p{\chi }_T}{\rho T{\chi }_{\rho }{\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}},$$(A.5a) $$c_V=\frac{c_p{\chi }_{\rho }}{{\mathrm{\Gamma }}_1}.$$(A.5b)

The adiabatic indices Γ₁ and Γ₃ are calculated as $${\mathrm{\Gamma }}_1=\frac{{\chi }_{\rho }}{1-{\chi }_T{\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}},$$(A.6a) $${\mathrm{\Gamma }}_3=1+{\mathrm{\Gamma }}_1{\mathrm{\nabla }}_{\mathrm{a}\mathrm{d}}.$$(A.6b)

The three remaining partial derivatives are $${\frac{\mathrm{\partial }s}{\mathrm{\partial }\rho }|}_T=\frac{c_V}{T},$$(A.7a) $${\frac{\mathrm{\partial }s}{\mathrm{\partial }T}|}_{\rho }=-\frac{p{\chi }_T}{T{\rho }^2},$$(A.7b) $${\frac{\mathrm{\partial }u}{\mathrm{\partial }\rho }|}_T=(1-{\chi }_T)\frac{p}{{\rho }^2}.$$(A.7c)

While the following three quantities are generally not relevant for standard giant planet evolution models, for completeness, we attempt to approximate them as well as possible given the current tabulated equations of state.

Calculating the mean molecular weight of the mixture requires knowledge of the number fractions of the molecules and ions involved in the mixture. Unfortunately, neither the Chabrier & Debras (2021) hydrogen-helium EoS nor the heavy-element tables include this information. However, these are tabulated in the Saumon et al. (1995) EoS, which we use to calculate the mean molecular weight of hydrogen (μ_X), helium (μ_Y) and their mixture (μ_X,Y) $$\frac{1}{{\mu }_{\mathrm{H}}}=\frac{2X}{1+x_{\mathrm{H}}++3x_{{\mathrm{H}}_2}},$$(A.8a) $$\frac{1}{{\mu }_{\mathrm{H}\mathrm{e}}}=\frac{3Y}{4(1+2x_{\mathrm{H}\mathrm{e}}++3x_{{\mathrm{H}\mathrm{e}}^{+}})},$$(A.8b) $$\frac{1}{{\mu }_{\mathrm{H},\mathrm{H}\mathrm{e}}}=\frac{1}{{\mu }_{\mathrm{H}}}+\frac{1}{{\mu }_{\mathrm{H}\mathrm{e}}},$$(A.8c)

where the x_i are the number fractions of the atoms, molecules or ions. For the heavy elements, we approximate their mean molecular weight using the atomic masses m_a of the involved chemical species, such that the total mean-molecular weight is, $$\frac{1}{\mu }=\frac{Z}{m_a}+\frac{1-Z}{{\mu }_{\mathrm{H},\mathrm{H}\mathrm{e}}}.$$(A.9)

Similarly, ignoring the heavy-element contribution and again using the Saumon et al. (1995) hydrogen-helium EoS, the mean number of free electrons per nucleon is, $$\frac{1}{{\mu }_{\mathrm{e}}}=\frac{x_{\mathrm{e},\mathrm{H}}}{{\mu }_{\mathrm{H}}}+\frac{x_{\mathrm{e},\mathrm{H}\mathrm{e}}}{{\mu }_{\mathrm{H}\mathrm{e}}},$$(A.10)

where the x_e,i are the number fractions of free electrons contributed by hydrogen or helium.

Finally, the inverse chemical potential of the free electrons η (degeneracy parameter) is calculated from the mean number of free electrons per nucleon and with the rational function approximation of the inverse Fermi-Dirac integral from Antia (1993).

Appendix B Effect of mixing on the inferred temperature profile

The assumed internal structure and whether convective mixing is included significantly affect the temperature profile within the planet. Figure B.1 compares the resulting temperature profiles for simulations with and without mixing after 10 Gyr. If mixing is ignored, composition gradients artificially suppress convection, resulting in interiors that are considerably hotter than they should be.

Fig. B.1 Comparison of the inferred temperature profiles after 10 Gyr for models with mixing, without mixing, and for a homogeneous profile with the same bulk metallicity as the corresponding composition gradient. All the models were run under the same assumptions as in Fig. 5.

Appendix C Atmospheric helium mass fraction with helium rain

As helium rains and helium droplets settle down to deeper layers, the planetary bulk composition has to be conserved. However, there is more than one way to ensure constant composition and it is unclear which treatment is most appropriate.

For example, let’s assume that the primordial composition of the atmosphere is X = 0.7, Y = 0.28, Z = 0.02, and at present, Y = 0.1. To keep X + Y + Z = 1, one can choose between the following options:

1. Assume that (X/Z)_primordial = (X/Z)_new and with having Y_new = 0.1, one can find X_new and Z_new.

2. If Y_new = 0.1, we can assume that 18% of helium that rained down is replaced by hydrogen, so that X_new = 0.88 and Z_primordial = Z_new = 0.02. In this case, however, the local hydrogen-to-metals ratio changes.

3. Assume that the missing 18% are of solar composition. In this case, however, “new” helium is brought to the atmosphere.

In our simulations with the instant_rain algorithm, we used the first option corresponding to a rain of pure helium. This preserves the ratio of hydrogen to heavy elements since it is the same in both the material from which the rain is formed and the material replacing the droplets (by convection from below). In that case, if we start with X = 0.7, Y = 0.28, and Z = 0.02 the atmospheric composition at present-day would be X = 0.875, Y = 0.1, and Z = 0.025.

We note that each option would lead to different results. In addition, reality is surely more complicated: Raindrops could also take down some hydrogen and heavy elements (for example, neon; see Stevenson & Salpeter 1977). In addition, a slight enrichment of heavy elements in the envelope may be expected as most of the species of interest (such as water) prefer metallic hydrogen over helium.

We suggest that future research should be dedicated to modeling the planetary evolution accounting for elemental mixing and phase separation self-consistently. This would require conserving the total mass of each species individually while ensuring that at every timestep, the phase diagram (which has multiple dimensions) is satisfied, determining X, Y and Z of the material being left behind as raindrops form and material is convected up from below.

All Tables

All Figures

	Fig. 1 Time evolution of two homogeneous planets in the (log ρ, log T) space with Z = 0.012 (dash-dotted purple lines) and Z = 0.20 (solid red lines) for times between 1 Myr and 4.5 Gyr. More transparent lines correspond to earlier times. The contours show the entropy, and the dashed black lines show the boundaries of the EoS.
In the text

	Fig. 2 Time evolution of two homogeneous planets in the (log p, log T) space with Z = 0.012 (dash-dotted purple lines) and Z = 0.20 (solid red lines) for times between 1 Myr and 4.5 Gyr. More transparent lines correspond to earlier times. The contours show the difference in density for the two compositions, and the dashed black lines show the boundaries of the EoS.
In the text

Fig. 3 Opacity profile of a 1 MJ planet with Teq = 200 K at about 4.5 Gyr. The pressure-temperature conditions were from the baseline model using the Freedman et al. opacity (blue line). The orange line includes grains with a scaling factor fgrains = 0.1. The effect of a single optically thick cloud deck at 10 bar with δc = 10 and κclouds,0 = 1 cm2/g is shown in the green line. The dashed purple line illustrates how an opacity window could be created due to a lack of alkali metals.

In the text

	Fig. 4 Sketch of the gentle_mixing algorithm. Rather than transitioning directly from the initial to the final heavy-element profile, the algorithm inserts several intermediate steps. The inset in the top right corner shows how the damping of mixing efficiency is accompanied by a reduction in time step.
In the text

	Fig. 5 Evolution of the composition profile for a 1 MJ planet and an entropy of 8 kB mu−1, starting from an extended composition profile (left) and a core-like structure (right).
In the text

	Fig. 6 Comparison of different helium rain implementations, showing the radius evolution (top) and the final helium mass fraction profile (bottom) of a homogeneous Saturn-mass model.
In the text

	Fig. 7 Radius evolution for a 1 MJ planet with Teq = 500 K for zero (top) and 5% bulk metallicity (bottom). The lines correspond to different models: orange and red for MESA with theT − τ and irradiated grey atmospheric models, green for GASTLI, and blue for CEPAM.
In the text

	Fig. 8 Distribution of helium (mass fraction) in present-day Jupiter. Shown are the results inferred by this work, Sur et al. (2025), and Howard et al. (2024).
In the text

	Fig. B.1 Comparison of the inferred temperature profiles after 10 Gyr for models with mixing, without mixing, and for a homogeneous profile with the same bulk metallicity as the corresponding composition gradient. All the models were run under the same assumptions as in Fig. 5.
In the text