© S. Markham et al. 2022

1 Introduction

When the molecular weight of a condensing volatile species exceeds the mean molecular weight of the gaseous atmosphere it inhabits, convective inhibition may occur, permitting stably stratified superadiabatic temperature gradients (Guillot 1995). Applications of convective inhibition to our Solar System has garnered renewed interest in recent years; it has been invoked to explain exotic meteorology on Saturn (Li & Ingersoll 2015), the underlying fluid dynamical theory has been further scrutinized (Friedson & Gonzales 2017; Leconte et al. 2017), and it has been applied to refine internal states and thermal evolution of Uranus and Neptune (Markham & Stevenson 2021). These models in general take as an input parameter the bulk mixing ratio of the condensing species at depth, assumed to be well mixed and dominated by hydrogen. The fundamental physics of the case where the mixture at depth is dominated instead by the condensing species, namely, a low-mass envelope overlying a global ocean, has not yet been considered. In this work we will approach the problem using thermodynamic arguments regarding the coexistence of multicomponent systems, generalizing the notion of convective inhibition to apply it to any phase mixture. The new insights in this paper are important because the case we consider, a low-mass envelope overlying a global ocean, is expected to correspond to super-Earth and sub-Neptune exoplanets.

Sub-Neptune and super-Earths were among the very first observed exoplanets, with Poltergheist and Phobetor found orbiting the pulsar Lich in 1992 (Wolszczcan & Frail 1992). The structure and composition of these bodies, intermediate in mass between Earth and Neptune with no analogs in the Solar System, have only gradually come into focus since. In the literature these bodies are currently referred to as either sub-Neptunes or super-Earths depending on the context; in this work we preferentially refer to them as “super-Earths” for clarity, although we are referring more generally to bodies in the mass range roughly between two and ten Earth masses that retain a gaseous envelope. The deluge of data from the Kepler mission represents the greatest advancement in our understanding of these bodies, making clear that planets do not have to become much larger than Earth before they begin to rapidly increase in apparent size (Borucki et al. 2011; Furlan et al. 2017) when defining size using the millibar pressure. These bodies’ mean densities can be described by a composition dominated by water (Seager et al. 2007), although the dearth of observed high-mass, high-density planets suggests many are likely to be a mixture of primarily silicates and hydrogen. A small amount of hydrogen can make an enormous difference in a planet’s apparent radius, because the millibar level can be many scale heights above a physical surface, and the combination of low molecular mass and high temperatures commonly observed in these bodies imply large scale heights. Canonical models usually involve an extended, convective gas envelope in thermal contact with a magma ocean (e.g., Ginzburg et al. 2016; Vazan et al. 2018). This model, which involves a silicate core with an overlying gas envelope, will be the focus of this paper. In particular, we discuss the implications of convective inhibition as a result of the condensation of silicate vapor. This paper is the first in a planned series of papers intended to fundamentally rethink our conception of these bodies, which perhaps behave much more like gas planets than terrestrial analogs with a gas envelope.

Recent observations of Jupiter and Saturn indicate both planets possess an extended diffuse core (Stevenson 2020; Mankovich & Fuller 2021; Dongdong 2019; Miguel et al. 2022), suggesting that significant metal enrichment may go alongside gas accretion during planetary formation (Ormel et al. 2021). In this paper we argue that if the initial envelope of super-Earths is substantially enriched in silicates, as recent formation models contend (Bodenheimer et al. 2018; Brouwers & Ormel 2020), this can meaningfully complicate our understanding of super-Earth structure and evolution. We argue that such enrichment leads naturally to a layer of static stability at depth such that the envelope convectively decouples from the core. In this work we demonstrate that convective inhibition of this kind is likely to occur on silicate planets that possess an envelope on the order on the order of 0.1–10% M_⊕. The phenomenon we describe may also apply to greater envelope masses, depending on currently unknown details regarding the critical mixtures of hydrogen and silicates. Planets with hydrogen envelopes in the mass range 0.1–10% M_⊕ are probably the most common in the universe (Borucki et al. 2011; Furlan et al. 2017).

We argue for a new model for sub-Neptune or super-Earth interiors illustrated in Fig. 1. We expect a hot, supercritical core beneath a cool gas envelope on the following basis. A layer of static stability at depth reduces super-Earth luminosity at early times such that the core loses heat inefficiently. Because of the expected low thermal transport efficiency under the relevant conditions (see Sect. 3.2), we expect the layer of static stability to be thin compared to the full extent of the atmosphere, even though the temperature change across the stable layer can be thousands of kelvins. In the thin stable layer limit, we can model a super-Earth’s steady state internal heat flow based on the mass of its envelope. We argue that the top of the core of these planets lies on the critical line of the hydrogen–silicate binary system at a temperature and pressure that exceed the corresponding values for the critical point of pure silicates, which we call the “critical endpoint” in this work. We expect the critical endpoint pressure level to be on the order of a couple kilobars (Xiao & Stixrude 2018), corresponding to an envelope mass of roughly 0.1% M_⊕. The corresponding base of the convective part of the envelope can be significantly cooler than the top of the core, such that the planet’s internal heat flux can be small even when the core is very hot (~10⁴ K). We predict that the luminosity of planets drops rapidly after isolation, but then remains nearly constant over geologically long timescales. This is an exploratory work focused primarily on describing a new model for super-Earth internal structure and thermal evolution, and constraining the relevant orders of magnitude. We intend to follow up this work with more detailed studies that address questions that naturally arise from this new model, including hydrogen pollution of the core and the process of formation.

This work follows the following structure: in Sect. 2, we extend the theory of convective inhibition to generalized thermodynamic contexts, using a water–hydrogen system (for which, unlike silicates, we have data) as an example to illustrate the relevant thermodynamics at play. In Sect. 3, we outline the internal structure of the super-Earths we model in this work. In Sect. 4, we calculate orders of magnitude for core cooling timescales and argue that under certain conditions the core can retain much of its primordial heat for billions of years. Then in Sect. 5 we identify several natural consequences of the super-Earth model we present in this work, but that lie beyond the scope of this work, to motivate collection of further data and follow-up papers in this series. Finally in Sect. 6 we summarize our findings, comment on observable consequences of our model, and identify the key data necessary to further advance our understanding of the most common class of planets in our galaxy.

Fig. 1 Cartoon representation of this paper’s findings. On the left is a representation of the canonical view of super-Earths (e.g., Ginzburg et al. 2016), with a frozen or molten silicate core in good thermal contact with a convective gas envelope. On the right we show our proposed model, with a layer of static stability between a (possibly hydrogen-polluted) supercritical fluid core and the silicate-vapor-enriched gas envelope.

2 Theoretical Foundations

In this section we aim to introduce the general notion of convective inhibition for unfamiliar readers, and to generalize previous theoretical work to apply in more exotic contexts. This paper follows the preestablished logic of convective inhibition by condensation in hydrogen atmospheres (Guillot 1995). Prior works focus on planetary atmospheres in the Solar System and assume that the liquid phase remains pure, and that the saturation vapor pressure follows Clausius-Clapeyron with constant latent heat. These approximations are appropriate in the context of low pressures and where the mixing ratio of the vapor is much less than unity. While these approximations are valid in the context of Solar System atmospheres, they break down at higher pressures and temperatures likely to be relevant on super-Earths. Therefore we wish to extend the logic to a general mixture of coexisting phases that remains applicable much deeper into the interior of planets. Of particular interest to this work is the case of a relatively low-mass hydrogen atmosphere with an infinite reservoir of volatile condensates (i.e., an ocean). This system will be of particular interested because observations indicate a large population of planets with masses intermediate between that of Earth and Neptune with radii far larger than can be explained with an iron-silicate composition (e.g., Furlan et al. 2017). If we adopt the prevailing interpretation, that the observed radii are due to hydrogen envelopes comprising about 1% of the planet’s mass, then we expect envelope basal pressures to exceed the critical endpoint pressure of silicates (1–2 kbar; Xiao & Stixrude 2018, M_atm ~ 10⁻³M_⊕). Under these conditions, it is no longer appropriate to approximate the coexistence of the condensate with the dry gas according to the saturation vapor pressure of the pure condensate (which we refer to as the ideal gas approximation in this work); one must use coexistence data. Because this scenario violates the assumptions underpinning prior derivations, a new more generalized framework is needed.

In this section we inspect three crucial differences in the regime of interest compared to prior works. In Sect. 2.1 we re-derive the phenomenon of convective inhibition in a framework that is robust to generalized thermodynamic environments. In Sect. 2.2 we seek to specify the top and bottom boundaries of a stable layer resulting from a critical mixture of a condensing species with gas. Then in Sect. 2.3 we motivate the necessity for Sect. 2.1 by indicating the importance of empirical coexistence curves for mixtures at supercritical pressures. In Sects. 2.2 and 2.3 we use data for a water/hydrogen system to illustrate the behavior with data, and note that the qualitative behavior for silicates should be similar but scaled to an order of magnitude higher temperatures and pressures.

2.1 Convective Inhibition in a Generalized Mixture

Convective inhibition for an ideal gas coexisting with a condensing species has been derived several times in the past (Guillot et al. 1995; Leconte et al. 2017; Friedson & Gonzales 2017). The criterion is often written in terms of mass ratio (1)

where μ_A and μ_B are the molar masses of the dry gas and condensing vapor respectively, and L is the latent heat of vaporization. This can be equivalently written in terms of the inhibition molar concentration (Markham & Stevenson 2021). (2)

where ϵ ≡ μ_B/μ_A. In this section we re-derive this criterion for a generalized phase mixture, making no explicit assumptions about the equation of state of the mixture. In the context of this work, we do this to apply the concept of convective inhibition to a critical (non-ideal) mixture; more generally, it can apply to any coexisting species of different densities. We put aside the question of the background temperature gradient for now, instead focusing purely on how a temperature increment affects density of an isobaric parcel in thermodynamic equilibrium. As we will demonstrate, this approach is sufficient to delineate convective inhibition, and recovers the same result as prior works in the ideal gas limit.

The coefficient of thermal expansion is (3)

where V is the molar specific volume of the material, and x refers to the composition. We postulate that convection operates as normal (heat rises) when the coefficient of thermal expansion α is positive. On the other hand, materials with a negative coefficient of thermal expansion will be stable against convection even for a superadiabatic temperature gradient. That is the reason why, for example, water in a pond will freeze from the top – the coefficient of thermal expansion for water is negative between its freezing point and about 4°C, where it is near zero. To assess when convection can be shut off in planetary interiors, we seek the point at which the coefficient of thermal expansion becomes zero. However, for most pure substances in the universe, the coefficient of thermal expansion is always positive; liquid water is a special case. For a pure ideal gas, for example, the coefficient of thermal expansion α = 1/T > 0 ∀ T. As has been demonstrated before, however, convection can be shut off for an ideal gas if its composition can be affected by temperature. Therefore we define the “effective” coefficient of thermal expansion as (4)

which accounts for the possibility that the composition and therefore the mean molecular weight μ of a fluid parcel may also depend on temperature. This applies to, for example, a gas parcel in contact with a volatile liquid ocean in thermodynamic equilibrium. In that case, increasing its temperature allows additional vapor molecules to enter the parcel, while decreasing its temperature will cause vapor molecules to condense and eventually rain out. It can also apply more generally to a variety of contexts relevant to planetary interiors, for example a eutectic system with a dense phase that freezes out, two partially miscible liquids of different densities, or an equilibrium chemical reaction. This work will specifically focus on the case of a critical mixture (see Sect. 2.3), although the derivation in this section is more general.

We now consider an isobaric mixture of two species A and B coexisting in a parcel of fluid, where the molecular mass of species B is greater than that of species A, μ_B > μ_A. The parcel contains a fixed mass of species A that does not undergo any phase transitions, in contact with an infinite reservoir of species B. Species B can enter or leave the system according to its equilibrium concentration.

In our system the molar concentration x ≡ N_B/(N_A + N_B) uniquely determines the composition of the parcel, and its mean molecular weight μ = μ_A (1 − x) + μ_Bx. Since N_A is fixed, x = 0 when N_B = 0 (pure substance A), and x → 1 as N_B → ∞ (when B becomes completely soluble). It follows directly from Eq. (4) that the effective coefficient of thermal expansion for this system is (5)

It should be noted that the first term on the right-hand side of Eq. (5) is equivalent to the classical coefficient of thermal expansion (Eq. (3)), while the remaining terms are corrections to account for changes in composition. The second term on the right-hand side accounts for the possible change in mean molecular spacing that occurs as a result of changing the composition, possibly due to new intermolecular forces between species A and B. The third and final term on the right-hand side accounts for the change in mean molecular weight associated with changing composition; this term can be negative if ϵ > 1 ⇐ μ_B > μ_A.

Convection is inhibited if α_eff → 0. Then rewriting Eq. (5), we obtain the generalized criterion for convective inhibition. (6)

If one specifies a differentiable miscibility x(T), a deviation from volume additivity , and a classical coefficient of thermal expansion α, one can solve Eq. (6) for T, the temperature at which convection will be shut off in equilibrium. Substituting this temperature into x(T) implies x_inh, the inhibition concentration¹, if a solution exists. If no solution within the domain 0 < x < 1 exists, then convective inhibition does not occur. As an example, we can solve for x_inh, for an ideal gas A in contact with a partially vaporized volatile liquid B we can specify , and (Clausius-Clapeyron), where L is the latent heat of vaporization for species B and R is the ideal gas constant. In that case we recover Eq. (2), identical to previous works.

The exact value of the inhibition concentration apparently depends on temperature, but we note for the convenience of the reader that it is on the order of 0.1–1% in the thermodynamic regimes of interest. For all saturated concentrations above the inhibition mixing ratio, superadiabatic temperature gradients will be Ledoux stable against convection.

We have so far not fully justified why it is sufficient to consider only the isobaric coefficient of thermal expansion. To do this we now perform simple linear stability analysis of a parcel in an environment that is in saturated thermodynamic equilibrium with an imposed temperature gradient . We want to assess whether a parcel displaced upward will be more or less dense than its surroundings. For the environment, (7)

For the parcel, (8)

where ∇_ab is the adiabatic temperature gradient that can be specified and may take different forms depending on the Grüneisen parameter of the materials, latent heat released, etc. A convective instability occurs if after some small upward displacement −δp, ρ_parcel < ρ_env. Therefore by subtracting Eqs. (8) from (7), the condition for a convective instability is (9)

(cf. Eq. (2) of Guillot et al. 1995 or Eq. (15) of Leconte et al. 2017). Therefore, to tell whether a superadiabatic temperature gradient ∇_T > ∇_ab will be stable against convection in an inviscid medium, it is sufficient to only calculate the sign of the isobaric effective coefficient of thermal expansion α_eff (Eq. (5)); the results are independent of the precise form of ∇_T and ∇_ab and also independent of the material’s compressibility.

One may question whether the kinetics and dynamics of these processes are relevant, and whether one can truly make the assumption that the system is in thermodynamic equilibrium. These are reasonable objections that have been investigated in detail in prior works (Friedson & Gonzales 2017; Leconte et al. 2017). These works find convective inhibition to be robust with finite settling timescales and imperfect molecular transport, although there are limits. Detailed microphysical modeling is beyond the scope of this work, and we encourage interested readers to refer to the above cited works for a more detailed analysis of these effects.

2.2 Boundaries of the Stable Layer

We seek the top and bottom boundaries of the stable layer. Evidently the inhibition mixing ratio at the top boundary depends only on the local gradient of the condensate–gas coexistence curve. Under the ideal gas approximation at constant pressure, then, x_inh becomes a function of temperature only. At higher pressures, an empirical coexistence curve must be used. Now that we have identified the mixing ratio at which convection is shut off, we wish to determine at which mixing ratio the convection reasserts itself. For example, in Markham & Stevenson (2021) this occurs after the mixing ratio reaches its “bulk” value, called q_max. However, this bulk value as such only applies in the limit where the system is dominated by gases. In the case we consider now with an infinite reservoir of condensing species, we must determine the mixing ratio at depth by other means. The solution depends on the thermodynamic conditions.

In this section, we refer to “subcritical” and “supercritical” conditions by referring to the critical endpoint of the pure condensate. We use the notation to refer to the critical endpoint of the pure condensing substance, to distinguish it from the critical curve (p_crit, T_crit)(x) which depends on the mixing ratio of condensate to hydrogen. The endpoint of the critical curve where x = 1 corresponds to the aforementioned critical endpoint.

We first consider the case where the hydrostatic pressure is subcritical. In this case, there exists a temperature for which the saturation vapor pressure approaches the hydrostatic pressure. Physically, this situation manifests as a stable layer truncated from below by a liquid ocean.

If the hydrostatic pressure is supercritical, then the behavior is subtler. In this case, the condensate abundance will increase with increasing temperature, stabilizing against convection, until exceeding the critical temperature, at which point the condensate and gas can mix in all proportions. When the two substances can mix in all proportions, condensation no longer occurs so that convective inhibition is no longer at play, and free convection can resume. Figure 2 shows the limits of the ideal gas approximation when dealing with supercritical conditions.

We see from Fig. 2 the behavior of a coexisting system in pressure environments exceeding the critical endpoint pressure . For marginally supercritical hydrostatic pressures, the ideal gas approximation is appropriate for determining the top of the stable layer, and the bottom of the stable layer is near the critical endpoint. However, we see for sufficiently high pressures (red curve in Fig. 2) that deviations away from the ideal case become important, suggesting somewhat lower temperatures at the top of the stable layer and higher temperatures at the top of the core. Hydrogen enrichment at the top of the core may then be significant. Indeed, hydrogen enrichment at the top of the core may be quite substantial, approaching a 1:1 molecular mixture at high pressures. Following figures will include the same symbols for the important intersection points.

Fig. 2 Illustration showing the stable layer boundaries in (x, T) space. Solid curves use data (Seward & Franck 1981, see Sect. 2.3), while dashed curves use the ideal gas approximation x = ps/p, diverging near the critical endpoint. Gray shaded areas are regions of parameter space wherein free convection is permitted. In the top-right corner, the gray region lies beyond the critical curve where phase separation does not occur. At the bottom of the figure, the gray region lies below the xinh curve described in Eq. (2). Blue, magenta, and red curves represent hydrostatic pressure environments of 300 bars, 1 and 2.5 kbars respectively. The solid points where the dashed curves intersect correspond to the top of the stable layer under the ideal gas approximation; the red asterisk shows the offset intersection using data. The open circles where the solid colored curves intersect the critical curve represent the top of the core.

2.3 Critical Coexistence of Rock and Gas

The approximation that x_s(T) = p_s(T)/p where p_s is the saturation vapor pressure and follows an Arrhenius relationship in temperature with constant latent heat is a fair approximation when . However, in cases where the saturation vapor pressure exceeds the hydrostatic pressure or approaches the critical pressure, as is of interest in this work, such an approximation is no longer valid, and empirical data are necessary to model the coexistence of phases as shown in Fig. 2. The solid curves in Fig. 2 show empirical coexistence curves of a gas-liquid mixture of water and hydrogen (Seward & Franck 1981; Bailey & Stevenson 2021). At low pressures, the liquid phase remains nearly pure, while more and more water enters the gas phase as temperature increases. Above the critical temperature of water ~650 K, the phases no longer separate and both species can coexist in a well-mixed supercritical fluid. Figure 2 shows the deviation of this empirical data compared to the assumption that , where L is the latent heat chosen to be 1.85 × 10¹⁰ erg g⁻¹ for that illustration, corresponding to the average value between the melting temperature and the critical endpoint temperature.

This critical curve in Figs. 3 and 4 interpolates between data from Seward & Franck (1981) and Bali et al. (2013). The shaded inhibition zone in Fig. 3 corresponds to the equivalent shaded areas in Figs. 4 and 5. We shade this particular region for convenience, because the behavior of the inhibition and critical curves are unknown at higher temperatures and pressures, although the reader should understand that in reality the inhibition zone should extend everywhere between the inhibition curve and the critical curve. This is a key point that warrants emphasis; in Sect. 4 we use the inhibition zone as shown in red in Fig. 4. That shaded zone bears a sharp point at the bottom, at which the inhibition temperature approaches the critical end-point temperature. However in reality, the critical temperature depends on composition, we actually expect the critical temperature to curve strongly toward higher temperatures as shown in Fig. 3.

We expect the basic thermodynamics to play out in much the same way for a mixture of silicates and hydrogen. We show roughly where we expect the silicate inhibition zone to be in Fig. 4, and will use this inhibition zone again in Fig. 5. For these figures we use the ideal gas approximation for silicates, with the saturation vapor pressure of pure silicates taken from ab initio simulations (Seward & Franck 1981). We use μ_B = 20 for silicates, corresponding to the mean value for an atomic disaggregation of SiO₂, while we use μ_A ~ 1.27 corresponding to a mixture of atomic hydrogen and helium with a 10% helium number fraction. Different choices for these values of μ_B and μ_A affect the placement of the inhibition curve by changing ϵ ≡ μ_B/μ_A; from Eq. (2) one can see that increasing e tends to decrease x_inh and vice versa.

Although we use the ideal gas approximation to delineate the rock inhibition zone in Figs. 4 and 5, in reality we expect the coexistence curves for the silicate–hydrogen system to look qualitatively similar to the water–hydrogen system, albeit with both pressure and temperatures increased by about an order of magnitude owing to the order-of-magnitude larger bond energies within silicate atoms compared to weaker hydrogen bonding between water molecules. Therefore, we expect the mixing ratio beneath the stable layer to depend on the pressure formation conditions. The mixing ratio immediately beneath the stable layer is set by the intersection of the solid black and colored curves in Fig. 2, marked with open circles. One should take note that the critical curve veers strongly to the left in Fig. 2, so that the behavior of the critical curve is very important when . Figure 2 is a major finding of this work, as it demonstrates that even with an unlimited abundance of condensates, a stable layer due to convective inhibition will still have well defined boundaries. Furthermore, in general our model permits hydrogen to exist below the stable layer. We discuss hydrogen pollution of the core in the conceptual sense in more detail in Sect. 5.1.

We must also consider chemistry particular to silicates that does not apply to ices that condense at temperatures cool enough that individual molecules do not dissociate. For silicates, unlike water, the atomic composition of the condensate in the condensed phase will not in general be identical to the atomic composition of the vapor phase (Xiao & Stixrude 2018). Furthermore, unlike a mixture of hydrogen and water, hydrogen can react chemically with silicates, complicating the approximation of pure substances in coexisting phases (e.g., Fegley & Lodders 1994). The chemistry is varied and complex, but likely one example is the equilibrium between silicate vapor and silane gas, for example, (10)

Detailed modeling of the complicated series of hundreds of simultaneous equilibrium chemical reactions is beyond the scope of this work (but has been modeled in prior works; see e.g., Schaefer & Fegley 2009), though we note its effect with respect to convective inhibition can still be accurately modeled using the generalized approach from Sect. 2.1².

Having now discussed in detail the complications arising from thermodynamic considerations of critical mixtures primarily using existing data for water to illustrate, we now will transition to discussing practical consequences for super-Earth internal structure and evolution in the following section. As shown in Fig. 4, Kepler 36c is an example of a planet whose early temperature structure intersects the rock inhibition zone, indicating that subsequent thermal evolution will be affected by convective inhibition. To attempt to quantify the effect convective inhibition will have on subsequent thermal evolution, we use the ideal gas approximations for convenience, with the inhibition zones demarcated by the red shaded region of Fig. 4. The boundaries of this region, the ideal gas inhibition curve, the ideal gas boiling curve, and the critical endpoint temperature, are chosen for convenience because they are values for which we have preliminary data for silicates (Xiao & Stixrude 2018). However, we hope in this section we have impressed upon the reader that these approximations are incomplete when considering a critical mixture, and our results should be regarded as notional orders of magnitude rather than quantitatively exact. Nevertheless, even without more precise data, there are already an abundance of impactful consequences of convective inhibition with an ocean.

Fig. 3 Temperature at the boundaries of the stable layer as a function of envelope basal pressure, using water. The blue and red curves correspond to the top and bottom of the stable layer respectively. As in Fig. 2, we use solid curves where we have data and dashed curves for the ideal gas approximation. We further use a dotted curve for the interpolated critical curve and for a notional (non-quantitative) altered inhibition curve that accounts for the non-ideal coexistence behavior of the water–hydrogen system. The intersection marker symbols from Fig. 2 (the blue, magenta, and red filed and open circles and the red star) are shown here for context.

Fig. 4 Examples of pressure-temperature profiles for planets in our Solar System and galaxy. If the planet mass is dominated by an ocean (infinite reservoir) with a low-mass hydrogen envelope, convection will be inhibited in the shaded region, due to water in the blue region and silicates in the red region. While the red region uses the ideal gas approximation, the blue region extends into higher pressure regimes; we use a color gradient and question mark to emphasize the behavioral uncertainty under these conditions. Data for Uranus are extrapolated from Lindal et al. (1987), the lupiter model is from Hubbard & Militzer (2016), and the Kepler 36c model shortly after isolation from Bodenheimer et al. (2018). Colored points correspond to the same symbols from Figs. 2 and 3, with added yellow square points corresponding to data from Bali et al. (2013).

Fig. 5 Envelope pressure–temperature profiles for a sample of envelope properties. The blue and red curves represent the temperature at the top and bottom of the stable layer, respectively; for a saturated gas mixture, convection is shut off between the blue and red curves. In the absence of data, this figure assumes xs = ps/p and does not consider the temperature dependence of latent heat or the compositional dependence of the critical curve. See Fig. 7 for a conceptual sketch of the temperature structure below the envelope.

3 Structure

Having established the theoretical foundations for convective inhibition by silicate vapor, we now seek to delineate under what circumstances these considerations are relevant, and to comment on the potential effects convective inhibition will have on super-Earth internal structure. In Sect. 3.1 we outline a simple envelope model. We find for a planet in a state of convective inhibition, that the mass of the envelope above the stable layer together with its equilibrium temperature uniquely defines its structure and luminosity. In Sect. 3.2 we place limits on the thickness of the stable layer. Finally, in Sect. 3.3 we speculate on the planet’s properties below the stable layer, and argue for why a more precise understanding must be predicated on so far unavailable data.

3.1 Envelope Structure

We consider the envelope structure under the following approximations: neglecting self-gravity of the atmosphere, an ideal gas equation of state, an adiabatic temperature gradient below the radiative-convective boundary (RCB), and a plane-parallel geometry. We further use the approximations for the coexistence of silicate vapor with hydrogen gas outlined in the previous section. Under the cold Murnahan equation of state, the radius of a silicate core scales as the quarter power of its mass (Seager et al. 2007) (11)

We model the saturation vapor pressure as (12)

where p_s is the saturation vapor pressure in GPa, and the parameters A and B are 11.8 and 45000 respectively, fitting ab initio simulations from Xiao & Stixrude (2018). Convective inhibition is initiated at pressure p₁ and temperature T₁ satisfying (13)

Under our approximations, the total pressure at the bottom of the envelope (or equivalently the top of the core) is where M_e is the total mass in the envelope. Conveniently, because (Eq. (11)), this can be rewritten to be independent of core mass and radius: (14)

The total pressure at the core includes both gas mass M_XY and some mass from silicate vapor, which can be non-negligible. The total mass of the envelope, then, is (15)

where ϵ is the molecular mass ratio of volatile vapor to dry air. In the following calculations we assume the stratosphere is cold enough that it contains negligible silicate vapor such that the lower bound of the right-hand side integrand of Eq. (15) becomes p_rcb instead of p = 0. We also assume p₁ ~ p_c (see Sect. 3.2). The above equation is independent of core mass. Therefore for a given envelope mass, the pressure at the top of the core is uniquely determined.

Because T₁ is uniquely determined for a given p₁ according to Eq. (13), this sets the adiabat for the convective part of the envelope. Highly irradiated super-Earths are expected to have effective temperatures very near their equilibrium temperatures (e.g., Ginzburg et al. 2016), and the temperature at the RCB will be near the effective temperature. Using the adiabatic relationship (16)

where is the adiabatic temperature gradient which is invariant if the Grüneisen parameter is invariant. Following Owen & Wu (2017), we take γ to be a constant 5/3. By combining these considerations, we can uniquely determine the RCB pressure by specifying the total mass M_XY of gas in the envelope. The intrinsic luminosity scales with the depth of the RCB (Ginzburg et al. 2016) (17)

where is the so-called “modified Bondi radius,” (18)

Equation (17) intuitively states that the planetary flux is inversely proportional to the optical depth of the RCB, i.e., , and is valid for large τ. Following Freedman et al. (2008), we assume the opacity scales linearly with pressure, so that the planetary flux diminishes rapidly with increasing depth of the radiative-convective boundary. We plot the planetary flux as a function of envelope mass in Fig. 6. Figure 6 compares two assumptions for atmospheric models. In both cases, the basal pressure is determined by the envelope mass. In the inhibited case, the corresponding basal temperature is T₁ from Eq. (13). In the uninhibited case the basal temperature is . These values are chosen so that the heat of the core is the same in both models.

The reader may notice that the 10⁻¹ M_⊕ curve in Fig. 5 intersects the inhibition curve very near the sharp point at the bottom of the labeled inhibition zone. That is the reason in the abstract we say that our results pertain to planets with envelope masses up to 10% M_⊕, because at higher pressures the adiabat does not intersect the labeled inhibition zone. However, we recall from Sect. 2 that the sharp point is in fact an artifact of our approximations, and that in reality we expect the critical curve to depart from the critical endpoint temperature, so that convective inhibition may still relevant to larger envelope masses. This approximation is also related to the fact that the uninhibited model flux approaches the inhibited model flux around 10⁻¹ M_⊕ in Fig. 6. Again, this is an artifact of our approximations, and the behavior would be different with a more realistic model for the inhibition and critical curves (see Fig. 3).

Figure 6 is centered around an envelope mass of 10⁻²M_⊕. We note for clarity that this value is in Earth masses, not as a fraction of the core mass. So, for example, a 5 M_⊕ super-Earth that is 1% gas by mass should be read as 5 × 10⁻² on Figs. 6 and 8. We further note that an atmospheric mass of ~10⁻²M_c is a realistic estimate for super-Earths based on hydrostatic equilibrium between the core and the gas disk during formation (see, e.g., Ginzburg et al. 2016).

The behavior of Fig. 6 shows a monotonically decreasing relationship for planetary flux as a function of envelope mass. Intuitively, one can understand the behavior this way: the pressure at the bottom of the envelope scales linearly with envelope mass. Under the thin stable layer approximation, we can then use the relationship from Eq. (6) to uniquely determine the temperature T₁ at the bottom of envelope such that x → x_crit (see Fig. 5). We note that in Fig. 5, the thermodynamic conditions at the base of the envelope, (p₁, T₁), depend on envelope mass only. Higher equilibrium temperature planets have deeper RCBs. Given the effective temperature, we can use the adiabatic relationship to place the RCB. Assuming fixed T₁, we can then use Eq. (17) to estimate that the luminosity (or flux, for fixed radius) should scale approximately as the inverse square of the envelope mass. In reality this is an overestimate; accounting for the fact that T₁ will increase as p₁ increases accounts for the somewhat shallower dependence of flux on envelope mass.

Figure 6 also shows that the flux of the planet is not too sensitive to the planet’s effective temperature, and that the flux decreases with increasing effective temperature. From Eq. (17), we see that , where the flux is the core-envelope internal heat flux (this proportionality holds when ∇_ab ~ 1/4). For fixed M_atm ⇒ (p₁, T₁) will likewise be fixed. Therefore we can infer the expected scaling for prcb from Eq. (16). Therefore, by combining Eqs. (17) and (16) using ∇_ab = 1/4, we then find the scaling relationship for flux at fixed envelope mass to be . This result is consistent with Fig. 6 that shows internal heat flux decreasing with increasing equilibrium temperature. We note that this result is sensitive to the adiabatic lapse rate, set by the Grüneisen parameter.

The critical pressure of silicates is about 1.4 kbar (Xiao & Stixrude 2018), corresponding to an envelope mass of about 10⁻³M_⊕. For much lower mass, the flux is sufficiently high that the core can cool quickly, and prior models subsequently work to a good approximation. Therefore for all cases of interest to us, the core will initially be in a supercritical state, rather than a magma ocean.

Fig. 6 Core-envelope heat flux as a function of envelope mass. Temperatures in the legend correspond to Teq. Symbols on the left of the plot show the contemporary internal heat fluxes of Earth, Neptune, and Jupiter for comparison.

3.2 Estimating the Thickness of the Stable Layer

We must challenge our assumption that the stable layer can be approximated as infinitesimally thin. Furthermore, we will challenge our assumption that the layers of the planet behave distinctly – either fully convective or fully stable – given the intrusion of eddy diffusivity. The thermal transport properties of the materials of interest under the relevant thermodynamic conditions are poorly constrained, but we can set reasonable estimates on their order of magnitude. We consider three thermal transport mechanisms of interest: conductive, radiative, and advective.

In the absence of thermal transport data at high temperatures and pressures, we resort to the ideal gas estimate . Thermal conductivity should be roughly independent of density, and the mean free path λ is inversely proportional to density. Using characteristic numbers, and using the critical temperature of silicates as an upper bound, we find the thermal conductivity to be k ~ 10⁵ g cm s⁻³K⁻¹.

Radiative heat transport in the limit where the mean free path is small compared the length scales of the problem (valid for high pressures), we can estimate the equivalent thermal conductivity due to radiative transfer as . The relevant parameters are reasonably well-constrained in this context, except the opacity. The thermodynamic conditions of interest are thousands of kelvins, set by silicate vapor pressure curves, and kilobars to hundreds of kilobars, set by the mass of the envelope. The high pressures imply a large amount of collisional opacity of hydrogen, while the temperatures are so high that partial ionization and highly opaque free electrons are expected. Furthermore in our case the stable layer is a mixture of hydrogen and silicates, so that the stable layer should be a cosmic mixture of many elements, not dominated by hydrogen or any other single substance that may have transparent windows in its transmission spectrum. We therefore expect this environment to be highly opaque. Following Freedman et al. (2008), the opacity for 1 kbar to hundreds of kbar (the pressure range of interest) should be significant, plausibly on the order of 1–100 cm² g⁻¹. To be cautious we assume a low opacity on the order of 1 cm² g⁻¹. In this case we obtain a characteristic thermal conductivity on the order of 10⁷g cm s⁻³ K⁻¹. If the material is more opaque, the corresponding equivalent thermal conductivity will be smaller, plausibly comparable to the microscopic thermal conductivity ~10⁵g cm s⁻³ K⁻¹. Our upper estimate is two orders of magnitude larger than the estimate for conduction, and would therefore dominate thermal transport.

Finally, we consider advective heat transport. Although the stable layer is called as such because it is stable against ordinary large-scale overturning convection, we cannot assume this layer is completely stagnant. We know statically stable layers of the Earth’s atmosphere can nevertheless involve considerable eddy diffusivity due to, for example, breaking gravity waves (Dornbrack 1998; Garcia & Solomon 1985). We can attempt to constrain the order of magnitude of this process, although it is highly uncertain. In order to move material upward, one must do work against gravity. We can assume as a first approximation that the eddy diffusivity behaves like thermal diffusivity, but acts on deviations away from adiabaticity. For heat flux on the order of 10³ erg cm² s⁻¹ (see Fig. 6), this predicts an eddy diffusivity of comparable order-of-magnitude to ordinary thermal diffusivity. Accounting for all these considerations, it appears radiative heat transport is most likely to dominate.

Using our rough upper bound for equivalent thermal conductivity of the system, we can estimate the thickness of the stable layer. From Fourier’s law, F = −k∇T, or L ~ k∆T/F. For a 1% M_⊕ envelope, we estimate that F should be on the order of 10³ erg cm⁻² s⁻¹, and use a characteristic value for the temperature difference ∆T ~ 10³ K (see Fig. 5). Therefore the thickness of the stable layer should be on the order of 100 km, or about <1% the expected radius of a super-Earth. The corresponding pressure drop would be on the order of a kbar, about 10% the pressure overlying the stable layer for a 1% M_⊕ envelope. For smaller atmospheric masses, the thermal flux will be greater and the density/opacity smaller, while for larger atmospheric masses the thermal flux can be smaller but the density/opacity larger. Therefore for our conservative estimates for relatively efficient radiative heat transport, the stable layer can be non-negligible in thickness and mass content, and this will add a small but non-negligible correction to Figs. 6 and 8.

3.3 Below the Envelope

We will keep our discussion of the planetary structure beneath the envelope largely conceptual and brief, as this preliminary exploratory work primarily focuses on the part of the planet above the stable layer (p₁ is Fig. 7). The envelope structure we propose here is very similar to prior models (e.g., Ginzburg et al. 2016). However, we will comment here on the structural assumptions that underpin Sect. 4, as well as comment on the shortcomings of these assumptions to be further addressed in Sect. 5.

The bottom of the stable layer truncates at the critical curve (see Fig. 2). For the purposes of Sect. 4, we refer to everything below the stable layer as the “core”, although this terminology may be misleading. Thermodynamics permits the layer beneath the stable layer to be highly enriched in hydrogen. The exact enrichment depends on formation conditions and the details of the coexistence curve under the relevant conditions, and there may be additional compositional gradients and boundaries (cf. Bodenheimer et al. 2018; Brouwers & Ormel 2020). Furthermore, the behavior of the core at depth depends on currently unavailable data about the physical properties of the supercritical rock/gas mixture under the relevant conditions. For example, the relative steepness of the critical curve and the adiabat can have profound implications for the expected interior structure (cf. Bailey & Stevenson 2021). We therefore refrain from advancing a detailed core model, which would be premature without the requisite data.

Nevertheless, we must employ some model. Having established the impossibility of an accurate model, we therefore opt for a simple model. Throughout this work, we have treated the core radius as scaling like the cold Murnahan equation of state. As the name implies, this equation of state requires corrections for a hot, supercritical core. However, as we argued in Sect. 2, there are currently too many unknowns to attempt a serious dynamic core model, as the results would depend on too many unsubstantiated assumptions about the behavior of critical mixtures of rocks and gas. Therefore, in Sect. 4 we model the core as a well mixed, fully convective incompressible magma. We apply this treatment as an end-member case for the purpose of placing a constraint on the orders-of-magnitude of the timescale. However, in reality the core is likely to be supercritical and significant quantities of gas may pollute the core. Section 5.1 qualitatively analyzes the possible effect of such core pollution on our results.

In our model we consider cases where the heat flow is fully balanced, that is, radiative-convective quasi-equilibrium. Therefore, the heat flow across the stable layer is balanced by heat flow into space. Because we regard the stable layer to be thin (p₁ = p_c), then, the atmospheric structure remains static until the core cools to below T₁, at which point convective inhibition ceases and convective thermal contact between the envelope and the core asserts itself. Subsequent cooling then progresses as prior theories (e.g., Ginzburg et al. 2016) suggest. Now we must compute the timescale on which the core cools from an initially hot (T ~ T_crit) state down to the temperature at the top of the stable layer.

Fig. 7 Sketch of the proposed temperature structures of silicate planets with gaseous envelopes described in Fig. 1. See Fig. 5 for explicitly quantified envelope pressure-temperature diagrams; this figure is only intended to be an illustrative sketch.

Fig. 8 Cooling timescales for a super-Earth computed using Eq. (20). The dotted curves show the cooling timescale without convective inhibition with the core temperature at . The solid curves include convective inhibition as described in this work.

4 Evolution

Previous works (e.g., Vazan et al. 2018) have identified that the existence of a gas envelope can considerably slow the thermal evolution process, such that a magma ocean evolutionary phase may persist for billions of years. Our findings here indicate that thermal evolution of the core could take place even more slowly, owing to the inefficient thermal transport through a stable layer. We begin by computing the cooling timescale of a body beginning in the state described in the previous section, with the top of core at the critical temperature and the bottom of the envelope at T₁ such that convection is inhibited. We assume in the radiative-convective quasi-equilibrium state that the envelope temperature profile is frozen, that is does not change in time. This occurs because, if the luminosity of the planet exceeds the heat transport from the core the envelope, then the convective envelope will cool, further thinning the stable layer, until the stable layer is sufficiently thin that the heat flow through the stable layer balances the heat flow out of the planet. In this section we argue that even if the planetary luminosity initially exceeds the heat flow between the core and envelope, the envelope will quickly relax into the equilibrium state. In radiative-convective quasi-equilibrium, (19)

The cooling timescale that describes how long the planet remains in this stabilized state, then, is (20)

We assume the heat capacity of the core agrees with the ideal gas limit (Bolmatov et al. 2013), consistent with prior works (e.g., Ginzburg et al. 2016). Assuming the luminosities computed in Fig. 6, we compute cooling timescales shown in Fig. 8 assuming a 5 M_⊕ core. The cooling timescale for the solid curves in Fig. 8 use Eq. (20). However, when considering a planet cooling without convective inhibition (dotted curves), we must account for the non-constant luminosity, as the luminosity of the planet will diminish as the atmosphere and core cool together. We can compute the core cooling time analytically by integrating Eq. (20) while accounting for a non-constant luminosity. We can rewrite the luminosity from Eq. (17) as (21)

As the atmosphere cools, the RCB becomes deeper. Under our assumptions, the relationship between p_rcb and the temperature at the bottom of the atmosphere T_btm is (22)

Therefore we have a relationship between the time-dependent top of core temperature, and the luminosity. We know the atmospheric state when and at the end state when T_btm = T₁. Therefore we can integration Eq. (20) by substitution to obtain (23)

We use this analytic expression for the dotted curve in Fig. 8.

The upper abscissa in Fig. 8 occurs because in this simplified model when M_atm ~ 10% M_⊕ (the bottom of the atmosphere in Fig. 5 is very close to the tip of the inhibition zone); therefore, the numerator in Eq. (20) becomes smaller. We use this simplifying assumption due to the absence of data, but in reality the tip we see in the silicate inhibition zone in Figs. 4 and 5 probably does not behave like that. In reality, we expect the actual critical curve for a silicate-hydrogen mixture to increase at high temperatures (see Sect. 2, Fig. 3).

The cooling times in Fig. 8 refer to the time it takes the core to cool sufficiently that convective inhibition is no longer important. These timescales are somewhat longer than the cooling timescales for the uninhibited case, by about a factor of 5 for a 1% M_⊕ hydrogen envelope. The calculations for Fig. 8 were performed assuming a 5 M_⊕ core, consistent with Fig. 6, but the cooling timescale results are actually independent of core mass under our assumptions. From Eq. (21), the luminosity of core scales linearly with the mass of the core. Because the heat content of the core also scales linearly with its mass, the cooling timescale is unaffected by this choice. One can also see that our analytic expression for the cooling time in the uninhibited case, Eq. (23), is independent of M_c for the same reason. This result is a consequence of our assumption for the equation of state of the core (Eq. (11)). If the mass radius scaling were to behave differently, the cooling timescale would no longer be independent of core mass. After the quantity we define as the cooling time has elapsed, the core and envelope can cool together, and one can subsequently employ existing evolution models (e.g., Ginzburg et al. 2016; Vazan et al. 2018; Chachan & Stevenson 2018).

The fact that for is an unphysical consequence of our treatment of the critical point as constant. A more accurate model must consider a composition-dependent critical curve, see Figs. 2 and 4. A more accurate model requires coexistence data of critical mixtures, which will tend to increase the critical temperature and pressure compared to a pure substance. We therefore urge readers to focus more on the concepts here than the numbers, as more data are needed.

We must briefly inspect the assumption that we can treat the envelope as initialized in the equilibrium state. An envelope 1% the mass of the core contains about 10% of the molecules. Therefore if the envelope were warm at early time, then the envelope could cool to a lower luminosity state on timescales about an order-of-magnitude shorter than the time it would take to cool the core and envelope simultaneously. Furthermore, if the initial state of the envelope is higher entropy than our identified equilibrium state, then thermal transport from the core will be tremendously inefficient across the stably stratified envelope (we quantify the early heat transport efficiency in the following paragraph). Therefore our estimates on cooldown timescales from Sect. 4 should be considered as lower bounds, because it may take even longer to cool if the envelope must first cool down to the equilibrium state.

At early times, the core will be hot and the internal heat flux of the planet high. We must therefore investigate the core cooling rate before the planet reaches global radiative-convective equilibrium. Assuming an initially adiabatic atmosphere, for example due to adiabatic compression when accreting isentropic disk gas, and if the envelope is initially saturated everywhere, then the bottleneck for core cooling will initially be the radiative heat loss due to an adiabatic gradient. The adiabatic temperature gradient under the ideal gas approximation for hydrogen will be ∇T ~ 10⁻⁵(M_c/M_⊕)^1/2 K cm⁻¹. Using our order-of-magnitude assumptions for opacity in the pressure regime of interest from Sect. 3.2, we can place constraints on the early heat flow to be on the order of 10⁰–10² erg cm⁻² s⁻¹. Comparing these values to Fig. 6, we see that the heat escape from the core at early times is much lower than the heat flux after achieving global radiative-convective equilibrium of the envelope, consistent with our assumptions.

5 Future Work

In this section, we highlight natural consequences of this new model for super-Earth internal structure. Each of these considerations are fruitful topics for follow-up works, and provide additional motivation for the acquisition of further experimental data. In Sect. 5.1, we discuss a major topic completely sidestepped in this work: the fact that supercritical hydrogen and silicates are expected to be miscible in all proportions. This may have important consequences for internal structure and evolution. Finally in Sect. 5.2, we comment on the importance of coupling our findings to a fully self-consistent and physics-based formation model. Each of these topics will be addressed in more detail in follow-up papers to this series.

5.1 Hydrogen Pollution in the Core

Recent formation models of super-Earths indicate a high metallicity region outside of the embryonic core polluted with hydrogen (Brouwers & Ormel 2020; Bodenheimer et al. 2018). Above the critical temperature and pressure, silicates and hydrogen can mix in all proportions. The exact nature of the mixture depends on thermodynamics, and the formation scenario. During formation, the mixing ratio depends on the ambient thermodynamic conditions (Fig. 2). At low temperatures, there will be two coexisting phases, one of which is nearly pure hydrogen-helium (some trace silicate vapor), while the other is nearly pure silicates (as droplets with trace dissolved hydrogen). Approaching the critical curve at high pressures, we expect the mixing ratio to approach 50%. Indeed according to Fig. 2 there is reason to think the stable layer should truncate at a temperature containing significant hydrogen enrichment. Therefore thermodynamically, significant hydrogen pollution of the core is allowed, but the amount and distribution depend on its formation conditions. The idea of hydrogen pollution into a magma core due to magma dissolution governed by Henry’s law (Chachan & Stevenson 2018), or due to a fugacity crisis (Kite et al. 2019) has already been studied in detail; in this work we suggest hydrogen pollution can be even more efficient, because hydrogen and silicates can mix in all proportions as a supercritical mixture.

The effect of hydrogen in the core depends on its distribution and quantity, but the direction of the effect is clear: greater hydrogen enrichment will inflate the core size relative to pure rock. Hydrogen pollution can also provide a buffer against atmospheric mass loss (cf. Chachan & Stevenson 2018), and an additional source of heat from gravitational settling as the core gradually degasses when it cools. Depending on its distribution, hydrogen may also produce Ledoux-stable compositional gradients that further insulate the core’s primordial heat (cf. Scheibe et al. 2021).

Understanding thermal evolution in detail is impossible without a better understanding of the coexistence of hydrogen and silicates under the relevant thermodynamic conditions. Nevertheless, the timescale for cooling can be estimated to first order from Fig. 8. The cooling of a hydrogen-polluted core will proceed more slowly than suggested by Fig. 8 for two reasons: first, because as the core degasses the RCB will become deeper, the luminosity will drop further, therefore cooling will become even less efficient. Second, the phase separation and settling of higher-density material at depth produces additional heat from gravitational energy. A detailed model of this generalized evolution will be the subject of a future work. However, we must emphasize that such an endeavor is currently premature. A detailed thermal evolution model done properly relies on a better understanding of the physical properties of hydrogen-silicate mixtures under the relevant conditions. The stable layer cannot be modeled properly without accurate thermal transport properties of the silicate/hydrogen mixture. The behavior of the core cannot be properly modeled without accurate supercritical compressibility data (for example, whether the adiabat is steeper or less steep than the critical curve). Thermal evolution cannot be properly modeled without accurate data on the coexistence curve for hydrogen/silicate mixtures. It is for these reasons that in this work we restrict ourselves to qualitative commentaries and order-of-magnitude estimations, rather than employing deceptively detailed dynamical models that rely on unsubstantiated assumptions. Additional ab intio quantum simulations of these materials under the relevant conditions are of great interest for understanding super-Earths.

5.2 Formation and Mass Loss

Our analysis so far has focused exclusively on the planet after formation, and assumes conservation of mass. We must include a brief discussion of the dependence of these results on formation conditions. The consequential elements to consider are the initial temperature of the core, the initial placement of the RCB, core cooling at early times, atmospheric loss at early times, the distribution of silicate vapor in the envelope, and possible pollution of hydrogen in the core. We will discuss each of these points individually.

The upper bound for the initial temperature of the core is the accretion energy. The gravitational binding energy is . Under our assumptions for the relationship between mass and radius for silicate cores, this corresponds to a temperature of approximately (24)

well above the critical temperature for silicate super-Earths. However, such extreme temperatures imply extremely efficient cooling according to the Stefan-Boltzmann law. Therefore, a detailed formation model is needed in order to determine the actual core temperature at early times. Prior modeling of super-Earth planets shows the core temperature can be maintained at high ~10⁴ K temperatures at early times after the emplacement of the hydrogen envelope (Bodenheimer et al. 2018; Brouwers & Ormel 2020). We point to these results from detailed formation studies to justify our assumption that the core can exist at or above the critical temperature of ~6.6 × 10³ K (Xiao & Stixrude 2018) at early times (see also Fig. 4). However, formation models suggest initial core temperatures actually in substantial excess of the critical endpoint temperature, indicating that some initial cooling must occur before before reaching the equilibrium state described in Sect. 3.

Indeed, if the core is initially hotter than the critical temperature, or if hydrogen accretions alongside silicates for some portion of the formation process, then this may lead to a situation in which the stable layer truncates according to our arguments in Sects. 2.2 and 2.3, and possesses an extended layer of hydrogen-enriched but silicate-dominated material. Such a scenario is not precluded by our previous analysis, and the basic arguments still apply. However, one must be careful about terminology. For example, the total mass of gas in the planet may then exceed what we have been calling the envelope mass. Furthermore, evolution would not proceed with the thermal structure of the envelope frozen; rather, the envelope would increase in mass as a function of time. In this case, the cooling timescales from Fig. 8 should be thought of as a lower bound, since this scenario includes additional heat sources from further gravitational settling, and longer cooling timescales as the envelope mass increases with time.

The RCB at early time will be determined by the disk conditions during formation. Formation models find this value to initialize near roughly 10 bars (Bodenheimer et al. 2018; Brouwers & Ormel 2020), varying somewhat depending on specific formation conditions. All equilibrium atmospheric profiles must have an RCB at or beneath its formation value, as is the case for higher-mass envelopes.

We further did not consider atmospheric loss at early times, a highly complex and poorly understood phenomenon that is likely a diagnostic phase of planetary formation and early evolution. This phenomenon has been investigated in detail in prior works, both atmospheric loss by photoevaporation (Owen & Wu 2017) or by core luminosity (Gupta & Schlichting 2019). These dynamics are beyond the scope of this analysis. We are not in this work performing detailed evolution models, but rather considering a planet in equilibrium in a snapshot in time. We do comment, however, that the core powered mass loss mechanism may be affected when considering convective inhibition, as heat flow from the core may be substantially reduced. We also note that hydrogen sequestered in a supercritical core can be a substantial source of additional hydrogen to buffer against envelope mass loss.

Our model assumes the envelope is saturated in silicate vapor, and that the core is pure silicates. The assumption of a saturated envelope is supported by formation models that explicitly consider silicate vapor (Brouwers & Ormel 2020; Bodenheimer et al. 2018). The mechanism for this is the simultaneous accretion of hydrogen and silicates, where at early times accretion is dominated by silicates with some gas pollution, and at later times dominated by gas accretion with silicate pollution ranging from pebbles to planetesimals vaporizing in the atmosphere before cooling and raining out, ensuring a saturated envelope. These formation models also find a so-called outer core polluted with hydrogen as a consequence of the stage of formation when some hydrogen is accreted alongside the continued accretion of silicates. At later evolutionary phases, this outer core may condense and rain out magma, effectively growing the inner, embryonic core (Brouwers & Ormel 2020). According to our model, we expect a polluted outer core to gradually undergo phase separation as described in Sect. 5.1.

6 Discussion

We have presented a new model for the interior and evolution of super-Earth class planets, likely to be the most common planets in the universe (Silburt et al. 2015). Here we have another application of convective inhibition that has been in the literature for decades (Guillot 1995) and has generated renewed interest in recent years (Leconte et al. 2017; Friedson & Gonzales 2017; Markham & Stevenson 2021). In this work we extend the arguments to apply them in the limit where the condensing species rather than the dry gas dominates in abundance. We find an extreme case of convective inhibition, wherein a hydrogen atmosphere with a layer of static stability can effectively insulate a core at very high temperatures (~10⁴ K) for geologic time, or potentially longer than the age of the universe.

This model should be thought of as exploratory, and further work is necessary. To begin with, more work needs to be done to demonstrate the general mechanism of convective inhibition in either a laboratory setting or a high-quality physics-based simulation as an emergent phenomenon. Details about our assumptions, including the assumption of saturation and thermodynamic equilibrium, must be assessed in a dynamical context. Furthermore, this model must be coupled with a realistic formation scenario, as the subsequent evolution may be sensitive to the planet’s initial conditions. The stable layer, which for low heat flows may not be extremely thin, must also be considered explicitly along with relevant thermal transport properties. Furthermore a serious investigation must explicitly account for the possibility of hydrogen pollution in the core, likely to be substantial based on both formation and thermodynamics arguments. Finally, we recapitulate that all of this analysis must rely on additional data on the critical coexistence of hydrogen and silicates, for example from ab initio quantum simulations.

Our model is of interest for the following reasons. First, we predict the luminosity of these planets will be low even at early times, plausibly between the luminosity of Earth and Neptune or less, and should be initially rather insensitive to the planetary system’s age for hundreds of millions to billions of years. This can potentially be tested via the direct imaging of father out exo-planets in the mass range of interest (e.g., Linder et al. 2019; Morley et al. 2017). Second, a supercritical core with a very low internal heat flow may have consequences for the dynamo and expected magnetic field of the planet. Third, the temperatures within the core may be so extreme that the contribution to the planet’s density due to heat cannot be neglected, and the required quantity of hydrogen to match observed exoplanet radii may be smaller than usually thought (cf. Bodenheimer et al. 2018). Along the same lines, sequestered hydrogen may further contribute to core inflation. Fourth, the bimodal distribution of exoplanet radii is a highly active area of research. This model may be of interest to researchers seeking to explain that observation, because current atmospheric escape models rely on interior models that neglect convective inhibition and the possibility that enormous reservoirs of hydrogen are permitted to dissolve into a supercritical core. Finally, this model for super-Earth interiors further demonstrates the importance of considering exoplanetary systems holistically using intuition arising from fundamental physics rather than analogs to more familiar worlds. Indeed, it is yet another reminder that our own planet may be atypical.

Acknowledgments

This work has been primarily funded by NASA FINESST grant number 80NSSC19K1520, as well as funding from the CNES postdoctoral program. We would also like to thank Yayaati Chachan for many insightful conversations, and the reviewer for highly constructive feedback.

References

All Figures

Fig. 1 Cartoon representation of this paper’s findings. On the left is a representation of the canonical view of super-Earths (e.g., Ginzburg et al. 2016), with a frozen or molten silicate core in good thermal contact with a convective gas envelope. On the right we show our proposed model, with a layer of static stability between a (possibly hydrogen-polluted) supercritical fluid core and the silicate-vapor-enriched gas envelope.

In the text

Fig. 2 Illustration showing the stable layer boundaries in (x, T) space. Solid curves use data (Seward & Franck 1981, see Sect. 2.3), while dashed curves use the ideal gas approximation x = ps/p, diverging near the critical endpoint. Gray shaded areas are regions of parameter space wherein free convection is permitted. In the top-right corner, the gray region lies beyond the critical curve where phase separation does not occur. At the bottom of the figure, the gray region lies below the xinh curve described in Eq. (2). Blue, magenta, and red curves represent hydrostatic pressure environments of 300 bars, 1 and 2.5 kbars respectively. The solid points where the dashed curves intersect correspond to the top of the stable layer under the ideal gas approximation; the red asterisk shows the offset intersection using data. The open circles where the solid colored curves intersect the critical curve represent the top of the core.

In the text

Fig. 3 Temperature at the boundaries of the stable layer as a function of envelope basal pressure, using water. The blue and red curves correspond to the top and bottom of the stable layer respectively. As in Fig. 2, we use solid curves where we have data and dashed curves for the ideal gas approximation. We further use a dotted curve for the interpolated critical curve and for a notional (non-quantitative) altered inhibition curve that accounts for the non-ideal coexistence behavior of the water–hydrogen system. The intersection marker symbols from Fig. 2 (the blue, magenta, and red filed and open circles and the red star) are shown here for context.

In the text

Fig. 4 Examples of pressure-temperature profiles for planets in our Solar System and galaxy. If the planet mass is dominated by an ocean (infinite reservoir) with a low-mass hydrogen envelope, convection will be inhibited in the shaded region, due to water in the blue region and silicates in the red region. While the red region uses the ideal gas approximation, the blue region extends into higher pressure regimes; we use a color gradient and question mark to emphasize the behavioral uncertainty under these conditions. Data for Uranus are extrapolated from Lindal et al. (1987), the lupiter model is from Hubbard & Militzer (2016), and the Kepler 36c model shortly after isolation from Bodenheimer et al. (2018). Colored points correspond to the same symbols from Figs. 2 and 3, with added yellow square points corresponding to data from Bali et al. (2013).

In the text

Fig. 5 Envelope pressure–temperature profiles for a sample of envelope properties. The blue and red curves represent the temperature at the top and bottom of the stable layer, respectively; for a saturated gas mixture, convection is shut off between the blue and red curves. In the absence of data, this figure assumes xs = ps/p and does not consider the temperature dependence of latent heat or the compositional dependence of the critical curve. See Fig. 7 for a conceptual sketch of the temperature structure below the envelope.

In the text

	Fig. 6 Core-envelope heat flux as a function of envelope mass. Temperatures in the legend correspond to Teq. Symbols on the left of the plot show the contemporary internal heat fluxes of Earth, Neptune, and Jupiter for comparison.
In the text

	Fig. 7 Sketch of the proposed temperature structures of silicate planets with gaseous envelopes described in Fig. 1. See Fig. 5 for explicitly quantified envelope pressure-temperature diagrams; this figure is only intended to be an illustrative sketch.
In the text

	Fig. 8 Cooling timescales for a super-Earth computed using Eq. (20). The dotted curves show the cooling timescale without convective inhibition with the core temperature at . The solid curves include convective inhibition as described in this work.
In the text