Jérémy Leconte
email: jeremy.leconte@lmd.jussieu.fr

Gilles Chabrier
email: chabrier@ens-lyon.fr

Layered convection as the origin of Saturn's luminosity anomaly

 . Existing evolution calculations, however, suggest that the energy released by this sedimentation process may not be sufficient to resolve the puzzle 9 . Here, we demonstrate using planetary evolution models that the presence of layered convection in Saturn's interior, generated, like in some parts of Earth oceans, by the presence of a compositional gradient, significantly reduces its cooling. It can explain the planet's present luminosity for a wide range of configurations without invoking any additional source of energy. This suggests a revision of the conventional homogeneous adiabatic interior paradigm for giant planets, and questions our ability to assess their heavy element content. This reinforces the possibility for layered convection to help explaining the anomalously large observed radii of extrasolar giant planets.

 . However, the effect of this gradient on the thermal evolution of the planet is usually neglected for sake of simplicity. When a vertical gradient of heavy elements is present in a fluid, the resulting mean molecular weight gradient (decreasing upward) can prevent large scale convection to develop by counteracting the destabilizing effect of the temperature gradient. The complex interaction between advection and diffusion of 1

heat and atomic concentration can trigger a hydrodynamical instability, called the double diffusive instability [START_REF] Stern | The salt-fountain and thermohaline convection[END_REF] , leading to a regime of double diffusive convection (also called semi convection). This process significantly affects heat and element transport, as observed in Earth oceans. The instability leads to either a state of homogeneous double diffusive convection where diffusive transport is only modestly enhanced by small scale turbulence, or to a state of "layered convection", with numerous, small convective layers separated by thin, diffusive interfaces corresponding to discontinuities in the composition of the fluid, in which transport is enhanced more significantly [START_REF] Radko | What determines the thickness of layers in a thermohaline staircase?[END_REF][START_REF] Rosenblum | Turbulent Mixing and Layer Formation in Double-diffusive Convection: Three-dimensional Numerical Simulations and Theory[END_REF] .

Numerical calculations show that, for relevant thermal and atomic diffusivities, both scenarios are possible, depending on the ratio of the compositional to superadiabatic temperature gradients [START_REF] Mirouh | A New Model for Mixing by Double-diffusive Convection (Semi-convection). I. The Conditions for Layer Formation[END_REF] . In planets, however, this latter quantity is not a free parameter but is imposed by the (measured) energy flux to be transported in the planet. It can be shown that this flux is too high to be transported by the diffusive regime of the double diffusive instability (said differently, the thermal gradient needed to transport the internal flux over the whole planet would be too high to be stabilized by the heavy element gradient [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] ). Under planetary conditions, the system will thus most likely settle in the regime of layered convection.

To study the impact of layered convection on giant planet thermal structure, we have developed an analytical model of the heat transport in such a medium [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] . Extending the usual mixing length formalism [START_REF] Hansen | Stellar Interiors. Physical Principles, Structure, and Evolution[END_REF] to layered convection, we can calculate the internal temperature gradient as a function of the local properties of the fluid and of a unique free parameter, namely the characteristic size of the convective layers, l, or equivalently the corresponding dimensionless mixing length parameter, α = l/H P , where H P is the pressure scale height in the fluid (see Methods). A first confirmation of the viability of this scenario for giant planet interiors has been given by the fact that structure models with layered convection matching all the observational mechanical constraints of Jupiter and Saturn (radius, gravitational moments), as well as the atmospheric helium and mean heavy element abundances, have been obtained for a rather wide range of mixing length parameters, namely α ∈ [10 -2 -4 × 10 -6 ] for Saturn and α ∈ [10 -2 -3 × 10 -5 ] for Jupiter [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] . These constraints suggest that layered convection is favored over homogeneous double diffusive convection (Supplementary Information). However, the question remains whether layered convection can also explain the present luminosity of the giant planets.

In order to answer this question, we have computed the thermal evolution of planets with layered convection. Because departure from adiabaticity can be significant in these interiors, usual isentropic evolution calculations 2 cannot be used. Instead, the evolution is computed by integrating the intrinsic luminosity, L int = -d E tot /dt, where the total energy E tot = E g + E th + E rot includes the gravitational, thermal and rotational energies (see Methods). The size of the convective/diffusive layers is assumed to have reached an equilibrium value [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF][START_REF] Radko | What determines the thickness of layers in a thermohaline staircase?[END_REF] and is thus kept constant. For all the models (i.e. the reference, homogeneous/adiabatic case and the semi convective ones, for any given α), the amount of heavy elements is kept constant and equal to the one that fulfills the gravitational moment constraints [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] (Supplementary Information). ) and intrinsic (T int ; dashed curves) temperature (see Methods) evolution in time for the adiabatic reference model (black) and three models with layered convection (α = 10 -2.5 , 10 -3 and 10 -3.5 from dark to light red). b) Zoom on present era. Dots show the observed effective temperature. At early ages the effective temperature of models with layered convection is lower due to inefficient convection. After a few hundred million years, these models become brighter due to the release of the excess of energy stored from the initial state.

As illustrated in Fig. 1, our results show that starting from the same total internal energy, the time required to release a given amount of this energy is significantly longer in the layered case than in the homogeneous adiabatic one. This does not imply, however, that the intrinsic luminosity of the object will necessarily be larger at any given time. On the contrary, at early ages, objects with a layered convection zone are far less luminous because of the reduced heat flux release. However, after some time, the decrease in luminosity imposed by layered convection is overpowered by the increase of internal energy to be released and planets with layered convection eventually become more luminous than the ones with adiabatic interiors (Fig. 2). In Saturn's case, this "crossing time" occurs before Saturn's present age so that layered convection yields a larger luminosity at 4.6 Gyr. As shown in Fig. 2 and Fig. 3, if the size of the convective layers is about 2-3×10 -3 pressure scale heights, this process properly leads to the currently observed effective temperature, radius, and gravitational moments of the planet without any additional energy source such as helium rains. The radius evolution of these models is shown in Supplementary Fig. 5. For this range of layer sizes, the present interior temperature of Saturn is only modestly affected compared to the adiabatic case (see Supplementary Fig. 6). The temperature during the early evolution, however, is much higher. ), ii) a scenario where only the inner 50% of the envelope exhibits layered convection (M LC ≈ 0.5 M; α = 2 × 10 -4 ), iii) a case with a layered convection shell between 0.5 and 0.7 M (M LC ≈ 0.2 M; α = 10 -4 ).

Fig. 2 suggests that Saturn models with convective/diffusive layers smaller than 2-3×10 -3 pressure scale heights will be too bright at the age of the Solar System. This is not necessarily true. If layered convection develops only within a restricted fraction of the planet, models with a lower convective efficiency (lower α) in the layered zone and with an efficient convection everywhere else can also reproduce Saturn's proper cooling timescale. Our results are valid for very different sizes of the layered convection zone. Without covering the whole parameter space, we illustrate this by showing an evolution track for Saturn where layered convection is restricted to the inner 50% in mass of the planet above the core, with α = 2 × 10 -4 , and an other track where layered convection occurs in a shell between 50% and 70% of the planet's mass (Fig. 3). This latter case shows that the layered convection region does not need to extend down to the core and could be present around the molecular/metallic transition region or near an immiscibility region in the gaseous envelope.

We have also tested the sensitivity of the model to the initial conditions. Indeed, scenarios of giant planet formation by core accretion do not provide so far a clear prescription for these initial conditions. Large uncertainties in various parameters (gas-to-dust ratio, opacities,...) and the lack of a proper treatment of the radiative shock during the gas runaway accretion phase, for instance, hamper a proper determination of the planet's initial energy content [START_REF] Pollack | Formation of the Giant Planets by Concurrent Accretion of Solids and Gas[END_REF][START_REF] Marley | On the Luminosity of Young Jupiters[END_REF] . Whereas these (unknown) initial conditions become inconsequential after about a few tens of million years for an adiabatic planet of Saturn's mass [START_REF] Marley | On the Luminosity of Young Jupiters[END_REF] , it takes significantly longer (up to a few billion years) for a planet where convection is inefficient, because of the lower luminosity and thus the longer thermal, so-called Kelvin-Helmholtz timescale (τ KH ∝ 1/L int ). Initial conditions thus have a stronger impact on the cooling timescale of planets with layered convection.

As shown in Fig. 4, however, these different initial conditions do not qualitatively affect our results, as there is always a possible trade-off between the initial energy content of the nascent planet (at the end of the runaway accretion shock) and the thickness or number of convective-diffusive layers. As seen in the figure, evolution of models with layered convection can adequately reproduce Saturn's presently observed luminosity with very low initial energy content provided sufficiently small convective/diffusive layers. Given the current uncertainties of formation models, thermal evolution of Saturn with layered, inefficient convection can thus meet the observational constraints for a wide variety of initial conditions. Interestingly enough, the same layered convection mechanism, which inhibits the heat flow, also provides a plausible explanation to the observed low luminosity of Uranus [START_REF] Podolak | Model of Uranus' interior and magnetic field[END_REF] . In this case, however, the "crossing time" discussed above should be larger than the age of the Solar System.

Because the present luminosity of Jupiter is smaller than the one that is predicted by homogeneous, adiabatic evolution models 2 (see Supplementary Fig. 7), observations do not seem to support the presence of layered convection over the whole planet at the present epoch. This does not preclude, however, the possibility for this process either to be present in part of the planet or to have occurred in the past, but the planet's gaseous envelope to have been homogenized early enough, due to the more vigorous convective flux. The flux emitted by Jupiter during its cooling history is large enough to redistribute a significant fraction of a massive initial core (of up to 15 -20 M ⊕ ) within less than 4.6 Gyr 21 . Layered convection in this planet could then have stopped after the erosion of the whole core, or at least of its soluble components [START_REF] Wilson | Solubility of Water Ice in Metallic Hydrogen: Consequences for Core Erosion in Gas Giant Planets[END_REF][START_REF] Wilson | Rocky Core Solubility in Jupiter and Giant Exoplanets[END_REF] . This would be consistent with both i) the large initial core mass needed to form the giant planet on the right timescale [START_REF] Stevenson | Formation of the giant planets[END_REF] and ii) Jupiter's present smaller core [START_REF] Guillot | The Interior of Jupiter[END_REF][START_REF] Saumon | Shock Compression of Deuterium and the Interiors of Jupiter and Saturn[END_REF][START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] . Jupiter and Saturn may thus have experienced quite different cooling histories and end up with dissimilar interiors.

Layered convection does not preclude either the possibility of a phase separation between hydrogen and either volatiles or helium in Saturn's interior. In fact, a phase separation favors the occurrence of layered convection [START_REF] Stevenson | The phase diagram and transport properties for hydrogenhelium fluid planets[END_REF] . Indeed, the strong compositional gradients (discontinuities) due to the pres- ence of an immiscibility region provide very favorable conditions for the development of the double diffusive instability and thus of layered convection. Our present treatment of layered convection does not take into account the energetic contribution of a possible phase separation. In that case, however, the energy release due to the phase separation does not need to be important, as the main role of immiscibility would be to trigger the double diffusive instability. Indeed, as both layered convection and phase separation increase the planet's cooling time, there is always a possible trade off between the respective contributions of these two processes to the planet's cooling rate. Note, however, that layered convection yields a hotter interior than the adiabatic evolution, favoring miscibility, especially in the past (see Supplementary Fig. 6). As a smaller immiscibility region, if any, is predicted for helium in Jupiter [START_REF] Lorenzen | Demixing of Hydrogen and Helium at Megabar Pressures[END_REF][START_REF] Morales | Phase separation in hydrogen-helium mixtures at Mbar pressures[END_REF] , this scenario would also explain the dichotomy revealed here between Jupiter and Saturn.

The cooling history of our Solar System giants, and of extrasolar giant planets in general, might thus be more complex than assumed by the conventional, simplistic homogeneous, adiabatic interior structure and evolution paradigm. The amazing diversity of the observed properties of extrasolar planets -in particular the anomalously large observed radii of hot jupiters [START_REF] Miller | The Heavy-element Masses of Extrasolar Giant Planets, Revealed[END_REF] that cannot be fully explained by the various proposed heating mechanisms [START_REF] Laughlin | On the Anomalous Radii of the Transiting Extrasolar Planets[END_REF][START_REF] Leconte | Is tidal heating sufficient to explain bloated exoplanets? Consistent calculations accounting for finite initial eccentricity[END_REF] -combined with the luminosity anomalies of our own Solar System giants, clearly suggest a significant revision of this paradigm and point to a broader, more complex picture of solar and extrasolar giant planet structure, composition and thermal evolution, with a direct impact on giant planet formation conditions.

Methods

Layered convection and thermal gradient

To study the impact of layered double diffusive convection on the interior thermal structure, we have developed an analytical model of the heat transport in such a medium [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] , which is briefly outlined below.

The fluid is assumed to be composed of a large number of small convective layers of height l, separated by thinner diffusive interfaces. The transport in the convective layers is described by a parametrization similar to the standard mixing length theory of convection, for which the mixing length is chosen to be equal to the height of the layers (l). The temperature gradient (∇ T ≡ ∂ ln T ∂ ln P ), hence the superadiabaticity in the gas, ∇ T -∇ ad , is given as a function of the flux to be transported, of the thermodynamic properties of the medium and of l.

In the diffusive interfaces, the thermal gradient is equal to the gradient needed to transport the whole intrinsic flux by diffusive processes alone (∇ d ). The thermal size of these diffusive interfaces, δ T , can be estimated by the fact that the convective timescale, given by the inverse of the Brunt-Väisälä frequency (τ c = N -1 T ) must be equal to the diffusive timescale in the interface (τ d = δ 2 T /κ T , where κ T is the thermal diffusivity of the medium). Thus, δ T = κ T /N T , and the only remaining free parameter is the height of the convective layers, or equivalently the mixing length parameter, α = l/H P , where H P is the pressure scale height in the fluid.

Finally, the mean temperature gradient ∇ T is a linear combination of the two aforementioned temperature gradients weighted by the relative size of the convective and diffusive regions, and depends on the size of the convective/diffusive cells [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] (i.e. on α).

Planetary model sequences and time integration

Thermal evolution tracks of a planet of a given mass M, composition (symbolically denoted by an array, X, encompassing all the information needed to know the abundances of all the materials at all depths, including the core), angular momentum (J), and thermal structure (parametrized here by the mixing length parameter α) yield relations of the type

L int ≡ 4πR 2 σ SB T 4 int = 4πR 2 σ SB T 4 eff -T 4 eq = L int (M, X, α, J,t), E tot ≡ E th + E g + E rot = E tot (M, X, α, J,t), R = R(M, X, α, J,t), (1)
where T int and T eff are the intrinsic and total effective temperatures, σ SB the Stefan-Boltzmann constant and R the radius of the planet. The mixing length parameter is assumed constant as convective/diffusive layers are assumed to have reached a state of equilibrium [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF][START_REF] Radko | What determines the thickness of layers in a thermohaline staircase?[END_REF] . Without any further guidance from tridimensional hydrodynamical simulations covering a large enough spatial and temporal domain, it is sensible to start with the most simple assumption. The equilibrium temperature (T eq ) is the temperature that the planet would have if its intrinsic luminosity were zero and represents the contribution to the solar flux through

σ SB T 4 eq ≡ (1 -A)L /(16πa 2 ), (2)
where L is the solar luminosity, a the orbital distance of the planet, and A its Bond albedo. Throughout this study, the values used for the equilibrium temperature are T eq = 109 K and 91 K for Jupiter and Saturn respectively. This corresponds to an absorbed luminosity of 5.14×10 16 W (1.11×10 16 W) and a Bond albedo of 0.37 (0.32) for Jupiter (Saturn). Within the observational uncertainties, those values are consistent with observations. The relations described by ( 1) can be derived with a grid of atmosphere models. These grids of atmospheric boundary conditions provide us with the temperature at the 10 bar level (or any reference level deeper than the photosphere), T 10 (T int , g), as a function of the intrinsic effective luminosity (T int ), and of the gravity, g ≡ GM/R 2 . The details of the functional form used to represent the atmospheric boundary conditions are detailed in the Supplementary Information.

Then, for any arbitrary value of the luminosity, T 10 is known, and the T -P profile is integrated inward by integrating the standard structure equations for a rotating body (Supplementary Information). Models are computed on a grid of luminosities. For each luminosity, the gravitational, thermal ( ũ being the specific internal energy) and rotational energy,

E g ≡ - Gm r dm, E th ≡ ũ dm, E rot ≡ 1 2 J 2 I , (3)
hence the total energy, are computed (see Supplementary Fig. 8). Rotational energy is computed assuming negligible angular momentum loss, and is found to represents only a small fraction of the total energy (≈ 2%). As the ANEOS equation of state provides us with ũ (the specific internal energy) for the volatiles, the thermal contribution of the core is included in our calculations. Finally, we are left with a tabulated form of functions of the type E tot = E tot (M, X, α, J, L int ), and time dependence can be retrieved by choosing an initial value for the luminosity or initial energy content and by integrating

-L int = ∂ E tot ∂ t M, X,α,J . (4)
∂ P ∂ m = - 1 4π Gm r4 + ω 2 6π r + ϕ ω (r), (6)
∂ r

∂ m = 1 4π r2 ρ , (7)
∂ T ∂ m = T P ∂ P ∂ m ∇ T , (8)
where m is the mass enclosed in the equipotential of mean radius r, ω the rotation rate of the planet, ϕ ω (r) a second-order correction due to the centrifugal potential. As discussed above, ∇ T is either equal to ∇ ad in the adiabatic case, or to ∇ T which depends on α in a layered convection zone. The angular velocity is given by the fact that J ≡ Iω, I being the angular moment of inertia computed for each model, is kept constant and equal to the present value. A closure equation is provided by the equation of state (EOS) of the mixture along the planet's interior profile. Such an EOS is generally given by the so-called ideal volume law for the mixture:

1 ρ = X ρ X + Y ρ Y + Z ρ Z , (9)
where X, Y and Z denote the mass fractions of H, He and heavy elements, respectively. For each value of α, the core mass and the mass fraction of heavy elements as a function of depth are considered constant in time and equal to the distribution that best matches the gravitational moments and observational constraints [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] . For adiabatic, homogeneous models, the core masses and uniform metal enrichment in the envelope are equal to ∼ 4 M ⊕ and Z = 0.11 in Jupiter and 25.5 M ⊕ and Z = 0.05 in Saturn. For models of Saturn with layered convection, the metal enrichments used for a given mixing length parameter α are shown in [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] and the core mass used is ∼ 20.5 M ⊕ . This difference in the core mass between the homogeneous model and the layered ones is responsible for the small energy difference remaining at low luminosity in Supplementary Fig. 8. Details about the procedure and the equation of states used can be found in ref. 13 and reference therein. Supplementary Fig. 7 shows adiabatic, homogeneous evolution tracks for Jupiter and Saturn computed with this method. As expected, in the homogeneous case, our calculations are in good agreement with previous results 2 : models cool to the observed effective temperature in ∼ 2.9 Gyr for Saturn, and ∼ 5.2 Gyr for Jupiter. As we use the same atmospheric model grids, the slight differences arising in the cooling times are probably due to the fact that our adiabatic models do not have exactly the same core mass as the most recent calculations 2 , and that we take into account the thermal contribution of this core.

Age Gyr

T int K

Figure 7: Effective (T eff ; solid curves) and intrinsic (T int ; dashed curves) temperature as a function of the planetary age for our adiabatic reference models of Jupiter (black) and Saturn (gray). Observed effective temperatures are also shown as dots. As seen, our calculations in the conventional homogeneous adiabatic case are in good agreement with previous results 2 : models cool to the observed effective temperature in ∼ 2.9 Gyr for Saturn, and ∼ 5.2 Gyr for Jupiter. , 10 -3 and 10 -3.5 sequences, respectively. As expected, for a given internal energy content, the intrinsic luminosity of models assuming layered convection decreases when α decreases and is significantly reduced compared with the adiabatic case. This increases the Kelvin Helmholtz timescale, τ KH ≡ GM 2 /(RL int ), i.e. the time needed for the object to loose the memory of its initial energy state.

Figure 1 :

 1 Figure1: Evolution of the total internal energy (E tot ≡ E g + E th + E rot ) with time for various Saturn models. The solid curve corresponds to the reference adiabatic model. The long dashed, dashed and dotted curves are models with layered convection with α = 10 -2.5 , 10 -3 and 10 -3.5 , respectively. Because of the reduced intrinsic luminosity caused by layered convection, these models cool more slowly and keep the memory of their initial energetic state much longer.

Figure 2 :

 2 Figure2: Cooling sequences of Saturn models with layered convection. a) Effective (T eff ; solid curves) and intrinsic (T int ; dashed curves) temperature (see Methods) evolution in time for the adiabatic reference model (black) and three models with layered convection (α = 10 -2.5 , 10 -3 and 10 -3.5 from dark to light red). b) Zoom on present era. Dots show the observed effective temperature. At early ages the effective temperature of models with layered convection is lower due to inefficient convection. After a few hundred million years, these models become brighter due to the release of the excess of energy stored from the initial state.

Figure 3 :

 3 Figure3: Impact of the size of the layered convection zone on Saturn's cooling sequence. a) Effective (T eff ; solid curves) and intrinsic (T int ; dashed curves) temperature evolution for the adiabatic model (black curve) and three scenarios with layered convection. b) Zoom on the present era. From darker to lighter red curves we have, i) our baseline scenario with layered convection present throughout the gaseous envelope of the planet (M LC = 0.8M; α = 2 × 10 -3 ), ii) a scenario where only the inner 50% of the envelope exhibits layered convection (M LC ≈ 0.5 M; α = 2 × 10 -4 ), iii) a case with a layered convection shell between 0.5 and 0.7 M (M LC ≈ 0.2 M; α = 10 -4 ).

Figure 4 :

 4 Figure4: Past luminosity of Saturn's model with layered convection. Evolutionary tracks of the intrinsic luminosity (in units of the presently measured luminosity) of Saturn models integrated backward in time from presently observed conditions with different mixing length parameters (α). All the models with α > 2 × 10 -3 have too short cooling timescales and cannot be integrated backward for 4.6 Gyr. Even for quite low initial energy contents, i.e. a low initial intrinsic luminosity, models with a sufficiently inefficient convection can reproduce the currently observed luminosity of Saturn.

Figure 5 :

 5 Figure5: Temporal evolution of the equatorial radius (in units of the measured present equatorial radius of Saturn, 60 268 km) of an adiabatic Saturn model (solid black curve), and of three models with layered convection (α = 3 × 10 -3 , 10 -3 and 3 × 10 -4 from dark to light red). The effect of the deformation due to the fast rotation is taken into account. The time needed for the adiabatic model fulfilling the constraints on the measured gravitational moments to cool down and contract to the present radius is too short.

Figure 6 :

 6 Figure 6: Internal temperature-pressure profile of our baseline model with layered convection (red curves; M LC ≈ 0.8 M; α = 2 × 10 -3 ) at three different ages (dotted: 0.1 Gyr; dashed: 1 Gyr; long dashed: 4.6 Gyr). The solid black curve represents the temperature profile inferred for the standard adiabatic model. For comparison, the 1 Gyr profile of an adiabatic evolution is shown by the gray dashed curve.

Figure 8 :

 8 Figure8: Intrinsic luminosity (L int ≡ 4πR 2 F int ) as a function of the total internal energy (E tot ≡ E g + E th + E rot ) for various sequences of Saturn models. The solid curve is the reference adiabatic model. The long dashed, dashed and dotted curves are the α = 10 -2.5 , 10 -3 and 10 -3.5 sequences, respectively. As expected, for a given internal energy content, the intrinsic luminosity of models assuming layered convection decreases when α decreases and is significantly reduced compared with the adiabatic case. This increases the Kelvin Helmholtz timescale, τ KH ≡ GM 2 /(RL int ), i.e. the time needed for the object to loose the memory of its initial energy state.

Table 1 :

 1 Parameters used in Eq. (5) to fit the atmospheric grids of ref. 2.When T 10 is known, the T -P profile is integrated inward by integrating the standard structure equations for a rotating body

		Jupiter	Saturn
	T mid (K)	224		225
	∆T (K)	25		25
		T < T mid T > T mid T < T mid T > T mid
	C	78.4	-283.	62.7	-1680.
	K	0.00550	198.	0.00202	831.
	T 0	-122.	143.	-97.1	141.
	β	-0.182	-0.114	-0.180	-0.0721
	γ	2.09	0.454	2.31	0.293
	3 Planetary structure integration

Acknowledgments J.L. thanks J. Fortney for making his atmospheric grids available to us in electronic format. The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 247060)

Author contributions

J.L. carried out analytical calculations, developed the model and performed the numerical simulations. G.C. suggested the idea and carried out analytical calculations. J.L and G.C. wrote the manuscript.

Supplementary Information

1 Layered convection or homogeneous double diffusive convection?

The transport efficiency derived from numerical simulations [START_REF] Rosenblum | Turbulent Mixing and Layer Formation in Double-diffusive Convection: Three-dimensional Numerical Simulations and Theory[END_REF][START_REF] Mirouh | A New Model for Mixing by Double-diffusive Convection (Semi-convection). I. The Conditions for Layer Formation[END_REF] for a homogeneous double diffusive medium (where transport is ensured by small scale turbulence) is equivalent, in our analytical model, to a medium in a state of layered convection with α ≈ 10 -8 -10 -9 (see Figure 3 of Leconte and Chabrier [START_REF] Leconte | A new vision of giant planet interiors: Impact of double diffusive convection[END_REF] ).

The fact that measured gravitational moments constrain α to be larger than 10 -6 ( 13 ) thus seems to favor the occurrence of layered convection over homogeneous double diffusive convection in Solar System giant planets. In addition, the evolutionary constraint on the luminosity discussed in the present study, which constrain α to be larger larger than 10 -4 (see Fig. 3), also seems to justify this hypothesis.

Atmospheric model grids

The relations described by (1) can be derived with a grid of atmosphere models. These grids of atmospheric boundary conditions provide us with the temperature at the 10 bar level (or any reference level deeper than the photosphere), T 10 (T int , g), as a function of the intrinsic flux exiting the planet (parametrized by the intrinsic effective luminosity, T int ), and of the gravity, g ≡ GM/R 2 .

We use the recent atmospheric grids derived specifically for Jupiter and Saturn by Fortney et al. [START_REF] Fortney | Self-consistent Model Atmospheres and the Cooling of the Solar System's Giant Planets[END_REF] . However, as a planet with layered convection can be quite hot (high internal temperature) but with a low luminosity (low effective temperature), our evolutionary tracks can start in the low g, low T int part of the parameter space that is usually not probed by adiabatic evolution models. The grids mentioned above thus do not cover this range and extrapolation in g had to be used. The rather low dependence of T 10 (T int , g) on gravity justifies this procedure.

In order to retrieve the behavior of T 10 (T int , g) in both the high and low T int regime (arbitrarily separated at T mid ), the numerical grids where fitted separately in the two regimes by functions of the form

A smooth function is recovered on the whole temperature range by interpolating linearly between the two functions in the range [T mid -∆T, T mid + ∆T ]. The 10 bar temperatures derived by this procedure do not differ by more than 2-3% with respect to the tabulated values. Values used for the fitting parameters are given in Table 1.