© ESO, 2012

1. Introduction

More than 500 planets have now been discovered orbiting stars outside our solar system, ranging from a few Earth masses to several Jupiter masses. Planets thus seem to be ubiquitous in nature. These discoveries raise fundamental questions about the inner composition, evolution, and origin of these bodies and about the fundamental properties of planets in general, including the ones in our own solar system. Characterizing their internal structure and composition, and from there, better understanding planet formation are major challenges in modern astronomy. Determining the heavy element content, for instance, provides key constraints to planet formation models, in particular on the efficiency of solid planetesimal accretion in the protoplanetary disk to build a planet embryo. While only the mass can be derived from observations for most detected planets, about 20% of these objects transit their parent star, allowing us to determine their radius. This gives us the mean density of the planet, constraining its gross composition (see Baraffe et al. 2010 for a recent review). Although providing an important diagnostic, this information is too limited to determine the element distribution within the planet, hence its precise compositional and thermal structure. Assuming planet formation is a universal process, one thus must turn to our solar system planets, in particular to the two gas giants Jupiter and Saturn, which encompass 92% of the planetary mass of the solar system, to derive more detailed interior models. Indeed, for our own giants, the gravitational moments have been determined with high accuracy by the various flyby missions entering their atmosphere and provide stringent constraints on their inner element distribution.

An important question, for instance, is whether the heavy elements present in giant planet interiors are located in a central core or are largely mixed into the hydrogen-helium (H/He) fluid envelope. In this case, a major issue is to determine whether convective mixing is efficient enough to yield a homogeneously mixed envelope or, alternatively, if giant planet interiors can exhibit a continuous compositional gradient. Exploring such a possibility is crucial to determining (i) the maximum amount of heavy elements compatible with observational constraints, and (ii) the efficiency of heat transport in giant planets. These two issues directly affect the planet mechanical (density), chemical (composition), and thermal (luminosity, temperature) structures at a given age, with major consequences on our understanding of planet formation and evolution.

Traditionally, giant planet models have always been based on two major specific assumptions about their inner profile, essentially for reasons of simplicity in the planet’s modeling (see e.g. Stevenson 1985). It is conventional to assume

• (i)

  that the inner structure of our giants consists of a few – generally 2 to 3 – superposed, well separated, homogeneous regions; that is, going from the planet’s center to the surface, a central solid rocky/icy core and a surrounding largely dominantly H/He gaseous envelope often split into an inner metallic region and an outer atomic/molecular one. The denser components are always supposed to have been accreted first or to have quickly settled into the center under the action of gravity (for planetesimals accreted after the runaway gas accretion);

• (ii)

  that planetary interiors are adiabatic, based on the need for giant planet heat flow to be transported by convection (Hubbard 1968).

All the present determinations of the internal – chemical, mechanical, and thermal – structures of the solar system planets, including their heavy material content, assume homogeneously stratified, adiabatic interiors (Saumon & Guillot 2004; Fortney & Nettelmann 2010).

Giant planet interiors, however, might depart from this conventional, simplified description, because of complex processes for which we lack an accurate description but which may very well be at play in real situations (Stevenson 1985). In this paper, we derive interior models for Jupiter and Saturn that relax the aforementioned preconceptions. Instead of the homogeneous layer assumption, we explore the possibility of a mixed, inhomogeneous solid-gas interior composition, leading to a heavy material gradient throughout the planet. Owing to the double diffusive instability (see Sect. 2.2), this gradient tends to suppress large-scale convection which can turn into either turbulent enhanced diffusion or layered convection. Because both these heat transport mechanisms are fairly inefficient compared to the usual convection, this compositional gradient thus leads to a significant departure from adiabaticity in the interior. As shown below, these models do fulfill the planet’s observational constraints, while leading to (i) a significantly increased metal content and (ii) significantly higher internal temperatures than inferred from homogeneously stratified adiabatic models.

This opens a new vision on planet structure, evolution, and formation efficiency. Such inhomogeneous interior profiles for solar system giant planets had briefly been suggested several decades ago by Stevenson (1985), but no attempt has ever been made since to derive consistent models and to verify whether such models would be consistent with the planet’s various observational constraints. This scenario has been revived recently in the context of extrasolar planets and has been shown to provide a possible or at least complementary explanation for the anomalously large observed radii of many of these bodies (Chabrier & Baraffe 2007). Indeed, as mentioned above, not only does an inhomogeneously stratified interior yield a different interior structure and global metal content, but it also decreases heat transport efficiency throughout the planet’s interior and thus affects its cooling, hence its mass-radius relationship at a given age, a crucial diagnostic in understanding (transiting) extrasolar planet structure and evolution. Since, as mentioned earlier, only loose constraints on the object’s internal composition are accessible for gaseous exoplanets, it is crucial to verify whether such unconventional internal structures are a viable possibility for our own giants. Furthermore, determining the maximum possible amount of heavy elements in Jupiter and Saturn and their distribution within the planet are important diagnostics of how our own solar system giants formed.

To address this issue, we first briefly review our current understanding of the double diffusive instability, and of the various regimes under which it can occur in Sect. 2. Then, in Sect. 3, we derive an analytical formalism, based on a standard parametrization of convection, similar to the mixing length formalism, which describes the global transport properties of an inhomogeneous convective/diffusive medium. In Sect. 4, we discuss the possibility of an equilibrium size for the convective/diffusive layers present in a semiconvective planet and derive analytical constraints for this equilibrium value. Finally, in Sect. 5, we derive semiconvective models of Jupiter and Saturn that are consistent with the gravitational moments and with the surface abundances measured by the Galileo and Cassini missions. This enables us to further constrain the number of possible convective/diffusive cells within these planets and to derive new values of the heavy element content and of the core masses for our gas giants.

2. Physical mechanisms leading to inhomogeneous density stratification

2.1. Compositional gradient

In the present calculations, in contrast to all previous planetary models, we consider a mass fraction of heavy material Z(m) at a depth r(m) within the planet (i.e. at the depth of the isodensity surface enclosing a mass m of the planet) that is continuously decreasing from the core to the surface, producing a compositional gradient within the gaseous envelope of the planet, (1)where is the pressure scale height.

The possible origin of such a compositional gradient is an open question (Stevenson 1985; Chabrier & Baraffe 2007). In the conventional scenario, all the accreted planetesimals are assumed to directly sink to the core and not to evaporate in the envelope, for sake of simplicity. In reality, however, incomplete mixing of large planetesimals or dissolution of a substantial fraction of volatiles and rocks from small solid bodies could occur in the envelope during the phase of planetesimal accretion on the nascent planet; a substantial amount of ice could also remain in the envelope (Iaroslavitz & Podolak 2007; Hori & Ikoma 2011).

The gradient might also stem from an only partial redistribution by small-scale convective motions of stably layered (soluble) constituents released by core erosion in the gas-rich envelope during the planet’s evolution (Stevenson 1982; Guillot et al. 2004), as seems to be supported for water by recent numerical simulations (Wilson & Militzer 2012). This could be enhanced by the immiscibility (phase separation) of an abundant enough material (e.g. helium, water) in the dominantly metallic-hydrogen envelope (Stevenson & Salpeter 1977a). This would change the dynamical properties of (double diffusive) convection near the regions of immiscibility and make the problem more complex. Therefore, to avoid extra complication, we only consider the occurence of this process for soluble material in the present paper.

At last, rapid rotation and/or strong magnetic fields, necessarily present in Jupiter and Saturn interiors , are known to hamper large-scale convection (Chabrier et al. 2007), possibly leading to imperfect mixing of heavy elements in part of the envelope.

2.2. Double diffusive convection

Convective systems in which (rapidly diffusive) heat and (slowly diffusive) composition have opposite destabilizing and/or stabilizing effects tend to develop inhomogeneous density profiles. The (de)stabilizing effect of heat can be quantified by the thermal gradient, , or more precisely by the super adiabaticity (∇_T − ∇_ad, where the derivative in ∇_ad is taken at constant specific entropy), and the one of the heavy elements by the gradient of mean molecular weight (μ), . The actual dynamical state of the medium depends on the value of the density ratio (2)where , ¹ (Stern 1960).

When both the mean molecular weight (due e.g. to a higher concentration of salt in salty water) and the temperature increase with height (∇_μ and (∇_T − ∇_ad) < 0), the compositional gradient is destabilizing, while the temperature gradient is stabilizing. This is referred to as the fingering case. In that case, convective instability develops when R_ρ < 1, which is equivalent to the Ledoux instability criterion. But even if 1 < R_ρ < 1/τ, where τ = D/κ_T is the ratio of solute (D) to thermal (κ_T) diffusivities, the slower diffusivity of elements compared to heat yields the so-called double-diffusive instability, which in turn leads to the formation of salt fingers and, sometimes, of thermo-compositional staircases, as observed in some parts of the oceans and in laboratory experiments (Turner 1967).

The opposite case, referred to as the diffusive case, corresponds to a fluid exhibiting a positive molecular weight gradient (∇_μ > 0). In that case, the fluid will be convectively unstable if this gradient is insufficient to stabilize the system against convective instability, i.e. if the Ledoux instability criterion, (3)is met². However, in the regions which are stable according to the Ledoux criterion, i.e. , but unstable according to the Schwarzschild criterion, (4)what happens is less clear, especially in the astrophysical context, where the very low values of the Prandlt number (Pr ≡ ν/κ_T, where ν is the kinematic viscosity) make direct numerical hydrodynamical simulations difficult. In particular, the exact nature of double diffusive convection, if it occurs, remains uncertain (homogeneous oscillatory convection or layered convection, i.e. uniformly mixed convective layers separated by thin diffusive interfaces characterized by a steep jump in the mean molecular weight). Analytical arguments (Radko 2003) and recent 3D hydrodynamical simulations (Rosenblum et al. 2011; Mirouh et al. 2012), however, seem to suggest the picture presented in Fig. 1. When the mean molecular weight gradient () decreases in a stable medium, homogeneous oscillatory convection, also called turbulent diffusion, first appears for (5)while well defined thermo/compositional layers start to develop when (6)Here corresponds to the point where the solute to heat buoyancy flux ratio ( ≡ γ^-1) stops decreasing when increases (see Radko 2003 for details). Its exact value, however, depends on the characteristics of the medium in a nontrivial way and is difficult to estimate (Rosenblum et al. 2011; Mirouh et al. 2012). For smaller ∇_μ gradients, the medium is unstable according to the Ledoux criterion, and the thermal forcing is strong enough to force large-scale overturning convection.

Various arguments seem to support, or at least not to exclude, the existence of layered convection under planetary conditions (Chabrier & Baraffe 2007). Conducting 3D hydrodynamics calculations over a wide domain of parameter space (Prandlt number and atomic to thermal diffusivity ratio) including the regime relevant for planetary interiors, Mirouh et al. (2012) always find a domain where γ^-1 decreases with , a necessary and sufficient condition for the layering instability, thus layer formation to occur. A central question is then the size of the layers, which is examined in Sect. 4. In any event, both homogeneous double-diffusive convection and layered convection – generically denominated as “semiconvection” in the following – are found to yield thermal and compositional fluxes that are significantly smaller than expected from standard convection. Indeed, the presence of diffusive interfaces strongly decreases the efficiency of heat transport compared with large-scale, adiabatic convection, leading in planet interiors to a significant departure from the usual adiabatic profile, as quantified below.

3. An analytical theory for layered convection

To investigate the impact of such strongly hampered convection on giant planet internal structure, we developed a simple subgrid model of layered convection. As illustrated in Fig. 2 and found in simulations (Rosenblum et al. 2011; Mirouh et al. 2012), we consider that a semiconvective zone consists of a large number, N_l, of well mixed convectively unstable layers of size l, separated by thin diffusive interfaces of thermal thickness δ_T, within which the large stabilizing compositional gradient completely inhibits convective motions.

3.1. Convective layers

Within each convective layer, the fluid is expected to follow the dynamics found in turbulent Rayleigh-Bénard convection within a cell of typical height equal to the size of the layer, l. By analogy with the mixing length formalism (see details in Appendix A), we define a dimensionless mixing length parameter by dividing l by the pressure scale height, α ≡ l/H_P.

In a laboratory or a numerical experiment, the efficiency of the convection is characterized by the thermal Nusselt number, (7)where by definition the total intrinsic flux (F_tot), the flux transported by diffusive processes (F_d), and the diffusive flux that would be present in a completely adiabatic zone () are given by (Cox & Giuli 1968; Hansen & Kawaler 1994) (8)where c_P is the heat capacity at constant pressure. It is found that, for high Rayleigh numbers, the Nusselt number follows a law of the type (9)where Ra _⋆  is a modified Rayleigh number³, which is the ratio of the strength of the thermal forcing to the one of the radiative losses (10)Here, g is the local gravity acceleration, and N_T the Brunt-Väisälä frequency.

As convection at very high Rayleigh numbers is difficult to study either experimentally or computationally, it is difficult to give precise values for a and C_L. For the bounded Rayleigh-Bénard problem, theoretical models suggest that the exponent of the convective flux law, a, could be equal to one third (Garaud et al. 2010)⁴. However, experiments done by Krishnamurti (1995) tend to show that this exponent could be lower, and as low as a = 0.2. On the other hand, for homogeneous Rayleigh-Bénard convection (without boundaries), Garaud et al. (2010) show that the regime predicted by the mixing length theory, i.e. (C_L = 1, a = 1/2; see Appendix A), is recovered. In the following, we thus consider and C_L = 1.

Equation (9) is sufficient to calculate the flux transported by convection once the super adiabaticity is known. To compute this, however, we must first define a quantity that can be computed a priori from the local thermodynamical properties of the medium and the total internal energy flux to be transported. Following Hansen & Kawaler (1994), this convective forcing can be defined by (11)Introducing ε_d ≡ ∇_d − ∇_ad, we rewrite Eq. (11) as Φ ≡ Φ₀   α⁴   ε_d where (12)It is clear from Eq. (11) that Φ₀ is a local constant of the medium, which characterizes its ability to transport energy by convection, independently of the mixing length or of the flux to be transported ( ∝ ε_d).

Then, from Eqs. (9) and (11), one sees that in a region where convection remains efficient enough, (13)which yields the super adiabaticity, (14)The range of super adiabaticity in the convective layers implied by this equation for the various possible exponents a is shown in Fig. 3 (pale red area). As seen, the uncertainty on a leads to a large dispersion on this super adiabaticity. In this high convective efficiency regime, we can further compute the mean convective flux which, by definition, is given by (15)As mentioned above, this scaling law only applies to the vigorous convection regime, i.e. Φ or Ra _⋆  ≫ 1. Thus, in terms of the layer height, convection remains efficient as long as α ≫ α_crit, where α_crit is defined such that Ra_⋆(α_crit) = 1, and, from Eq. (10), the critical layer size is (16)where d is the lengthscale of the fastest growing mode of the linear instability, (e.g. Baines & Gill 1969). Analytical and numerical arguments show that the size of the fastest growing layers is equal to 10–100   d ≳ l_crit (Radko 2003; Rosenblum et al. 2011). The efficient convection regime is thus appropriate in the planetary domain, where Pr ~ 10^-2 − 10^-1 (Chabrier & Baraffe 2007).

For conditions prevailing in the interior of the actual Jupiter (Saturn), the mean thermal diffusivity is κ_T ~ 5 × 10^-5 m² s^-1 (Potekhin 1999), Φ₀ is equal to Φ₀ ≈ 3 × 10³³ (9 × 10³²) and ε_d ≈  10 (10; see also Fig. 5), so that α_crit = 2 × 10^-9 (3 × 10^-9). In the following, all the order-of-magnitude estimates done throughout the text use these values.

3.2. Interfaces

In the interfaces of thermal size δ_T, overturning convection is inhibited by the strong jump in molecular weight. However, these regions do not need to be in the fully diffusive regime (), but can also be in the oscillatory convection regime characterized by Eq. (5). Indeed, as shown by Radko (2005) the condition , where is the inverse density ratio within the interface, is a sufficient criterion not only to ensure the stability of the interface itself, but also of the whole stack of layers that would otherwise merge into one homogeneous layer (see Sect. 4.1 for details).

Therefore, the interface is very likely either in a stable diffusion state or in a state of weakly turbulence enhanced diffusion. Rosenblum et al. (2011) show that in this regime Nu_T ≲ 2, meaning that most of the energy is transported by pure diffusion. We thus assume that the thermal gradient to be used in these regions is given by the gradient needed to transport the whole outgoing energy flux by diffusion: (17)The thermal diffusivity κ_T encompasses the contribution of all diffusive processes. If diffusion is ensured by photons, as generally in most astrophysical objects, the diffusive thermal gradient reduces to the so-called radiative gradient, ∇_rad (e.g. Hansen & Kawaler 1994). However, in the deep interior of giant planets and in degenerate bodies, density can be high enough for the electrons to become degenerate enough to efficiently conduct thermal energy (see Stevenson & Salpeter (1977b) and Chabrier & Baraffe (2007) for the characteristic radiative and conductive opacities under Jovian planet conditions). Along this paper, we use the generic denomination ∇_d for the diffusive temperature gradient, keeping in mind that diffusion is now due to electronic or atomic motions, with a characteristic thermal diffusivity κ_T. In the present calculations, we use the conductive thermal diffusivities calculated by Potekhin (1999).

3.3. Mean thermal gradient

Once we have calculated the thermal gradient in the convective zones of size l, and in the diffusive interfaces of size δ_T, we need to determine the mean properties of a whole stack of convective-diffusive cells. As discussed in Chabrier & Baraffe (2007), since the convective plumes must be fed by the diffusive interfaces, the thermal convective () and diffusive () time scales should be similar in each respective layer. Therefore, (18)Not surprisingly, the size of the interface, δ_T, is related to the lengthscale of the most unstable mode of the linear instability by δ_T =    d/Pr^1/4. In addition, comparing the mean kinetic energy of an upwelling eddy with the potential energy barrier created by the negative buoyancy in the diffusive interface, we see that the above condition also entails that convective overshooting can be neglected (see Chabrier & Baraffe 2007 for details).

Knowing δ_T/l enables us to compute the mean thermal gradient to be used to compute the planet’s structure (see Sect. 5.1) (19)Substituting ∇_T and δ_T/l by their expressions in Eqs. (14) and (18), and developing the mean gradient to first nonvanishing order, we find that (20)or equivalently (21)As expected, two terms appear in the mean super adiabaticity, which is shown in Fig. 3 (dark gray region). The first one is due to the temperature jumps at each interface and is controlled by the convective overturning timescale that determines the size of these interfaces. The second one is due to the super adiabaticity in the convective layers themselves, which is enhanced by their small size. The relative contributions of these two terms only depend on the properties of the convection, characterized by the exponent a.

Interestingly, if the impact of a on the super adiabaticity in the convective layers is significant, the impact on the mean super adiabaticity of the whole stack of layers,  ⟨ ∇_T ⟩  − ∇_ad, is rather small (see Fig. 3). This can be understood as follows. For the smallest value of the exponent considered here, a = 0.2, convection is very inefficient and most of the mean super adiabaticity is contained in the convective layers (upper contour of gray region). In that case, Eq. (19) is dominated by the second term,  ⟨ ∇_T ⟩  ≈ ∇_T.

When the exponent is increased, the super adiabaticity in the convective layers (∝Φ ^{− a/(1 + a)}) drops dramatically, but this entails a growth of the convective timescale, and thus of the thickness of the interfaces (see Eq. (18)). For efficient convection, , this effect eventually dominates, and the mean super adiabaticity again increases with a. This yields a rather narrow region for the mean super-adiabaticity within the considered range of exponent variation.

Interestingly, to first order, the a = 1/2 and a = 0.2 cases, i.e. the two extreme values considered here, yield the same mean superadiabaticity. Therefore, the model presented here for semiconvection depends only weakly on the precise prescription chosen for convection. In the following, we show results using C_L = 1 and a = 1/2, while keeping in mind that for α ≫ α_crit these results would be the same as the one obtained in the a = 0.2 case.

3.4. Mean solute gradient and flux

Disregarding convective overshoot, the solute flux (F_Z) is determined by the transport at each interface. As for heat, the turbulent transport of solute is expected to be very inefficient, and we assume purely diffusive processes (see Sect. 3.2). The impact of a more efficient transport is discussed in the next section. In this regime, the solute flux is thus given by (22)where the subscript I describes quantities determined within the interface, indicating that the solute flux F_Z is determined by the element-fraction variation rate at each interface. When neglecting the small inhomogeneities in the convective layer, the interfacial and mean gradients are related by (23)where  ⟨ ∇Z ⟩  describes the mean gradient of the heavy-element mass fraction, i.e. the value of dZ/dr averaged over several convective/diffusive cells, and δ_Z is the length over which compositional jump occurs. To constrain the value of δ_Z, two limiting arguments can be used. On one side, equating the diffusive timescales for heat and solute leads to (24)On the other hand, considering that the compositional jump must cover the whole interface to ensure its mechanical stability, one would expect that (25)The present paper focuses on the impact of double-difffusive convection on the planet’s mechanical and thermal structure. The impact on the evolution will be addressed in a forthcoming paper (see Sect. 6). Therefore, this uncertainty on the heavy element flux, discussed in the next section, does not have a significant impact on the results presented in the present study. Determining the solute transport properties in the regime of layered convection more precisely, however, will be central to evolutionary calculations. 3D hydrodynamical simulations in a realistic parameter range are thus strongly needed (Mirouh et al. 2012).

The term ∇_μ,I, which is the relevant quantity when computing the equilibrium condition of the interface, is given by (26)The precise value of depends on the precise chemical composition of the considered heavy element, but is typically around unity, the value we use in numerical estimates below. We also use D ≈ 10^-7 − 10^-8 m² s^-1 (i.e. τ ≈ 10^-2 − 10^-1), appropriate for giant planet interior conditions (Stevenson & Salpeter 1977b).

4. Layer size: analytical point of view

4.1. Existence of an equilibrium height

The determination of the size of the layers, when layered convection is occurring, is a complex task. The problem is rendered even more difficult by the tendency for the small layers that initially form to merge into thicker layers. The question is then to know if and when layer merging ends.

In the fingering case, i.e. when the solute is the destabilizing actor, Radko (2005) derived a criterion for the linear stability of thermo-compositional staircases against merging. In this picture, an initially inhomogeneous medium starts forming layers if its density ratio (R_ρ) is below the density ratio, R_min, for which the heat to solute buoyancy flux ratio ( ≡ γ) stops decreasing when R_ρ increases. Then, layers are unstable and merge as long as the density ratio within the interface between them (R_I) is less than R_min. Then, as R_I increases with the layer height (under some conditions), it eventually reaches R_min and the merging process stops. At this stage, an equilibrium height is reached by the convective layers and the staircase is fully equilibrated.

The linear stability analysis of Radko (2005) can be redone in the diffusive case, by simply accounting for the need for the signs of the various density and temperature differences and of the fluxes to be changed. Then, after some lengthy but straightforward algebra, one can show that layers are unstable and merge when (27)which defines . As for the fingering case, this linear stability analysis shows that, if layers form, their merging will stop when their height reaches a finite equilibrium value.

4.2. Global constraints

Although precisely estimating this equilibrium value is difficult, some strong limits on the layer size can be derived theoretically. We consider a stack of layers extending over a zone of size L, and define a global gradient  ⟨ ∇Z ⟩  ≈  −ΔZ/L, where ΔZ is the difference between the mass ratios of heavy element at the bottom and at the top of the semiconvective zone.

On one hand, the mean molecular weight gradient in all the interfaces, ∇_μ,I, must be high enough to satisfy the stability criterion discussed in Sects. 2.2 and 3.2. This implies that (28)Substituting ∇_μ,I by its expression in Eq. (26) and δ_Z by using Eqs. (24) and (18), the criterion (28) for a planetary scale (L ≈ R _p) semiconvective zone reads as (29)where α_min denotes a lower limit for the layer size. Under the present conditions in the interiors of our gas giants (see end of Sect. 3.1 and Table 1), this yields , with the less restrictive constraint being obtained for , as summarized in Table 1. However, as the existence of layers allows to be as high as (1 + Pr)/(τ + Pr) ≈ 10¹−10², these constraints could be severely tightened, as shown in Table 1. This simply confirms that, in order for layered convection to be stable, convective cells must remain larger than the diffusive interfaces, and the medium is always in the convective regime, α ≫ α_crit.

On the other hand, the solute gradient within the planet will be homogenized on a typical timescale (30)Using Eq. (22)–(24), and taking  ⟨ ∇Z ⟩  ≈ ΔZ/R _p, Eq. (30) becomes (31)Therefore, to avoid a complete homogenization of giant planet interiors in less than 5 Gyr, their present age, α must be less than α_max ≈ 1 × 10^-4 for a = 0.2 and 1.4 × 10^-3 for a = 0.5. However, as mentioned in Rosenblum et al. (2011), turbulent transport due to the double-diffusive instability can yield compositional Nusselt numbers, Nu_μ⁵, around 2 to 4, yielding even stronger constraints, as summarized in Table 1. These low values of α justify a posteriori the approximation of continuous thermal and heavy element profiles when considering the planet’s entire internal structure. However, layered inhomogeneities could be dynamically regenerated over time. In that case, layered convection will be an ongoing process in the planet’s interior.

Considering the possibility that the ratio of the compositional to the thermal size of the interface does not scale as the square root of the ratio of the diffusivities, so that Eq. (25) should be used instead of Eq. (24), yields an increase in both α_min and α_max by a factor .

This analysis shows that, for the age of the solar system and for the conditions prevailing in gas giant interiors, a range of layer sizes exists for which ongoing layered convection is a viable mechanism. According to our estimate, this range is relatively large and corresponds to (32)The uncertainty on the lower and upper bounds are due to our poor knowledge of the behavior of , and of the solute flux at low Prandlt number, respectively. Even given these uncertainties, that convection is always in a regime of relatively high Rayleigh number appears to be a robust prediction. To make an attempt to overcome these limitations, we examine in Sect. 5 how observational data can narrow this possible domain of α, by further constraining the degree of super adiabaticity, and thus the size of the layers, in our solar system gas giants.

4.3. An alternative scenario

We stress that the aforementioned constraints only apply to layered convection and do not preclude the possibility that, under some conditions, double diffusive convection may manifest itself in the form of homogeneous double diffusive convection and act like a turbulent diffusion (Rosenblum et al. 2011). However, for this to happen, criterion (28) must be verified in the inefficient convection regime, δ_T ≫ l (see Appendix A). The following criterion must then hold (33)As and , by definition, a zone of turbulent diffusion cannot extend over the entire planet’s scale unless ∇_d ≲    ∇_ad, i.e. if the whole object is nearly diffusive in the first place.

In an object with a heat flux high enough to be convectively unstable, the semiconvective zone must then be confined to a fraction of the planet, in particular, but not necessarily, near an immiscibility region or a phase transition for instance. In this case, the total size of the zone must verify (~1/10 in Jupiter), a condition for which a large enough jump in the heavy element mass fraction can be enough to stabilize the whole zone against convection and to open a diffusive buffer in the interior (where ∇_T ≈ ∇_d). From the global point of view of the planet, this would nearly act as a discontinuity in composition, temperature, and entropy. This possibility that such a diffusive-like buffer exists in the interior of our gas giants, as a consequence of double-diffusive instability, must be kept in mind.

In the following, however, we do not consider this scenario any further and only consider the effect of a planetary scale layered-convection zone.

5. Numerical results for solar system giant planets

In this section, we examine whether the presence of semiconvection in Jupiter’s and Saturn’s interiors can be consistent with the various available observational constraints. We first derive homogeneous interior reference models in Sect. 5.2. Then, in Sect. 5.3, we incorporate our model for layered convection into the standard method used to compute interior structure models of rotating gaseous planets (presented in Sect. 5.1) and determine the area of the composition/layer size space parameter which is consistent with observed gravitational moments and surface abundances.

5.1. Hydrostatic equilibrium and figures of the planet

Solar system giant planets are rapidly rotating bodies (the period of rotation is about 10 h), with the centrifugal potential representing about 10% of the gravitational potential. This modifies the hydrostatic equilibrium condition between the pressure gradient and the gravitational force in the interior, which now writes as (34)where (35)and (36)denote respectively the gravitational and centrifugal potentials, with differential rotation   ω  _p(ξ), where ξ is the distance from the position r to the rotation axis, and G the gravitational constant. In the present study,   ω  _p is assumed to be constant and given by the magnetospheric rotation rate. Because of the symmetry of the centrifugal potential with respect to both the rotation axis and the equatorial plane, surfaces of equal densities for these objects are supposed to be generalized ellipsoids of revolution whose exact shape is given by (37)where is the mean radius of the equipotential, P_2n are the usual Legendre polynomials, θ is the colatitude, and the s_2n are a set of figure functions. The last can be derived using the theory of figures for rotating bodies (Zharkov & Trubitsyn 1978), and must be solved iteratively with the set of perturbed 1D hydrostatic equilibrium equations where m is the mass enclosed in the equipotential of mean radius , ω _p the rotation rate of the planet, a second-order correction due to the centrifugal potential, which depends on the figure functions. As discussed in Sect. 3, the prescription to be used for ∇_T is determined by the energy transport processes.

The departure from sphericity of the isodensity surfaces results in a perturbation of the external gravity field V_G(r,θ) that writes as where r is the radial distance from the center of the planet, M the mass of the planet, R_eq the equatorial radius, θ the colatitude, P_2n are Legendre polynomials of order 2n, and J_2n denote the gravitational moments, which can be computed once the figure equations have been solved. The measured gravity moments provide stringent constraints on the density profile and the possible layering within these planets.

As, in practice, Legendre polynomial expansions are truncated at a given order n, 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: (43)where X, Y and Z denote the mass fractions of H, He and heavy elements, respectively. For the H/He fluid, the most widely used EOS is the Saumon-Chabrier-vanHorn EOS (Saumon et al. 1995; SCvH). For the heavy material, we have used the “Rock” EOS of Hubbard & Marley (1989) for silicates and the “Ice” ANEOS equation of state (Thompson & Lauson 1972) for volatiles (CH₄, NH₃, H₂O). The impact of the differences between various EOS’s on exoplanet structure and evolution has been explored in Baraffe et al. (2008).

Once such EOSs, P [ρ(X_i)] , have been specified, structure models with various compositions are calculated by solving iteratively the aforementioned hydrostatic equilibrium condition for a rotating body and the third-order level-surface theory (Zharkov & Trubitsyn 1978) to obtain a model that reproduces the observed values of the radius, R_eq, and gravitational moments J₂ and J₄ measured by the Pioneer and Voyager missions (see Table 2).

5.2. Reference homogeneous models

In conventional giant planet models, the abundances of heavy elements are chosen to be constant in the gaseous H/He envelope, with a possible discontinuity at the transition between the molecular and metallic regions (Chabrier et al. 1992). Under the actual conditions found in Jupiter and Saturn, the thermal gradient that would be needed to transport the whole flux by diffusive processes, ∇_d, is always larger than ∇_ad and the whole interior is convective in the homogeneous case according to the Schwarzschild criterion. Since convection is very efficient (see Sect. 3), the super adiabaticity needed to transport the outgoing energy is in the range of 10^-8 − 10^-9, so that the structure can be solved by setting ∇_T = ∇_ad in Eq. (40).

To have a reference case, we use the formalism described in Sect. 5.1 to derive homogeneous, adiabatic interior models representative of the usual two-layer composition. Since we use the interpolated SCvH EOS, we do not consider the effect of a plasma phase transition, and we are thus left with only two free parameters, namely the core mass (M_c) and the metal mass fraction in the gaseous envelope (Z_env)⁶. The temperature, density, and pressure profiles of our most representative homogeneous models of Jupiter and Saturn are shown in Fig. 4 (solid curves). These are composed of a solid core of mass M_c = 3.9 and 25.6   M_⊕ surrounded by a H/He gaseous envelope with a constant metal fraction Z_env = 0.11 for Jupiter and Z_env = 0.05 for Saturn (summarized in Table 3). These reference models yield interior enrichment that is consistent with previous determinations (Chabrier et al. 1992; Saumon & Guillot 2004; Guillot 2005).

5.3. Inhomogeneous models

We now derive semiconvective, inhomogeneous interior models for Jupiter and Saturn. We stress that all these models are consistent, within the observational uncertainties, with the measured gravitational moments of Jupiter and Saturn(see Table 2; Campbell & Synnott 1985; Campbell & Anderson 1989).

An additional constraint on the outermost value of the compositional gradient is provided by the surface abundance of heavy elements in the planets measured by the 1995 Galileo Entry Probe mission. Indeed, elemental abundances of the atmospheres of solar giant planets are observed to differ significantly from each other and from the solar composition, because they are enriched by a factor  ~2–4 and  ~2–8 with respect to the Sun’s atmosphere for Jupiter and Saturn, respectively, as shown in Table 2 (Guillot 2005). Moreover, the planet’s total mean abundances of H and He ( and ) must recover the values of the protosolar nebula, i.e. .

In our calculations, the adjustable parameters to fulfill all these constraints are chosen to be the mass of the core (M_c), the mean heavy element mass fraction in the gaseous envelope (), and the global compositional variation in the envelope (ΔZ_env, the difference between the metal mass fraction just above the central core and the one in the atmosphere). To assess the robustness of our results with respect to the equation of state chosen to describe the thermodynamics of the heavy material, we derived several sets of models for which the composition of the core varies from pure ice to pure rock.

5.3.1. Number of layers

From a macroscopic point of view, an important quantity describing layered convection is the number of convective-diffusive layers, N_l. This number is roughly equal to the ratio of the size of the semiconvective zone, comparable to the planet’s radius, R _p, if this zone extends over the whole planet, to the height of a typical convective/diffusive cell, l + δ_T. As shown in Sect. 4, in the regime of interest, l + δ_T ≈ l, and N_l ≈ R _p/l. Because H_P ≈ R _p in the deep interior, the number of layers in the planet is thus approximately equal to N_l ~ α^-1, and in the following we always refer indifferently to either α or (44)As shown in Sect. 4, N_l is constrained to lie within the range 10^2 − 4 ≤ N_l ≤ 10⁶⁻⁹. Note that, given the small size of the diffusive-convective layers compared with the size of the planet, the discontinuous (staircase-like) temperature and composition profiles can be approximated well by continuous mean thermal and compositional gradients ( ⟨ ∇_T ⟩  and  ⟨ ∇_μ ⟩ , respectively) to determine the planet’s global structure, as illustrated in Fig. 2.

This possible range of numbers of layers is further constrained by our numerical calculations, which show that, in order to reproduce our giant planet observational constraints, no more than  ~2.5  ×  10⁵ layers (α_min ≈ 4 × 10^-6) can in reality be present in Saturn and  ~3  ×  10⁴ (α_min3 × 10^-5) in Jupiter (see Table 1). Indeed, more layers leads to such high temperatures in the interior that the induced mean density decrease cannot be counterbalanced by an increase in the heavy element mass fraction compatible with the observed surface abundances. This is because the larger the number of layers, the smaller the size of each convective cell, reducing the maximum height a convective eddy can travel to transport heat before being stopped by the negative buoyancy present in the diffusive interface.

A large number of layers thus decreases convective heat (and composition) transport efficiency. This leads to an increase in the mean super adiabaticity, as portrayed in Fig. 5, which in turn immediately implies a rise in the internal temperature, as illustrated in Fig. 4. It is important to stress that super adiabaticity is the physical quantity most directly constrained by the data. Thus, whereas the allowed range of number of layers (or equivalently of layer sizes) may slightly depend on the model used to parametrize semiconvection (see Sect. 3), the allowed range of super adiabaticity displayed in Fig. 5 should remain weakly affected.

5.3.2. Enhanced heavy material enrichment

The pressure, density, and thermal profiles obtained in the most extreme semiconvective case compatible with the observational constraints discussed above are shown in Fig. 4 (dashed curves). As seen in the figure and as expected from the above discussion, the non-adiabatic envelope profile obtained in the semiconvective case yields substantially higher internal temperatures than the usual adiabatic assumption, since heat and material redistributions are partly inhibited by diffusive processes. The pressure and density profiles, on the other hand, are barely affected, since they are strongly constrained by the gravitational moments. As a result, at basically fixed density profile, a higher temperature profile must be compensated for by a larger amount of heavy material within the envelope. This is illustrated in Fig. 6, where we show the abundance profiles, as calculated in Appendix B, corresponding to semiconvective models with different numbers of layers. The bottom curve in each panel corresponds to models with 1000 layers, while the other curves correspond to a gradually increasing number of layers.

Therefore, to compensate the radius increase (density decrease) due to the hotter interior, semiconvection yields a significantly larger total metal content compared with conventional homogeneous models. This can be seen in Fig. 7, which shows the amount of heavy elements in the core and envelope for the various cases discussed here, as summarized in Table 3. For Saturn, up to 50   M_⊕ of the heavy elements could be present in the planet, while for Jupiter the heavy material content could reach 63 M_⊕. This corresponds to about 25 and 10 times the solar abundances, respectively⁷. Since these values only depend on the allowed amount of super adiabaticity, they should not strongly depend on our modeling of diffusive/convective transport, as mentioned above. In contrast, the maximum amount of heavy elements compatible with the observational constraints for the homogeneous, adiabatic models, is about 30   M_⊕ for Saturn and 40 M_⊕ for Jupiter, in agreement with previous studies (Saumon & Guillot 2004).

Semiconvection not only increases the global metal content, but it also yields a completely different distribution of heavy elements. While the global enrichment of the planet is increased in the inhomogeneous models, the mass of the central core is decreased, as heavy elements are primarily redistributed in the gaseous envelope. In the case of Saturn, the vertical spread in core mass at fixed number of layers illustrated in Fig. 7 is obtained when varying the core composition from pure ice (top) to pure rock (bottom). In Jupiter the inferred core mass is too low for the equation of state to make a significant difference. One could wonder why the homogeneous case is not continuously recovered when α tends toward 1. This slightly counter-intuitive effect arises because completely homogeneous models (central core plus a fully homogeneous envelope) cannot in general reproduce both the observed J₂ and J₄, at least when using the SCvH EOS (Chabrier et al. 1992; Saumon & Guillot 2004). Thus, if we relax the constant Z condition in the envelope, the presence of a compositional gradient and of a smaller core appears to be the best solution for reproducing observational data, even in the absence of any additional superadiabaticity.

For Jupiter, models can be found that match the gravitational moments without the presence of a central, completely differentiated core. These are shown as red dots on the bottom right of Fig. 7). Such cases yield an atmospheric metallicity   Z  _atm ~ 4−5   Z_⊙. That the possible erosion of the core mass would have been more efficient in Jupiter than in Saturn might stem from the larger energy flux available in Jupiter (Guillot et al. 2004).

6. Prospect for giant planet evolution

While the aim of present study is to examine and to constrain the properties of transport by semiconvection inside Jupiter and Saturn at the present time, the impact of nonadiabatic interiors on the cooling of these planets, and of giant planets in general, remains to be explored. This requires more cumbersome evolutionary calculations, for which initial conditions will be crucial, as will be explored in a forthcoming study. Evolution will add an additional constraint, namely that the planet cooling rates yield the correct properties at the age of the solar system, which will put more stringent constraints on the range of possible layer sizes/numbers. Without going into such detailed calculations, however, the following points can be mentioned.

6.1. Merging of the layers

As mentioned in Sect. 4, soon after they form, layers are expected to merge, leading to thicker layers, until the layer height reaches either the planet size, yielding a standard adiabatic interior, or an equilibrium value, which is the semiconvective case considered here. Numerical simulations by Radko (2005) and Rosenblum et al. (2011) tend to show that the equilibration timescale of the staircase is much shorter than the typical timescale for the evolution of the planet.

Thus, if an equilibrium height is reached, as discussed in Sect. 4.1, this rather quick equilibration timescale suggests that the layer size should remain roughly constant during the evolution, or change slowly with the mean properties of the medium (e.g. the mean molecular weight gradient). Because precise prescription for the height of the layers is still lacking and demands a more precise knowledge of the heat transport properties of layered convection under astrophysical conditions, it seems reasonable, as a first guess, to use a constant size, whose value has been constrained in the earlier sections, throughout the evolution.

6.2. Initial heavy element distribution

Conventional models based on fully adiabatic thermal profiles notably lead to cooling times about 15% longer than the age of the solar system for Jupiter (Fortney et al. 2011). In principle, the hotter non-adiabatic internal structures suggested in the present paper will prolong the cooling and thus worsen the problem.

However, in the case of the erosion of an initially large core, part of the gravitational work will be spent eroding the core and mixing the material upward and will thus not contribute to the total luminosity, thereby quickening the cooling. All these effects must be properly accounted for to determine the appropriate cooling timescale.

In addition, if Jupiter and Saturn initial cores were allowed to be relatively large ( ≳ 10   M_⊕), the corresponding high surface density of solids in the protosolar nebula will quicken the formation timescale in the conventional core accretion scenario, helping solve the related formation timescale problem (Pollack et al. 1996). Finally, since, in the present scenario, some of the ablated material from the accreted planetesimals during the planet’s early formation stages is still distributed throughout the envelope, this will (i) reduce the heating due to gravitational energy release produced by the infalling planetesimals on the planet embryo and (ii) increase the envelope’s mean molecular weight. Both effects will cause the protoplanet to contract more quickly, shortening again the planet’s formation timescale in the conventional core accretion scenario (Pollack et al. 1996). A correct exploration of the impact of inhomogeneous interiors upon giant planet history thus necessitates an investigation of the consequences not only on the thermal evolution but also on the formation process.

7. Conclusion and perspective

In this paper, we have first developed an analytical formalism to describe layered convection, based on a standard parametrization of convection similar to the MLT formalism. This formalism allows quantitative determination of the expected number of diffusive layers or, equivalently, of the average characteristic mixing-length parameter, in a semiconvective planet interior characterized by a given total flux and a given thermal (and compositional) diffusivity. Furthermore, this formalism allows exact determination of the characteristic thermal gradient in the presence of double-diffusive convection, and thus of the related amount of super-adiabaticity within the planet’s interior.

Using this formalism, we computed semiconvective interior models of Jupiter and Saturn. We showed that a stratified internal structure for solar system gaseous giants with a compositional gradient of heavy material extending over a substantial fraction of the planet is a viable hypothesis, because such models can fulfill all the observed gravitational and atmospheric constraints for these planets. This new possibility differs from the conventional description of giant planet interiors, assumed to be composed of two main superposed, well-identified layers of homogeneously distributed material, namely a solid core surrounded by a dominantly gaseous H/He envelope. The consequences of the present giant planet interior description are multiple:

• (i) our Jovian planets might be significantly more enriched inheavy elements than previously thought;

• (ii) their interior temperature, thus heat content, might be much higher than usually assumed;

• (iii) the inner temperature profile could depart significantly from the usually assumed adiabatic profile.

We stress that these conclusions do not depend on the precise model used to describe double diffusive convection. Besides directly affecting our conventional vision of giant planet mechanical, compositional, and thermal structures, these results have profound impact on our understanding of planet formation and cooling properties. Indeed, the revised possible maximum amount of heavy material bears direct consequences on the determination of the efficiency of solid planetesimal accretion during planet formation in the protoplanetary nebula, suggesting an early and efficient capture of planetesimals for ours and, probably also, extrasolar giant planets. Moreover, the greater heat content and the departure from adiabaticity, as well as the possibility of significant core erosion from an initially large core, directly affect the planet cooling histories. Departure from adiabaticity, in particular, implies less efficient heat transport, which is a direct consequence of the inhibited convective motions owing to a persistent compositional gradient, and thus a lower heat flux output rate, than assumed in the conventional approach.

These results open a new window, and raise new challenges, to our present understanding of planet structure, formation, and evolution. Important is that the viability of such stratified interior models for our solar system gas giants applies directly to the case of extrasolar planets, reinforcing the possibility that such a lower heat flux output could at least partly explain the anomalously large radius of several transiting “hot Jupiters” (Chabrier & Baraffe 2007). Indeed, it seems that invoking an extra source of (tidal, kinetic, or magnetic) energy dissipation in these object interiors can not completely solve this “radius anomaly” puzzle and that an alternative or complementary process is necessary (Laughlin et al. 2011). Unconventional, inhomogeneous non-adiabatic planetary interiors, as suggested in the present study, might provide the missing piece of the puzzle.

---

Acknowledgments

We would like to thank our referee D. J. Stevenson for his sharp questions that considerably enlarged the scope of this study. 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)

References

Appendix A: The case of the mixing length theory

Here, we briefly discuss the particular case of the mixing length theory (MLT; Hansen & Kawaler 1994). This case can be recovered in the efficient convection regime by the more general model presented in Sect. 3, by setting C_L = 1 and a = 1/2, but we show below that the MLT formalism can be extended to the inefficient convection regime.

Hansen & Kawaler (1994) showed that in the Pr ≪ 1 regime, recasting their Eq. (5.60) in the notation of Sect. 3, the modified Rayleigh number verifies (A.1)with (A.2)where σ   N_T is the growth rate of a convective eddy and the inverse of the convective time. In general, these equations can be solved numerically. It can be easily verified that in the limit of efficient convection, i.e. Ra _⋆  ≫ 1, σ tends toward unity, meaning that radiative losses are negligible and that the convective time tends toward N_T. In this limit, a simple expansion of Eq. (A.1) yields and as expected in the standard MLT formalism (Hansen & Kawaler 1994).

In the 1 ≫ Ra _⋆  ≫ Pr regime, however, developing Eq. (A.2) yields Then, because the convective time is dramatically increased, Eq. (18) rewrites in that case as (A.3)The size of the diffusive interfaces thus grows until the convective layers eventually turn into a completely diffusive medium for which  ⟨ ∇_T ⟩  = ∇_d.

Appendix B: Computation of the heavy element content

The mean molecular weight gradient needed to stabilize the fluid against large-scale convection can be caused by an inhomogeneous distribution of both helium (Y(m)) and metals (Z(m)) in the hydrogen (X(m))-rich medium. In practice, both gradients can be present at the same time and either compete or contribute constructively.

In our model, we consider an ideal mixture of heavy elements within an H/He envelope whose H/He mass ratio is kept constant and equal to its value in the protosolar nebula, (H/He)_proto. This implies (B.1)everywhere in the planet’s gaseous envelope. We are then left with only one degree of freedom. Following previous calculations (Chabrier et al. 1992), in order to have a flexible determination of the metal enrichment and a thermodynamically consistent EOS in the gaseous phase, we approximate the metal mass fraction by an effective helium mass fraction (Y′) in the H/He EOS. For the core, the metal mass fraction is correctly described by the appropriate water and silicate EOS mentioned in the text.

The various element mass fractions and the corresponding metal enrichment are thus inferred from the relation (B.2)which gives Z at each depth along a given model P-T profile (Chabrier et al. 1992). The hydrogen and helium mass fractions are then derived using (B.3)and (B.4)In this simple model, a Z gradient thus necessarily yields a competing inhomogeneous helium distribution within the planet. Because the mean molecular weight of an H/He mixture at fixed temperature and pressure only depends on Y/(X + Y), only the Z variations need to be considered to compute ∇_μ in our simplified model. The Z(m) profile is then integrated to obtain the total amount of heavy elements mixed in the gaseous layers for each planetary model, as portrayed in Fig. 7. In the most general case, with an intrinsic inhomogeneity of the helium distribution,

caused for instance by its immiscibility in metallic hydrogen, both the Y and Z gradients would have to be properly calculated.

All Tables

All Figures