Layered convection as the origin of
Saturn’s luminosity anomaly
As they keep cooling and contracting, Solar System giant planets radiate more energy than they receive from the Sun. Applying the first and second principles of thermodynamics, one can determine their cooling rate, luminosity, and temperature at a given age. Measurements of Saturn’s infrared intrinsic luminosity, however, reveal that this planet is significantly brighter than predicted for its age. This excess luminosity is usually attributed to the immiscibility of helium in the hydrogenrich envelope, leading to "rains" of heliumrich droplets. Existing evolution calculations, however, suggest that the energy released by this sedimentation process may not be sufficient to resolve the puzzle. 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.
Many arguments suggest the existence of compositional gradients in giant planet interiors, as a consequence of either their formation process or their cooling history. 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 heat and atomic concentration can trigger a hydrodynamical instability, called the double diffusive instability, 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.
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. 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). 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. Extending the usual mixing length formalism 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, , or equivalently the corresponding dimensionless mixing length parameter, , where 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 for Saturn and for Jupiter. 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 cannot be used. Instead, the evolution is computed by integrating the intrinsic luminosity, , where the total energy includes the gravitational, thermal and rotational energies (see Methods). The size of the convective/diffusive layers is assumed to have reached an equilibrium value 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 (Supplementary Information).
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 23 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.
Fig. 2 suggests that Saturn models with convective/diffusive layers smaller than 23 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 , 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 (gastodust 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. Whereas these (unknown) initial conditions become inconsequential after about a few tens of million years for an adiabatic planet of Saturn’s mass, 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, socalled KelvinHelmholtz timescale (). 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 tradeoff between the initial energy content of the nascent planet (at the end of the runaway accretion shock) and the thickness or number of convectivediffusive 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. 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 (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 ) within less than 4.6 Gyr. Layered convection in this planet could then have stopped after the erosion of the whole core, or at least of its soluble components. This would be consistent with both i) the large initial core mass needed to form the giant planet on the right timescale and ii) Jupiter’s present smaller core. 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. Indeed, the strong compositional gradients (discontinuities) due to the presence 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, 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 that cannot be fully explained by the various proposed heating mechanisms  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, which is briefly outlined below.
The fluid is assumed to be composed of a large number of small convective layers of height , 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 (). The temperature gradient (), hence the superadiabaticity in the gas, , is given as a function of the flux to be transported, of the thermodynamic properties of the medium and of .
In the diffusive interfaces, the thermal gradient is equal to the gradient needed to transport the whole intrinsic flux by diffusive processes alone (). The thermal size of these diffusive interfaces, , can be estimated by the fact that the convective timescale, given by the inverse of the BruntVäisälä frequency () must be equal to the diffusive timescale in the interface (, where is the thermal diffusivity of the medium). Thus, , and the only remaining free parameter is the height of the convective layers, or equivalently the mixing length parameter, , where is the pressure scale height in the fluid.
Finally, the mean temperature gradient 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 (i.e. on ).
Planetary model sequences and time integration
Thermal evolution tracks of a planet of a given mass , composition (symbolically denoted by an array, , encompassing all the information needed to know the abundances of all the materials at all depths, including the core), angular momentum (), and thermal structure (parametrized here by the mixing length parameter ) yield relations of the type
(1) 
where and are the intrinsic and total effective temperatures, the StefanBoltzmann constant and 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. 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 () is the temperature that the planet would have if its intrinsic luminosity were zero and represents the contribution to the solar flux through
(2) 
where is the solar luminosity, the orbital distance of the planet, and its Bond albedo. Throughout this study, the values used for the equilibrium temperature are K and K for Jupiter and Saturn respectively. This corresponds to an absorbed luminosity of 5.14W (1.11W) 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 (Planetary model sequences and time integration) 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), , as a function of the intrinsic effective luminosity (), and of the gravity, . 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, is known, and the 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,
(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 (). 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 , and time dependence can be retrieved by choosing an initial value for the luminosity or initial energy content and by integrating
(4) 
References and Notes

Pollack, J. et al. A calculation of Saturn’s gravitational contraction history. Icarus. 30, 111128 (1977).

Fortney, J. J., Ikoma, M. Nettelmann, N., Guillot, T. & Marley, M. S. Selfconsistent Model Atmospheres and the Cooling of the Solar System’s Giant Planets. Astrophys. J. 729, 32 (2011).

Salpeter, E. E. On Convection and Gravitational Layering in Jupiter and in Stars of Low Mass. Astrophys. J. 181, L83L86 (1973).

Stevenson, D. J. & Salpeter, E. E. The dynamics and helium distribution in hydrogenhelium fluid planets. Astrophysical Journal Supplement Series. 35, 239261 (1977b).

Hubbard, W. B. & DeWitt, H. E. Statistical mechanics of light elements at high pressure. VII  A perturbative free energy for arbitrary mixtures of H and He. Astrophys. J. 290, 388393 (1985)

Pfaffenzeller, O., Hohl, D. & Ballone, P. Miscibility of Hydrogen and Helium under Astrophysical Conditions. Phys. Rev. Lett. 74, 25992602 (1995).

Lorenzen, W., Holst, B. & Redmer, R. Demixing of Hydrogen and Helium at Megabar Pressures. Phys. Rev. Lett. 102, 115701 (2009).

Morales, M. A. et al. Phase separation in hydrogenhelium mixtures at Mbar pressures. PNAS. 106, 1324 (2009).

Fortney, J. J. & Hubbard, W. B. Phase separation in giant planets: inhomogeneous evolution of Saturn. Icarus. 35 (2003).

Stevenson, D. J. Cosmochemistry and structure of the giant planets and their satellites. Icarus. 62, 415 (1985).

Podolak, M., Hubbard, W. B. & Stevenson, D. J. Model of Uranus’ interior and magnetic field. In Uranus, University of Arizona Press, 2961, (1991).

Chabrier, G. & Baraffe, I. Heat Transport in Giant (Exo)planets: A New Perspective. Astrophys. J. Lett. 667, L81L84 (2007).

Leconte, J. & Chabrier, G. A new vision of giant planet interiors: Impact of double diffusive convection. Astron. & Astroph. 540, A20 (2012).

Stern, M. E. The saltfountain and thermohaline convection. Tellus. 12, 172+ (1960).

Radko, T. What determines the thickness of layers in a thermohaline staircase? J. Fluid Mech. 523, 7998 (2005).

Rosenblum, E. P., Garaud, P., Traxler, A. & Stellmach, S. Turbulent Mixing and Layer Formation in Doublediffusive Convection: Threedimensional Numerical Simulations and Theory. Astrophys. J. 731,66+ (2011).

Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L. & Wood, T. S. A New Model for Mixing by Doublediffusive Convection (Semiconvection). I. The Conditions for Layer Formation. Astrophys. J. 750, 61 (2012).

Hansen, C. J. & Kawaler, S. D. Stellar Interiors. Physical Principles, Structure, and Evolution. SpringerVerlag Berlin Heidelberg New York (1994).

Pollack, J. B. et al. Formation of the Giant Planets by Concurrent Accretion of Solids and Gas. Icarus. 124, 62 (1996).

Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P. & Lissauer, J. J. On the Luminosity of Young Jupiters. Astrophys. J. 655, 541549 (2007).

Guillot, T., Saumon, D. & Stevenson, D. J. The Interior of Jupiter. In Jupiter. The Planet, Satellites and Magnetosphere.

Wilson, H. F. & Militzer, B. Solubility of Water Ice in Metallic Hydrogen: Consequences for Core Erosion in Gas Giant Planets. Astrophys. J. 745, 54 (2012).

Wilson, H. F. & Militzer, B. Rocky Core Solubility in Jupiter and Giant Exoplanets. Phys. Rev. Lett. 108, 1101 (2012)

Stevenson, D. J. Formation of the giant planets. Pla. Space Sci., 30, 755764 (1982).

Saumon, D. & Guillot, T. Shock Compression of Deuterium and the Interiors of Jupiter and Saturn. Astrophys. J. 609, 11701180 (2004)

Stevenson, D. J. & Salpeter, E. E. The phase diagram and transport properties for hydrogenhelium fluid planets. Astrophys. J. Suppl. Ser. 35, 221237 (1977a).

Miller, N. & Fortney, J. J. The Heavyelement Masses of Extrasolar Giant Planets, Revealed. Astrophys. J. Lett. 736, L29 (2011)

Laughlin, G., Crismani & M., Adams, F. C. On the Anomalous Radii of the Transiting Extrasolar Planets. Astrophys. J. Lett. 729, L7+ (2011)

Leconte, J., Chabrier, G., Baraffe & I., Levrard, B. Is tidal heating sufficient to explain bloated exoplanets? Consistent calculations accounting for finite initial eccentricity. Astron. & Astroph. 516, A64+ (2010).
Correspondence and requests for materials should be addressed to J.L.
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/20072013 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 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 (see Figure 3 of Leconte and Chabrier).
The fact that measured gravitational moments constrain to be larger than () 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 (see Fig. 3), also seems to justify this hypothesis.
2 Atmospheric model grids
The relations described by (Planetary model sequences and time integration) 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), , as a function of the intrinsic flux exiting the planet (parametrized by the intrinsic effective luminosity, ), and of the gravity, .
We use the recent atmospheric grids derived specifically for Jupiter and Saturn by Fortney et al.. 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 , low 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 had to be used. The rather low dependence of on gravity justifies this procedure.
In order to retrieve the behavior of in both the high and low regime (arbitrarily separated at ), the numerical grids where fitted separately in the two regimes by functions of the form
(5) 
A smooth function is recovered on the whole temperature range by interpolating linearly between the two functions in the range . The 10 bar temperatures derived by this procedure do not differ by more than 23% with respect to the tabulated values. Values used for the fitting parameters are given in Table 1.
3 Planetary structure integration
When is known, the profile is integrated inward by integrating the standard structure equations for a rotating body
(6)  
(7)  
(8) 
where is the mass enclosed in the equipotential of mean radius , the rotation rate of the planet, a secondorder correction due to the centrifugal potential. As discussed above, is either equal to in the adiabatic case, or to which depends on in a layered convection zone. The angular velocity is given by the fact that , 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 socalled ideal volume law for the mixture:
(9) 
where , and 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. For adiabatic, homogeneous models, the core masses and uniform metal enrichment in the envelope are equal to and in Jupiter and and in Saturn. For models of Saturn with layered convection, the metal enrichments used for a given mixing length parameter are shown in and the core mass used is . 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: models cool to the observed effective temperature in Gyr for Saturn, and 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, and that we take into account the thermal contribution of this core.