© The Authors 2024

1 Introduction

Ground-based infrared observations from Low (1966) revealed that Jupiter and Saturn emit more energy than they receive from the Sun. Their observed luminosities cannot be explained solely by the continued cooling and contraction of these planets and it was suggested that both planets possess an additional internal heat source. The origin of this heat source was proposed to be the differentiation of helium, due to a phase separation between hydrogen and helium (H-He; e.g., Smoluchowski 1967; Salpeter 1973). As helium becomes immiscible with metallic hydrogen (Stevenson 1975; Stevenson & Salpeter 1977b) present in the interior of giant planets, helium droplets can form and sink (known as helium rain), providing an additional source of thermal energy. Homogeneous evolution models which do not include change of composition with time, have consistently found that Jupiter’s cooling time is roughly in line with the age of the Solar System (Hubbard 1969). With the age of the Solar System being 4.56 ± 0.10 Gyr (Connelly et al. 2012), the inferred cooling times of 4.2 to 5.3 Gyr for Jupiter (Graboske et al. 1975; Hubbard 1977; Saumon et al. 1992; Guillot et al. 1995) indicate that helium de-mixing may not have started at all or started only recently (Stevenson & Salpeter 1977a). On the contrary, similar models for Saturn predict significantly shorter cooling times, between 2 and 2.7 Gyr (Pollack et al. 1977; Saumon et al. 1992; Fortney & Nettelmann 2010), suggesting a strong excess of luminosity. This, in turn, suggests that in Saturn helium de-mixing is significant.

Nevertheless, the excess luminosity observed for Jupiter and especially Saturn originally motivated the idea of helium de-mixing; however, this does not serve as undeniable evidence for helium rain occurring in their interiors (Stevenson 2020). For example, the presence of a heavy-element gradient, already suggested by Stevenson (1985) and supported by recent interior models constrained by Juno and Cassini measurements (Wahl et al. 2017; Debras & Chabrier 2019; Mankovich & Fuller 2021; Miguel et al. 2022; Militzer et al. 2022; Howard et al. 2023) could also significantly affect the thermal evolution of giant planets. In this scenario, heat would initially be trapped in the deep interior with the composition gradient acting as a thermal barrier. It could then be released later, which would explain such an excess in luminosity (Leconte & Chabrier 2013; Vazan et al. 2016).

The strongest indication of helium rain in Jupiter is still the in situ measurements of the Galileo probe (von Zahn et al. 1998), which showed a slight depletion of helium (Y_atm = 0.238 ± 0.005) compared to the protosolar value (Y_proto = 0.27) (Asplund et al. 2021) and also neon depletion, which could be sequestered in the falling helium droplets (Wilson & Militzer 2010). Saturn’s atmospheric helium content remains uncertain and needs to be measured by a dedicated atmospheric entry probe (Mousis et al. 2016; Fortney et al. 2023).

Meanwhile, progress in improving our understanding of the H-He phase diagram has enabled the calculation of inhomogeneous evolution models, thereby accounting for helium rain. Such calculations from Hubbard et al. (1999) and Fortney & Hubbard (2003) based on early H-He phase diagrams (Hubbard & Dewitt 1985; Pfaffenzeller et al. 1995) confirmed that helium rain could explain Saturn’s measured luminosity. Extensive ab initio calculations were then conducted to model the behaviour of H-He mixtures at high pressures. Lorenzen et al. (2009, 2011) performed molecular dynamic simulations based on density functional theory (DFT-MD) on a large range of H-He mixtures, using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation (XC) functional and the ideal entropy of mixing approximation. The resulting phase diagram was employed in evolution calculations of Jupiter (e.g., Nettelmann et al. 2015; Mankovich et al. 2016) and Saturn (e.g., Püstow et al. 2016). The ideal entropy of mixing approximation was avoided by Morales et al. (2009, 2013), leading to lower de-mixing temperatures. Schöttler & Redmer (2018) provided the latest theoretical phase diagram, using the van der Waals density functional (vDW-DF) which led to even lower de-mixing temperatures. Recently, Mankovich & Fortney (2020) searched for consistent solutions for Jupiter’s and Saturn’s evolution by slightly shifting the phase diagram in temperature calculated by Schöttler & Redmer (2018). This was only possible by assuming an enhanced Bond albedo for Saturn. Surprisingly, a laser-shock experiment (Brygoo et al. 2021) yielded much higher de-mixing temperatures than the theoretical predictions (see Fig. 1).

Understanding giant planets requires a good comprehension of the H-He phase diagram, yet it remains very uncertain given the important discrepancy between theoretical calculations and single experiments. In this study, we applied the existing phase diagrams to the evolution of Jupiter and Saturn. We then searched for a coherent solution between the evolution of Jupiter and Saturn since the same fundamental properties of hydrogen and helium are at play. We find that this can be achieved by a temperature shift of the existing phase diagrams and, as we discuss below, this could provide hints about the location and shape of the exact phase diagram. Such an analysis may also be able to constrain the internal structure of these planets.

Fig. 1 Comparison of immiscibility curves from H-He phase diagrams. The thick solid lines show the phase curves for a solar mixture from ab initio calculations of Lorenzen et al. (2011) and Schöttler & Redmer (2018) and from the laser-shock experiment of Brygoo et al. (2021). The thin solid lines show the same phase curves but with temperature offsets of -1250, +350, and -3850 K, respectively. These offsets correspond to the ones used in Sect. 3.1.1, which are required to fit Jupiter’s evolution to Galileo’s measurement of its atmospheric helium content. The dashed lines correspond to the adiabats of Jupiter and Saturn today (see Sect. 3 for details). As the region of de-mixing was considered only at pressures greater than 1 Mbar in our models, the region of the Brygoo2021 phase diagram below 1 Mbar was not used (shown with the dotted red line).

2 Methods

2.1 H-He phase diagrams

In this study, we used the theoretical phase diagrams of Lorenzen et al. (2011) (hereafter, Lorenzen2011) and Schöttler & Redmer (2018) (SR2018), as well as a constructed phase diagram based on experimental data Brygoo et al. (2021) (Brygoo2021). These three phase diagrams are shown (for a H-He mixture in protoso-lar proportions) in Fig. 1. These phase diagrams are then coupled to the evolution calculations of Jupiter and Saturn. This requires us to have access to the phase curves at several helium fractions, x_He, and in various P − T conditions. Therefore, we linearly interpolated the data from Lorenzen et al. (2011) and Schöttler & Redmer (2018) to obtain the de-mixing temperatures on a rectangular grid with 301 values of x_He between 0 and 1 and 301 values of P between 1 and 24 Mbar. The Brygoo et al. (2021) data provides only one curve at x_He = 0.11 (Y = 0.33). We constructed a more complete phase diagram (on the same grid as Loren-zen2011 and SR2018) in the following way. We used the SR2018 phase diagram and estimate the temperature differences at each pressure value between the curve at x_He = 0.11 and the curves at all other x_He values. Next, we applied these temperature differences to the single curve from Brygoo et al. (2021). Figure 1 also shows the phase curves of the Lorenzen2011, SR2018 and Brygoo2021 phase diagrams shifted by temperature offsets of T_offset = −1250, +350 and −3850 K, respectively.

These offsets must be applied to infer a depletion of helium in the atmosphere of Jupiter that is consistent with the Galileo measurement by the evolution model (more details in Sect. 3.1). The value of these offsets depend on the used H-He equation of state. Given the high discrepancy between the three original phase diagrams employed here, it is appropriate to allow for temperature offsets when studying the evolution of Jupiter and Saturn (e.g., Nettelmann et al. 2015; Mankovich et al. 2016; Püstow et al. 2016; Mankovich & Fortney 2020).

2.2 Calculation of the helium profile

Applying H-He phase diagrams to the evolution of giant planets allows for calculating the helium distribution throughout their evolution. We follow the procedure of an earlier work, which is aptly described for instance, Nettelmann et al. (2015) and Mankovich et al. (2016) and also illustrated in Figure 2 of Mankovich & Fortney (2020). The main ideas are summarized below.

Starting from an initial “hot and puffed” state, giant planets cool down and their pressure-temperature profile eventually intersects with a region of the phase diagram, which predicts de-mixing between hydrogen and helium (see Fig. 1). Helium droplets start to form in the intersecting region and then sink to deeper regions. Helium can then mix again in the deeper region if de-mixing is not expected. Such a mechanism is possible due to the hierarchy on different timescales (Stevenson & Salpeter 1977b; Mankovich et al. 2016), namely: the diffusion timescale of the helium droplets (~10⁻¹ s), timescale of convective mixing (~yr), and the timescale corresponding to the long-term thermal planet evolution (~Gyr). The region of the planet where phase separation occurs is depleted in helium until it achieves exact saturation, namely, this means that it has reached a new equilibrium abundance where no further separation occurs. As the planet continues to cool down, there might be a new equilibrium abundance and, thus, helium can be depleted further. The region above phase separation will also be depleted as it supplies helium to the region undergoing phase separation when convective mixing occurs. The inner part of the planet is progressively enriched in helium.

In the following, we describe how we proceeded with our evolutionary models. At every time step and in each layer of our model, the temperature T(P, x_He) is compared to the predicted de-mixing temperature, T_de–mix(P, x_He). As we also allowed for a temperature offset on the phase diagram, if T(P, x_He) < T_de–mix(P, x_He) + T_offset, de-mixing occurs. We calculated the helium concentration for which both terms of the inequality are equal to one another and we find the new helium abundance, x_He, after de-mixing (for exact saturation). The region undergoing phase separation yields one minimum value, x′_He,min, at a pressure denoted as P₀. As our evolution timestep (~ 10-100 Myr) is larger than the convective mixing timescale, we assign x′_He,min to every layer with P ≤ P₀. The shape of the phase diagram then naturally builds a helium gradient in the phase separating region, where P > P₀. The helium lost by de-mixing is redistributed uniformly in the region below the location where phase separation occurs. If de-mixing occurs in the very deep interior near the planetary center or if enough helium sinks, a helium ocean forms (i.e., a region dominated by helium).

2.3 Evolution model setup

We calculated the evolution of Jupiter and Saturn using CEPAM (Guillot & Morel 1995) with our key objective being to find evolution models that reproduce the measured atmospheric helium abundance at present day. We solved the classical structure equations (see, e.g., Guillot 2005; Helled 2018; Miguel & Vazan 2023; Helled & Howard 2024) as follows: $$\frac{\partial P}{\partial r}=-\rho g,$$(1) $$\frac{\partial m}{\partial r}=4\pi r^2\rho ,$$(2) $$\frac{\partial T}{\partial r}=\frac{\partial P}{\partial r}\frac{T}{P}{\nabla }_T,$$(3) $$\frac{\partial L}{\partial r}=-4\pi r^2\rho T\frac{\partial S}{\partial t},$$(4)

where P is the pressure, r is the radius, ρ is the density, g is the gravitational acceleration, m is the mass, T is the temperature, ${\nabla }_T=\frac{\text{d}\ln T}{ \text{d}\ln P}$ is the temperature gradient, L is the intrinsic luminosity, and S is the specific entropy (noting that we neglected the rotation). Equation (4) is particularly relevant in the context of evolution calculations. The change in composition over time due to helium differentiation (Sect. 2.2) is accounted for in the entropy time derivative term.

We describe here the setup of our baseline models (presented in Sect. 3.1). We modeled the planets based on a central core, made of 50% rocks and 50% ices, surrounded by a homogeneous H-He envelope that follows protosolar proportions. All the heavy elements are concentrated in the core and the core mass should be interpreted as a total heavy-element mass (see Mankovich et al. 2016). The assumption of a pure H-He envelope is consistent with the existing phase diagrams calculations that correspond to pure H-He mixtures with no heavy elements (more discussion in Sect. 4).

Our baseline models use a core mass of 30 M_⊕ for Jupiter and 20 M_⊕ for Saturn, which correspond to the total heavy-element mass inferred from recent interior models based on gravity and seismology data (Mankovich & Fuller 2021; Miguel et al. 2022; Howard et al. 2023). Different core masses have been tested and results are presented in Sect. 3.2.3. The assumed structure of our models is simplified; the need to consider more realistic distributions of heavy elements is discussed in Sect. 4. The helium gradient formed by helium de-mixing is kept adiabatic in our models. However, a super-adiabatic temperature gradient may remain stable against convection in the presence of a composition gradient (Ledoux 1947). Such a possibility is discussed in Sect. 3.2.4.

Prior evolution models of Jupiter and Saturn with helium rain used outdated equations of state (EOSs), which do not properly account for H-He interactions (e.g., Nettelmann et al. 2015; Mankovich et al. 2016; Püstow et al. 2016). Mankovich & Fortney (2020) used the MH13 EOS (Militzer & Hubbard 2013), which includes non-ideal mixing effects; however for a constant composition (Y = 0.245). Our baseline models use the CMS19+HG23 EOS, which combines pure hydrogen and helium tables from Chabrier et al. (2019) with non-ideal mixing effects from Howard & Guillot (2023). The specificity of this EOS is that non-ideal mixing effects depend on composition. We recall that the corrections on density or entropy are defined as ΔV(X, Y) = XYV_mix and ΔS (X, Y) = XYS _mix, where X and Y are the mass fractions of hydrogen and helium and V_mix and S_mix are tabulated quantities from Howard & Guillot (2023).

Furthermore, we used the atmospheric models of Fortney et al. (2011). As in Mankovich & Fortney (2020), we recalculated the Jupiter table to make it consistent with the Cassini Bond albedo measurement (Li et al. 2018). We also accounted for the evolution of the Sun’s luminosity by interpolating linearly from 0.7 L_⊙ at 0 Gyr to 1.0 L_⊙ at 4.56 Gyr. The effects of the EOS and of the atmospheric model on evolution calculations are presented in Sects. 3.2.1 and 3.2.2, respectively.

Figure 2 summarizes the procedure we followed. Once the input parameters of the model (planetary mass, core mass, and temperature offset) were set, we calculated the planetary evolution, using a chosen EOS, atmospheric model, and phase diagram. At the end of the evolution, we compared the calculated atmospheric helium mass fraction at t = 4.56 Gyr (the current-age of the planets) to the measured value (Y_atm = 0.238 ± 0.005 from Galileo, von Zahn et al. 1998). If it is not within the uncertainty bounds, we modified the temperature offset and run again the evolution calculation until the measured value is obtained. The radius and the effective temperature as a function of time were also calculated, along with the output parameters. As we present in Sect. 3.2.3, the radius can be adjusted by changing the core mass (or atmospheric model). In Appendix B, we show the results of the evolution models when we used the constraint on the 1 bar temperature instead of the planetary age. The measured values of effective temperature (T_eff), radius (R), age, and atmospheric abundance of helium (Y₁) are similar to those in Mankovich & Fortney (2020) and we reproduce their table here (see Table 1).

Fig. 2 Flowchart of our evolution calculations. The input parameters Mp,Mcore,Toffset, 𝓡p are the planetary mass, the core mass, the temperature offset, and the density ratio, respectively (see Sect. 3.2.4). CEPAM calculates the atmospheric helium mass fraction, the planetary radius, and the effective temperature as a function of time. If the inferred atmospheric helium mass fraction at the present time (4.56 Gyr) does not fall within the uncertainty range of the measured value by Galileo, we updated Toffset and repeated the procedure.

3 Results

3.1 Baseline models

3.1.1 Jupiter

We first ran the evolution models for Jupiter, coupled with the original Lorenzen2011, SR2018, and Brygoo202l phase diagrams. We identified the required temperature offsets and applied them to the original phase diagrams, so that evolution calculations yield an atmospheric helium abundance that is in agreement with the Galileo measurement (Y_atm = 0.238 ± 0.005). The results are shown in Fig. 3. The original Lorenzen2011 phase diagram leads to too much de-mixing at t = 4.56 Gyr, while the SR2018 indicates that de-mixing has not started yet. These results are consistent with previous models from Nettelmann et al. (2015), Mankovich et al. (2016) and Mankovich & Fortney (2020). Interestingly, the Brygoo202l phase diagram yields Y₁ = 0.04 after only 1.5 Gyr. In this case, de-mixing starts already 500 Myr after Jupiter’s formation and we would expect an overly large depletion of helium in Jupiter’s outer envelope in the present era. The Galileo probe measured only a slight depletion in Jupiter’s atmosphere compared to protosolar. We therefore conclude that the original Brygoo2021 phase diagram appears to be very inconsistent with Jupiter’s evolution.

By shifting the Lorenzen2011, SR2018 and Brygoo2021 phase diagrams with adequate offsets (−1250, +350 and −3850 K, respectively), we can fit Jupiter’s evolution to Galileo’s measurement. These offsets are determined with an accuracy of approximately ± 50 K, as variations within this range ensure that the inferred atmospheric helium abundance stays within the uncertainty bounds of Galileo’s measurement. We stress again that the values of these offsets hold only in conjunction with the chosen setup and notably the applied EOSs. In comparison, Nettelmann et al. (2015), using the Lorenzen2011 phase diagram, only required T_offset = −200 K to match Jupiter’s observed atmospheric helium abundance. However, they used the older SCvH EOS (Saumon et al. 1995), which leads to a warmer planetary adiabat by about 1000 K. Similarly, Mankovich et al. (2016) inferred an offset of the same order of magnitude. Using the SR2018 phase diagram and the MH13 EOS, Mankovich & Fortney (2020) required T_offset = 539 ± 23 K. The H-He EOS is crucial as it dictates the interior temperature; hence, it affects the intersection with the phase diagram.

Here, we only intended to constrain the phase diagrams (using T_offset) by fitting the outer envelope helium abundance yielded by Jupiter’s evolution at t = 4.56 Gyr to Galileo’s value. Still, we found an effective temperature for Jupiter that is relatively close to the measured value. We found T_eff lower than the measured value by up to 1.5 K depending on the applied phase diagram. Homogeneous models that do not account for helium rain predict T_eff lower by 2.7 K. Such homogeneous models predict a cooling time of ~4.0 Gyr. For the three shifted phase diagrams, we find a cooling time slightly lower than the age of the planet, by only about ~0.2 Gyr, which is smaller than the uncertainty due to some modeling aspects (Sect. 3.2). The onset of helium differentiation occurs at ages between 3.5 and 3.7 Gyr, in agreement with Nettelmann et al. (2015), at 3.7 Gyr, and Mankovich & Fortney (2020), 3.6–4.0 Gyr. The helium gradient in Jupiter remains small (from Y₁ ~ 0.238 to Y₂ ~ 0.28). It covers a narrow range of pressures and a small portion of the planet (from about 1–3 Mbar), in line with Nettelmann et al. (2015) and Mankovich et al. (2016). It should be noted that not only the location of the phase diagram is important but also its shape. The Brygoo2021 phase diagram exhibits a plateau at P > 2 Mbar and yields a steeper helium gradient compared to the two other phase diagrams.

3.1.2 Saturn

The same fundamental properties of hydrogen and helium are at play in Jupiter and Saturn. Therefore, to model the evolution of Saturn, we then applied the constrained phase diagrams, with the offsets required by Jupiter’s evolution in line with the Galileo measurement. The results are shown in Fig. 4.

As expected, because Saturn’s interior is colder than Jupiter’s, helium de-mixing is more pronounced. We found a depleted atmosphere (and hence outer envelope) with Y1 values of 0.16, 0.13, and 0.14 for the shifted Lorenzen2011, SR2018, and Brygoo2021 phase diagrams, respectively. These values are consistent with the recent estimates from Cassini, namely, with the lower bound of Koskinen & Guerlet (2018) (0.16–0.22) and the upper bound of Achterberg & Flasar (2020) (0.075–0.13). The models from Nettelmann et al. (2015) yielded similar Y₁ values (0.12–0.15) by shifting the Lorenzen2011 phase diagram and using the SCvH EOS while Mankovich & Fortney (2020) found values below 0.10 by shifting the SR2018 phase diagram and using the MH13 EOS (see Sect. 3.2.1).

Right above the heavy-element core, a region with Y = 1 (referred to as the helium ocean) was formed, due to helium sinking and accumulating there. For simplicity, we assume Y = 1 in this helium ocean but He-rich phases of Y = 0.8–0.9 are predicted by the Lorenzen2011 and SR2018 phase diagrams and the presence of heavy elements, not accounted here, would also suggest a region, where Y is not exactly 1. Models from Püstow et al. (2016) and Mankovich & Fortney (2020) also harbour a helium ocean, where the mass fraction of helium is about 0.9 and 0.95, respectively. The extent of this helium ocean depends on the chosen phase diagram, but in all cases, it is confined within the range of 20–30% of Saturn’s total mass. A significant helium gradient, from Y₁ ~ 0.1 to Y₂ ~ 1, covers a wide range of pressures (from about 1–9 Mbar) and a large part of the planet (30–70% by mass). The width of the helium gradient is within the value found by Mankovich & Fortney (2020), namely: 2–14 Mbar. However, it also depends on the presence and the size of a pure heavy-element core (see Sect. 3.2.3). As for Jupiter, the steepness of the helium profile depends on the shape of the phase diagram, with Brygoo2021 yielding a different helium distribution than the Lorenzen2011 and SR2018 ones that are more comparable.

The obtained effective temperatures are in quite good agreement with the measurements. We found T_eff values lower than the measured value by 2.0 and 1.9 K for the shifted Loren-zen2011 and Brygoo2021 phase diagrams, respectively. The shifted SR2018 phase diagram leads to T_eff values that are larger by 0.9 K. Homogeneous models predict T_eff lower by about 10 K and therefore a cooling time that is off by about 2 Gyr (as mentioned in Sect. 1). The shifted Lorenzen2011 and Brygoo2011 phase diagrams lead to a cooling time slightly lower than the age of the planet, by 0.5 Gyr at most. On the contrary, the shifted SR2018 phase diagram yields a cooling time slightly larger, by about ~0.l Gyr. Consequently, the goal of finding a consistent solution for the evolution of both Jupiter and Saturn by applying the same H-He phase diagram is nearly attained. We found the onset of helium differentiation at about 1.3 Gyr quite independently of the applied phase diagram. This is in agreement with the predicted onset at 1-2 Gyr from Püstow et al. (2016) and at 1.5Gyraccording to Mankovich &Fortney (2020).

The potential formation of a helium ocean or helium core was already suggested by Stevenson & Salpeter (1977a) and recent models from Püstow et al. (2016) and Mankovich & Fortney (2020) have supported its presence in Saturn. Figure 5 shows Saturn’s evolution at different times (corresponding to the SR2018 case with T_offset = 350 K). The helium ocean appears −500 Myr after the onset of helium de-mixing. It then keeps growing throughout the evolution. The present-day structure of Saturn is therefore likely to harbour such a helium ocean. We argue that interior models designed to fit Saturn’s gravity field and seismology data should account for the presence of a helium ocean. The thermal and electrical properties of a helium ocean have been investigated by (Preising et al. 2023). We suggest that future studies investigating Saturn’s magnetic field generation and interior dynamics should consider the existence of such a helium ocean.

Fig. 3 Evolutionary calculations of lupiter. Top panels: present-day helium mass fraction as a function of pressure. The pressure range corresponds to the relevant region where de-mixing takes place. Bottom panels: effective temperature as a function of age. Columns (from left to right) correspond to results where we applied the Lorenzen2011, SR2018 and Brygoo2021 phase diagrams, respectively. For comparison, black dashed lines show homogeneous evolution calculations, namely with no de-mixing. Blue solid lines show evolution calculations coupled with the original phase diagrams, namely with no temperature offset. For Brygoo2021, the blue line is dotted and shows the distribution at t = 1.5 Gyr. The orange solid lines show results where a temperature offset, Toffset, was applied to the corresponding phase diagram in order to fit the abundance of helium measured by Galileo (represented by the black dot on top panels). The black dot in the bottom panels corresponds to the measured effective temperature and the errorbar refers to the uncertainty in age.

Fig. 4 Evolutionary calculations of Saturn. Description is same as for Fig. 3. The shaded areas correspond to a “helium ocean", a region above the core where the mass fraction of helium reaches unity. The homogeneous evolution calculations represented by black dashed lines stop because the atmospheric models from Fortney et al. (2011) do not cover lower effective temperatures.

3.2 Model sensitivity to key ingredients

3.2.1 Equation of state

In Sect. 3.1, we showed results using the CMS19+HG23 EOS. The non-ideal mixing effects due to interactions between hydrogen and helium are incorporated and depend on composition (see Sect. 2.3). With helium differentiation, it is crucial to properly assess the thermodynamical quantities as the composition changes with time. The helium mass fraction can span a wide range of values within the interiors of giant planets; in particular, in Saturn (Y ∈ [0.1, 1]). Previous calculations used the MH13 EOS (Mankovich & Fortney 2020). Such an EOS also includes non-ideal mixing effects, but for a constant composition (Y = 0.245). Hence, we ran evolution calculations using the MH13 EOS and compared them to our results with CMS19+HG23. The comparison is shown in Fig. 6. We note that this is a specific case of the SR2018 phase diagram (with T_offset = 350 K). Results using the other phase diagrams are shown in Appendix A.

When we first tried to constrain the H-He phase diagrams by shifting them to fit Jupiter’s evolution to the Galileo measurement, we found that the required temperature offsets are similar (different by less than 150 K). However, the effective temperatures obtained for Jupiter are significantly different from those obtained with CMS19+HG23. Before the onset of de-mixing, CMS19+HG23 yields a colder interior and therefore leads to a shorter cooling time of only about 0.2 Gyr compared to the case with the MH13 EOS. Nonetheless, once de-mixing starts, the discrepancy between the two EOSs increases. Figure 6 shows that while de-mixing starts at about the same age, the slope of the effective temperature with respect to age is different. The break of the curve from the homogeneous evolution is less pronounced with CMS19+HG23. Therefore, it seems that accounting for the composition-dependence of non-ideal mixing effects leads to a lower gain of energy from helium differentiation. The cooling time is thus shortened by about 1 Gyr. For Saturn, CMS19+HG23 yields larger Y₁ values (by 0.03 at most). As a result, models using MH13 are expected to underestimate the current helium content in Saturn’s atmosphere.

Fig. 5 Growth of the helium ocean in Saturn. The calculations come from the results of Saturn’s evolution using the SR2018 phase diagram with Toffset = 350 K (Sect. 3.1.2). Top panel: Helium mass fraction at different ages. The black dashed line corresponds to the protosolar value at Y = 0.27. At t = 2 000 Myr, the helium ocean is already there and starts growing. Bottom panel: temperature-pressure profiles of Saturn’s envelope at these ages. Thick portions of the profiles indicate the helium ocean region. Phase curves from the SR2018 phase diagram with Toffset = 350 K are shown with grey lines for Y = 0.33, 0.27, 0.18, 0.1 and 0.05. Grey-shaded areas correspond to the de-mixing regions.

3.2.2 Atmospheric model

Next, we investigated the sensitivity of the results to the used atmospheric models. We compare the results using the atmospheric model from Fortney et al. (2011) to calculations using the Marley et al. (1999) model. The latter is based on older opacity data and does not account for the evolution of the Sun’s luminosity. The comparison is shown in Fig. 7. We found that the Marley etal. (1999) atmospheric model leads to shorter cooling times, by atleast0.5 Gyr for Jupiter and up to almost1 Gyr for Saturn. The constrained phase diagram applied to Saturn’s evolution leads to a larger mass fraction of helium in the atmosphere, namely: by 0.02. The helium core is thus extended across a smaller range of pressures as less helium sinks to the deep interior. This comparison shows the significant effect of the atmospheric model on the evolution calculations of giant planets.

Fig. 6 Effect of the EOS on the evolution of Jupiter and Saturn. We compare the CMS19+HG23 EOS (Chabrier et al. 2019; Howard & Guillot 2023) and the MH13 EOS (Militzer & Hubbard 2013). The calculations use the SR2018 phase diagram with Toffset = 350 K (Sec 3.1). Top panel: the present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Bottom panel: effective temperature as a function of age. The black dashed lines show homogeneous evolution calculations. The black dots show the measured effective temperature and the errorbar corresponds to the age uncertainty.

3.2.3 Core mass

As mentioned in Sect. 2.3, our baseline models used a core mass M_c of 30 M_⊕ for Jupiter and 20 M_⊕ for Saturn which is consistent with current interior models of the two planets but is still uncertain. Here we investigate the sensitivity of the results to the assumed core masses. We therefore also run calculations with a core mass of 10 M_⊕ for both planets. Using the SR2018 phase diagram and applying T_offset = 350 K, Jupiter’s model with a core mass of 10 M_⊕ yields Y₁ = 0.237 while the model with 30 M_⊕ yields Y₁ = 0.235 (Fig. 8, left panel). Thus, a similar offset value allows both cases to match Jupiter’s evolution to the observed atmospheric helium abundance. This suggests that our results for the temperature offsets are insensitive to the chosen core mass. For Saturn, models with M_c = 10 M_⊕ and 20 M_⊕ also yield similar values of Y₁. However, for both planets, models with different core masses yield different planetary radii (Fig. 8, right panel). This shows that the observed planetary radius can be reconciled by adjusting the core mass, without considerably changing the inferred temperature offset and atmospheric helium abundance. We also show the inferred 1 bar temperature of our models. Our Jupiter models agree with the lower bound of the measured T_1bar, while our Saturn models agree with the upper bound of the measurement. The values of T_1bar obtained when applying the other phase diagrams are consistent with the measurements (see Appendix B). In Appendix B, we show how T_1bar (rather than time) can also be used to constrain the planetary evolution.

For Jupiter, we find similar results for the helium distribution at t = 4.56 Gyr for the two core mass values (Fig. 9). However, for Saturn the inferred upper bound of the helium ocean is different. Decreasing M_c by 10 M_⊕ leads to a present-day helium ocean from P > 6.3 Mbar instead of 7.5 Mbar. This is because the presence of a more massive heavy-element core leads to higher pressure values in the planetary center. As a result, we can conclude that the heavy-element distribution in Saturn affects the location of a helium ocean. In addition, the effective temperatures yielded by calculations with a smaller core are slightly lower (by less than 1 K, in agreement with Mankovich et al. 2016), as these models start colder because they include a smaller mass of heavy elements. The evolution with a core mass of 10 M_⊕ agrees with Jupiter’s mean radius while a core mass of between 10 M_⊕ and 20 M_⊕ would agree with Saturn’s mean radius. Fitting the planets’ radii would hence require a core mass lower than realistic total heavy-element mass inferred from interior models based on the measured gravitational field (and seismology in the case of Saturn) data (Mankovich & Fuller 2021; Miguel et al. 2022; Howard et al. 2023), especially for Jupiter. This test demonstrates that adjusting the core mass (i.e., the heavy-element content) could yield the measured radius without affecting significantly the required T_offset and inferred Y₁ and T_eff values at present-day.

Fig. 7 Effect of the atmospheric model on the evolution of Jupiter and Saturn. We compare the atmospheric models of Fortney et al. (2011) and Marley et al. (1999). The calculations use the SR2018 phase diagram with Toffset = 350 K (Sect. 3.1). Top panel: present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Bottom panel: effective temperature as a function of age. The black dashed lines correspond to homogeneous evolution calculations while the black dots show the measured effective temperature and the errorbar gives the uncertainty in terms of age.

Fig. 8 Comparison of models with different core masses but when applying the same temperature offset. The calculations use the SR2018 phase diagram with Toffset = 350 K. Displayed quantities are at t = 4.56 Gyr. Left panel: atmospheric helium mass fraction as a function of the temperature offset. The horizontal dashed lines show the upper and lower bounds of the observed value of Y1. Right panel: planetary adius as a function of the 1 bar temperature. The black errorbars show the observed radii and 1 bar temperatures of Jupiter and Saturn (see Fig. B.1 for details about these measured values). The dashed lines show the possible range of models with core masses that fall between the two values we considered here.

Fig. 9 Effect of the core mass on the evolution of Jupiter and Saturn. The calculations use the SR2018 phase diagram with Toffset = 350 K. Top panel: present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Middle panel: effective temperature as a function of age. Black dots show the measured effective temperature. Bottom panel: planet’s radius as a function of age. The black dashed lines correspond to homogeneous evolution calculations. Black dots show the measured radii and the errorbar the uncertainty in age.

3.2.4 Super-adiabaticity

As helium rains out, a helium gradient forms in the interiors of Jupiter and Saturn. This helium gradient could inhibit convective mixing. Double-diffusion (Rosenblum et al. 2011; Mirouh et al. 2012; Wood et al. 2013) may occur, resulting in a super-adiabatic temperature gradient. We parameterize the level of super-adiabaticity in the region of the helium gradient using the density ratio, R_ρ. The temperature gradient can then be defined as (see also Mankovich & Fortney 2020; Nettelmann et al. 2015): $${\nabla }_{\text{T}}={\nabla }_{\text{ad}}+R_{\rho }B,$$(5)

with $$B=\frac{{\chi }_{\rho }}{{\chi }_{\text{T}}}{\left(\frac{\text{d\hspace{0.17em}}\ln \rho }{\text{d\hspace{0.17em}}\ln Y}\right)}_{P,T}{\nabla }_Y,$$(6)

where $${\chi }_{\rho }={\left(\frac{\text{d\hspace{0.17em}}\ln P}{ \text{d\hspace{0.17em}}\ln \rho }\right)}_{T,Y};{\chi }_{\text{T}}={\left(\frac{\text{d\hspace{0.17em}}\ln P}{ \text{d\hspace{0.17em}}\ln T}\right)}_{\rho ,Y};{\nabla }_Y=\frac{\text{d\hspace{0.17em}}\ln Y}{ \text{d\hspace{0.17em}}\ln P}.$$(7)

Therefore, the expected change in temperature due to super-adiabaticity depends on the helium gradient and also on the EOS as χ_ρ and χ_T are inferred from the EOS tables. We use R_ρ = 0.05 and R_ρ = 0.1 for Jupiter and R_ρ = 0.004 and R_ρ = 0.008 for Saturn and simulate their evolution with these values. This time, we used the Lorenzen2011 phase diagram (with T_offset = −1250 K), for two reasons. First, it allowed us to find solutions for Saturn. Second, it ensured numerical stability as calculating the helium distribution becomes more challenging when considering super-adiabaticity. Furthermore, super-adiabatic models for Saturn have been smoothed with splines to remove spurious data and focus on the main trend of the evolution.

The simulation results are shown in Fig. 10. As expected, we find that a super-adiabatic region accelerates the cooling of Jupiter and hence yields lower Y₁ values (see, e.g., Mankovich et al. 2016). Our baseline model (adiabatic) already yielded an effective temperature slightly lower than the measured one, as a result, the inclusion of a super-adiabatic temperature gradient due to the helium gradient does not help reconcile the measurement. However, given the sources of uncertainties discussed previously (EOS, atmospheric model), we cannot exclude the existence of such a super-adiabatic region in Jupiter. Previous works suggested so (Fortney & Hubbard 2003; Nettelmann et al. 2015; Mankovich et al. 2016). In particular, Mankovich & Fortney (2020) derived a modest super-adiabaticity (R_ρ ~ 0.05, although the number comparison is not straightforward because R_ρ depends on the EOS and the Y gradient) using MH13 instead of CMS19+HG23, which may have overestimated Jupiter’s cooling time. When using MH13 (see Appendix A), we found T_eff larger than the measured value, allowing us to also infer R_ρ that is on the same order of magnitude as in Mankovich & Fortney (2020). However, the improved CMS19+HG23 EOS prevents us from finding models with a super-adiabatic gradient, as this would lead to lower T_eff. The presence of heavy elements in the envelope in the models of Mankovich & Fortney (2020) also leads to a hotter interior. For Saturn (as discussed above) helium rain begins earlier. We find that super-adiabatic models first accelerate the cooling once de-mixing begins. However, after a certain time (~2 Gyr), the planetary cooling and contraction are delayed. This is because heat is initially stored in the planetary deep interior, which leads in the early stages after the onset of de-mixing to a lower T_eff compared to the adia-batic case (e.g., Leconte & Chabrier 2013). Heat is then released in the outer part of the planet, which leads to a larger T_eff at present. As a result, these Saturn evolution models coupled with Lorenzen2011 and that include super-adiabaticity are more consistent with the measured effective temperature. On the other hand, the SR2018 phase diagram led to slightly higher T_eff with adiabatic models (see Fig. 4). In this case, super-adiabatic models would appear less consistent with the effective temperature measurement. We note that Mankovich & Fortney (2020) found that super-adiabatic (R_ρ < 0.065) and adiabatic models are equally likely for Saturn. We confirm that it is hard to draw clear conclusions regarding the magnitude of super-adiabaticity in Saturn’s interior.

Figure 11 shows the temperature-pressure profiles of the adiabatic and super-adiabatic evolution models of Jupiter and Saturn. One can hence relate the density ratio R_ρ to the increase of temperature in the planetary deep interiors. At t = 4.56 Gyr, the temperature in Jupiter’s outer region is lower when R_ρ = 0.1 compared to the adiabatic case. On the other hand, in Saturn’s outer region, the temperature at the same age is higher when R_ρ = 0.008 compared to the adiabatic case. This confirms the opposite behaviour of Jupiter and Saturn: super-adiabatic models result in a lower T_eff for Jupiter today, while they yield a higher T_eff for Saturn today. Nevertheless, constraining the exact level of super-adiabaticity in Jupiter and Saturn remains challenging with such evolutionary models. The presence of heavy elements in the envelope and the existence of a fuzzy core would also affect the long-term thermal evolution of the planets. Measuring the helium content in Saturn’s atmosphere by a probe could help constrain the evolution of Saturn and put limits on the magnitude of the super-adiabaticity of the helium gradient.

Fig. 10 Effect of super-adiabaticity on the thermal evolution of Jupiter and Saturn. We compare calculations with different values of Rρ (see Eq. (5)). The simulations use the Lorenzen2011 phase diagram with Toffset = −1250 K. Top panel: present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Bottom panel: effective temperature as a function of age. The black dashed lines correspond to homogeneous evolution calculations, while the black dots show the measured effective temperature and the errorbar the uncertainty in age.

Fig. 11 Temperature-pressure profiles of the envelopes of Jupiter and Saturn at different ages, for adiabatic and super-adiabatic interiors. We compare calculations shown on Fig. 10. Dash-dotted lines show T − P profiles of our baseline models, assuming adiabaticity (Rρ = 0). Solid lines show T − P profiles for Rρ = 0.1 for Jupiter (top panel) and Rρ = 0.008 for Saturn (bottom panel).

4 Discussion

4.1 Limitations

While our study presents a step forward towards a better understanding of the evolution of Jupiter and Saturn with helium rain, it also has some limitations. For simplicity, we assumed that all the heavy elements are concentrated in the planets’ cores and the envelopes consist of a pure H-He mixture. Such a heavy-element distribution is oversimplified. First, homogeneous envelopes (even with Z ≠ 0) are at odds with updated interior models of the planets that infer the presence of heavy-element gradients in the deep interiors of both Jupiter and Saturn (e.g., Wahl et al. 2017; Debras & Chabrier 2019; Mankovich & Fuller 2021; Militzer et al. 2022; Howard et al. 2023). Second, the atmospheres of Jupiter and Saturn are expected to be enriched in heavy elements by a factor of ~3 and ~5–10 compared to pro-tosolar for Jupiter and Saturn, respectively (see e.g., Guillot et al. 2023, and references therein). However, the simplified assumed internal structure allows us to consistently couple our evolution models in which the planetary envelopes consist of pure H-He, consistent with available phase diagrams that were constructed for pure H-He mixtures. Since our models include simplifications, the inferred required temperature offsets, T_offset, and the predicted mass fractions of helium in the outer (Y₁) and inner (Y₂) parts of the planetary interior should be taken as indicative of a trend.

We note that constraining the exact level of super-adiabaticity in the formed helium gradient in Jupiter and Saturn remains hard with evolution calculations. Our models suggest that Saturn’s helium gradient could be super-adiabatic. However, the potential presence of a heavy-element gradient could significantly affect the thermal evolution of both planets (Leconte & Chabrier 2013; Vazan et al. 2016). This should be studied further and implemented in future evolution models of Jupiter and Saturn. Theoretical work on convective inhibition (Guillot 1995) could also help to constrain the properties of the super-adiabatic temperature gradient in the phase separating region. Nettelmann et al. (2015) investigated the influence of convection inhibition due to H-He phase separation on Jupiter’s temperature profile and heat transport and constrained the height of semiconvective layers. Future studies could consider magnetic field measurements, which seem to support the presence of a stably stratified layer near the region where phase separation occurs (Connerney et al. 2022; Moore et al. 2022).

Finally, at the moment, the helium atmospheric abundance in Saturn is uncertain. In this work, we found that using the MH13 EOS (Militzer & Hubbard 2013) instead of CMS19+HG23 (Chabrier et al. 2019; Howard & Guillot 2023) leads to Y₁ values reduced to 0.10. A probe to Saturn (Mousis et al. 2016; Fortney et al. 2023) that will measure the helium abundance in the atmosphere could test our model predictions, and constrain Saturn’s evolution and its internal structure.

4.2 Microphysics and atmospheres

The evolution models we presented strongly rely on the EOS which is uncertain. Using a state-of-the-art EOS for H-He yields colder planets. Accounting for non-ideal mixing effects, which depend on the helium fraction (CMS19+HG23 EOS, Howard & Guillot 2023) leads to an energy gain from de-mixing lower than what is obtained with mixing effects for constant composition (MH13 EOS, Militzer & Hubbard 2013). The cooling time can be reduced by several hundred million years in the case of Jupiter when de-mixing has started. The atmospheric model also has an important effect. Using updated atmospheric models, accounting for the evolution of the Sun’s luminosity and Jupiter’s revised Bond albedo, we find longer cooling times than previous atmospheric models, on the order of a few hundred million years. Further investigations of the physical and chemical processes occurring in giant planet atmospheres and the implementation of more realistic atmospheric models (e.g., Chen et al. 2023) will allow for a better understanding of the thermal evolution of giant planets.

An improved understanding of the H-He phase diagram is clearly desirable to further constrain the planetary evolution. At the moment, there is still a large discrepancy between various theoretical calculations and the single experiment result. While our results for the evolution of Jupiter and Saturn provide hints on the location and shape of the exact phase diagram, further explorations of the immiscibility between hydrogen and helium are needed, across a wide range of pressures and temperatures, and for several mixture ratios. In addition, future ab initio simulations and experiments that incorporate heavy elements are required. For example, the presence of helium, even in small quantities, delays the metallization of hydrogen (Vorberger et al. 2007; Mazzola et al. 2018; Helled et al. 2020). The inclusion of heavy elements in theoretical calculations and experiments could therefore influence the phase separation between hydrogen and helium. This should be investigated and implemented in evolution models.

5 Summary and conclusions

We presented thermal evolution models of Jupiter and Saturn accounting for helium rain. We applied theoretical (Lorenzen et al. 2011; Schöttler & Redmer 2018) H-He phase diagrams and our construction of a phase diagram based on the experimental data of Brygoo et al. (2021) and investigated their effect on the evolution of both planets. Our models of Jupiter and Saturn assumed a central dense core surrounded by a H-He envelope consistent with current H-He phase diagrams. Overall, the main conclusions from this study depend on our assumptions (e.g., the applied H-He EOS) and can be summarized as follows:

1. An offset in temperature of −1250, +350, and −3850 K, respectively, is required for the phase diagrams of Lorenzen et al. (2011), Schöttler & Redmer (2018) and Brygoo et al. (2021) to find agreement between Jupiter’s evolution and the Galileo measurement of its atmospheric helium abundance (Y_atm = 0.238 ± 0.005).

2. When applying the same offset in temperature to the evolution of both Jupiter and Saturn, our models of both planets are in rather good agreement with their measured effective temperature (by less than 2 K). This suggests that the exact H-He phase diagram may be obtained by shifting the existing ones with our estimated offsets in temperature (see Fig. 1).

3. The temperature shift implied by the experimental phase diagram (Brygoo et al. 2021) is substantial (−3850 K), suggesting that the original experimental data is very inconsistent with the evolution of both planets. It would lead to a substantial depletion of helium (Y₁ ~ 0.04) in Jupiter’s atmosphere in less than 2 Gyr.

4. Demixing starts at about 3.5 Gyr for Jupiter and at about 1.3 Gyr for Saturn. Helium rain leads to very different interiors between Jupiter and Saturn. Only slight de-mixing (Y₁ = 0.238) in a small part of Jupiter (between 1 and 3 Mbar) occurs, whereas it leads to the presence of a large helium gradient (from Y₁ ~ 0.13 to Y₂ ~ 1) on a large proportion of Saturn (between 2 and 9 Mbar).

5. Saturn is likely to have a helium ocean which forms rapidly in a few hundred million years after the onset of He de-mixing. We find Y ~ 1 in this region, due to simplifying assumptions; however, more realistic values are likely to be around 0.9–0.95 (Püstow et al. 2016; Mankovich & Fortney 2020). The extent of this helium ocean depends on the size of the pure heavy-element core.

6. Saturn’s atmospheric helium mass fraction is between 0.13 and 0.16, consistent with recent Cassini’s estimates, namely: the lower bound of Koskinen & Guerlet (2018) (0.160.22) and the upper bound of Achterberg & Flasar (2020) (0.075-0.13). However, measurements of Saturn’s atmospheric helium content cover a broad range of values (see Fig. B.1) and a more accurate determination is still required.

Our evolution models confirm and are in agreement with previous results. In this work, we find consistent values for the atmospheric helium abundance and similar ages for the onset of helium rain. These findings demonstrate the robustness of our evolution calculations. Overall, evolution models serve as significant complements to static interior models and can unveil important information on the planetary structure, for example, the presence of a helium ocean in Saturn. Future evolution models, including both heavy-element gradients (Vazan et al. 2018; Müller et al. 2020) and H-He phase separation will provide new insights into the evolution and internal structure of giant planets.

Our comprehension of giant planet interiors will also benefit from an improved understanding of the behaviour of hydrogen and helium at high pressures and temperatures, which yields accurate phase diagrams and EOSs. In addition, measuring Saturn’s atmospheric composition by an entry probe is needed to determine the He abundance in Saturn’s atmosphere, constrain the magnitude of He rain, and provide additional clues on the H-He phase diagram and therefore on its formation and evolution history. Furthermore, the detection of warm and cold giant exoplanets allows us to apply our knowledge from the Solar system’s giant planets to characterize giant planets around other stars (Müller & Helled 2023). Finally, further synergies between planetary science, high-pressure physics, and exoplanetary science, as well as advancements from observations, laboratory experiments, and theoretical and numerical calculations will be essential to reveal the nature of gas giant planets in our planetary system and beyond.

Acknowledgements

We thank the referee for valuable comments which helped improve the manuscript. We acknowledge support from SNSF grant 200020_215634 and the National Centre for Competence in Research ‘PlanetS’ supported by SNSF. We thank the Juno team for insightful discussions. We also thank Christopher Mankovich for helpful discussions.

Appendix A Comparison of evolution calculations with different EOSs

We compare in Figs. A.1 and A.2 our results with the MH13 (Militzer & Hubbard 2013) and CMS19+HG23 (Chabrier et al. 2019; Howard & Guillot 2023) EOSs for the evolution of Jupiter and Saturn, respectively.

Fig. A.1 Evolutionary calculations of Jupiter using the CMS19+HG23 and MH13 EOSs. Description is same as for Fig. 3.

Fig. A.2 Evolutionary calculations of Saturn using the CMS19+HG23 and MH13 EOSs. Description is same as for Fig. 4.

Appendix B Constraining evolution with 1 bar temperature

In this work, we used time to constrain our evolution models, seeking solutions that match the measured atmospheric helium abundance at t = 4.56 ± 0.1 Gyr. This time-based approach has been employed in previous studies (Fortney & Hubbard 2003; Mankovich et al. 2016; Püstow et al. 2016; Mankovich & Fortney 2020). However, the 1 bar temperature T_1bar, which serves as an outer boundary condition and defines the entropy in the planetary envelope, can also be used to constrain evolution (Nettelmann et al. 2015; Nettelmann et al. 2024). Indeed, as the planet cools, T_1bar decreases and one can track the planetary evolution by looking at the decrease in T_1bar.

First, we assess whether our baseline models (see Sect. 3.1) agree with the measured 1 bar temperatures. For Jupiter, Galileo measured T_1bar = 166.1 ± 0.8 K (Seiff et al. 1998). However, this was measured in a hot spot; the recent reanalysis of the Voyager data yielded 167.3 ± 3.8 K and 170.3 ± 3.8 K at 0°N and 12°S respectively (Gupta et al. 2022). For Jupiter, our baseline models predict values of the 1 bar temperature of 163.0, 163.0 and 162.6 K for the shifted Lorenzen2011, SR2018 and Brygoo2021 phase diagrams respectively. Our baseline models of Jupiter are therefore in line with the lower bound of the measured T_1bar. For Saturn, Voyager measured T_1bar = 135 ± 5 K (Lindal et al. 1985). We find values of 134.5, 141.0 and 135.2 K for the three shifted phase diagrams respectively, in agreement with the measurement and its upper bound.

We next investigate how the 1 bar temperature constraint affects our conclusions. Figure B.1 shows the inferred atmospheric helium mass fraction as a function of T_1bar. Our baseline model using SR2018 with T_offset = 350 K yielded for Jupiter T_1bar = 163.0 K (at t = 4.5 Gyr), in line with the lower bound of the measured 1 bar temperature. Applying an offset of T_offset = 650 K allows the evolution model to match the upper bound of the measurement (blue dashed line). Using T_1bar rather than the time-based constraint shows that our estimates of the temperature offsets can vary by up to about 300 K. We recall that the uncertainty on the measured atmospheric helium abundance affects those offsets by only ±50 K. Applying the phase diagram with T_offset larger by 300 K to Saturn’s evolution could yield an inferred atmospheric helium abundance different by 0.05.

Similarly, the planetary radius could also be used to constrain evolution. However, the radius does not evolve much during the last billion years (see Fig. 9). Also, since the inferred planetary radius mainly depends on the heavy-element content, matching the radius would only need to adjust the core mass in our models. As shown in Figs. 8 and 9, we find the same values of T_offset to fit the measured atmospheric helium abundance, when using a core mass of 10 or 30 M_⊕ in Jupiter. Therefore, a small radius mismatch in our baseline models does not affect our conclusions.

Fig. B.1 Inferred atmospheric helium mass fraction as a function of the 1 bar temperature. Solid lines correspond to our baseline models for Jupiter and Saturn, using SR2018 with Toffset = 350 K. Dashed lines show similar evolution models but with Toffset = 650 K. The blue error-bar corresponds to the Galileo measurement of Yatm (von Zahn et al. 1998) and the Voyager re-estimation of T1bar (Gupta et al. 2022). Red errorbars (from top to bottom) correspond to estimates from Conrath & Gautier (2000),Koskinen & Guerlet (2018),Achterberg & Flasar (2020), Conrath et al. (1984) and the 1 bar temperature from Lindal et al. (1985). Black dots indicate specific ages in the evolution.

References

All Tables

All Figures

Fig. 1 Comparison of immiscibility curves from H-He phase diagrams. The thick solid lines show the phase curves for a solar mixture from ab initio calculations of Lorenzen et al. (2011) and Schöttler & Redmer (2018) and from the laser-shock experiment of Brygoo et al. (2021). The thin solid lines show the same phase curves but with temperature offsets of -1250, +350, and -3850 K, respectively. These offsets correspond to the ones used in Sect. 3.1.1, which are required to fit Jupiter’s evolution to Galileo’s measurement of its atmospheric helium content. The dashed lines correspond to the adiabats of Jupiter and Saturn today (see Sect. 3 for details). As the region of de-mixing was considered only at pressures greater than 1 Mbar in our models, the region of the Brygoo2021 phase diagram below 1 Mbar was not used (shown with the dotted red line).

In the text

Fig. 2 Flowchart of our evolution calculations. The input parameters Mp,Mcore,Toffset, 𝓡p are the planetary mass, the core mass, the temperature offset, and the density ratio, respectively (see Sect. 3.2.4). CEPAM calculates the atmospheric helium mass fraction, the planetary radius, and the effective temperature as a function of time. If the inferred atmospheric helium mass fraction at the present time (4.56 Gyr) does not fall within the uncertainty range of the measured value by Galileo, we updated Toffset and repeated the procedure.

In the text

Fig. 3 Evolutionary calculations of lupiter. Top panels: present-day helium mass fraction as a function of pressure. The pressure range corresponds to the relevant region where de-mixing takes place. Bottom panels: effective temperature as a function of age. Columns (from left to right) correspond to results where we applied the Lorenzen2011, SR2018 and Brygoo2021 phase diagrams, respectively. For comparison, black dashed lines show homogeneous evolution calculations, namely with no de-mixing. Blue solid lines show evolution calculations coupled with the original phase diagrams, namely with no temperature offset. For Brygoo2021, the blue line is dotted and shows the distribution at t = 1.5 Gyr. The orange solid lines show results where a temperature offset, Toffset, was applied to the corresponding phase diagram in order to fit the abundance of helium measured by Galileo (represented by the black dot on top panels). The black dot in the bottom panels corresponds to the measured effective temperature and the errorbar refers to the uncertainty in age.

In the text

	Fig. 4 Evolutionary calculations of Saturn. Description is same as for Fig. 3. The shaded areas correspond to a “helium ocean", a region above the core where the mass fraction of helium reaches unity. The homogeneous evolution calculations represented by black dashed lines stop because the atmospheric models from Fortney et al. (2011) do not cover lower effective temperatures.
In the text

Fig. 5 Growth of the helium ocean in Saturn. The calculations come from the results of Saturn’s evolution using the SR2018 phase diagram with Toffset = 350 K (Sect. 3.1.2). Top panel: Helium mass fraction at different ages. The black dashed line corresponds to the protosolar value at Y = 0.27. At t = 2 000 Myr, the helium ocean is already there and starts growing. Bottom panel: temperature-pressure profiles of Saturn’s envelope at these ages. Thick portions of the profiles indicate the helium ocean region. Phase curves from the SR2018 phase diagram with Toffset = 350 K are shown with grey lines for Y = 0.33, 0.27, 0.18, 0.1 and 0.05. Grey-shaded areas correspond to the de-mixing regions.

In the text

Fig. 6 Effect of the EOS on the evolution of Jupiter and Saturn. We compare the CMS19+HG23 EOS (Chabrier et al. 2019; Howard & Guillot 2023) and the MH13 EOS (Militzer & Hubbard 2013). The calculations use the SR2018 phase diagram with Toffset = 350 K (Sec 3.1). Top panel: the present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Bottom panel: effective temperature as a function of age. The black dashed lines show homogeneous evolution calculations. The black dots show the measured effective temperature and the errorbar corresponds to the age uncertainty.

In the text

Fig. 7 Effect of the atmospheric model on the evolution of Jupiter and Saturn. We compare the atmospheric models of Fortney et al. (2011) and Marley et al. (1999). The calculations use the SR2018 phase diagram with Toffset = 350 K (Sect. 3.1). Top panel: present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Bottom panel: effective temperature as a function of age. The black dashed lines correspond to homogeneous evolution calculations while the black dots show the measured effective temperature and the errorbar gives the uncertainty in terms of age.

In the text

Fig. 8 Comparison of models with different core masses but when applying the same temperature offset. The calculations use the SR2018 phase diagram with Toffset = 350 K. Displayed quantities are at t = 4.56 Gyr. Left panel: atmospheric helium mass fraction as a function of the temperature offset. The horizontal dashed lines show the upper and lower bounds of the observed value of Y1. Right panel: planetary adius as a function of the 1 bar temperature. The black errorbars show the observed radii and 1 bar temperatures of Jupiter and Saturn (see Fig. B.1 for details about these measured values). The dashed lines show the possible range of models with core masses that fall between the two values we considered here.

In the text

Fig. 9 Effect of the core mass on the evolution of Jupiter and Saturn. The calculations use the SR2018 phase diagram with Toffset = 350 K. Top panel: present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Middle panel: effective temperature as a function of age. Black dots show the measured effective temperature. Bottom panel: planet’s radius as a function of age. The black dashed lines correspond to homogeneous evolution calculations. Black dots show the measured radii and the errorbar the uncertainty in age.

In the text

Fig. 10 Effect of super-adiabaticity on the thermal evolution of Jupiter and Saturn. We compare calculations with different values of Rρ (see Eq. (5)). The simulations use the Lorenzen2011 phase diagram with Toffset = −1250 K. Top panel: present-day helium mass fraction as a function of pressure. The black dot shows the Galileo measured value. Bottom panel: effective temperature as a function of age. The black dashed lines correspond to homogeneous evolution calculations, while the black dots show the measured effective temperature and the errorbar the uncertainty in age.

In the text

	Fig. 11 Temperature-pressure profiles of the envelopes of Jupiter and Saturn at different ages, for adiabatic and super-adiabatic interiors. We compare calculations shown on Fig. 10. Dash-dotted lines show T − P profiles of our baseline models, assuming adiabaticity (Rρ = 0). Solid lines show T − P profiles for Rρ = 0.1 for Jupiter (top panel) and Rρ = 0.008 for Saturn (bottom panel).
In the text

	Fig. A.1 Evolutionary calculations of Jupiter using the CMS19+HG23 and MH13 EOSs. Description is same as for Fig. 3.
In the text

	Fig. A.2 Evolutionary calculations of Saturn using the CMS19+HG23 and MH13 EOSs. Description is same as for Fig. 4.
In the text

Fig. B.1 Inferred atmospheric helium mass fraction as a function of the 1 bar temperature. Solid lines correspond to our baseline models for Jupiter and Saturn, using SR2018 with Toffset = 350 K. Dashed lines show similar evolution models but with Toffset = 650 K. The blue error-bar corresponds to the Galileo measurement of Yatm (von Zahn et al. 1998) and the Voyager re-estimation of T1bar (Gupta et al. 2022). Red errorbars (from top to bottom) correspond to estimates from Conrath & Gautier (2000),Koskinen & Guerlet (2018),Achterberg & Flasar (2020), Conrath et al. (1984) and the 1 bar temperature from Lindal et al. (1985). Black dots indicate specific ages in the evolution.

In the text