The Runaway Greenhouse on Sub-Neptune Waterworlds

Ever since the pioneering work of Ingersoll (Ingersoll 1969), elaborated on in important ways by Kasting (1988) and Nakajima et al. (1992), the runaway greenhouse phenomenon has played a central role in thinking about planetary evolution. In particular, it is broadly taken as defining the inner limit of the conventional liquid-water habitable zone. Most attention to date has focused on terrestrial planets, defined as planets with a rocky surface, but there is considerable current interest in the class of planets called sub-Neptunes, which may well be the most common type of planet in our galaxy (Bean et al. 2021). Hycean worlds are sub-Neptunes with an H₂-dominated atmosphere over a temperate liquid water ocean and have been proposed as a novel form of habitable world (Madhusudhan et al. 2021). Such planets are essentially a variant of the hydrogen-supported habitable states discussed by Pierrehumbert & Gaidos (2011), with more emphasis placed on application to sub-Neptune waterworlds. The radiative results in Pierrehumbert & Gaidos (2011) have recently been confirmed by Mol Lous et al. (2022), who also show that habitable conditions can persist for a billion years or more. Hycean worlds would succumb to a runaway greenhouse (or never condense an ocean in the first place) if they are above the runaway threshold. Interest is further sparked by the detection of water in the hydrogen-dominated atmosphere of K2-18b (Benneke et al. 2019; Tsiaras et al. 2019), which is estimated to be very near the threshold instellation for runaway.

In this article we address the issue of how a water-dominated sub-Neptune equilibrates if its instellation exceeds the runaway greenhouse threshold and the discontinuous switch in planetary characteristics that is expected for sub-Neptunes below versus above the threshold. The work reported in this article has close affinities with studies of sub-Neptune interior structure (Mousis et al. 2020; Aguichine et al. 2021; Nixon & Madhusudhan 2021) and of atmospheric inflation in the postrunaway state on terrestrial planets by Turbet et al. (2019). Our work has the modest goal of bridging these two classes of study by showing how the interior structure of a water-rich sub-Neptune depends on instellation and on clarifying the way a runaway greenhouse scenario plays out on a sub-Neptune. There has been a great deal of previous work on the interior structure of sub-Neptunes (see Nixon & Madhusudhan 2021 for a fairly comprehensive review of this body of work), and our results do not add particularly to what is already known about the range of interior structures and their effects on the mass–radius relation, apart perhaps from illustrating how the AQUA equation-of-state database (Haldemann et al. 2020) makes such calculations simple. However, little of that work connects the interior structure to the planetary energy budget constraints imposed by a self-consistent radiative-convective atmosphere; without such coupling, a connection cannot be made with the runaway greenhouse phenomenon. The work of Mousis et al. (2020) and Aguichine et al. (2021) on sub-Neptunes with supercritical water interiors is a notable exception. The states treated there are essentially the same as the superrunaway states treated in the present work, but the connection to the runaway greenhouse phenomenon is only mentioned in passing in Mousis et al. (2020) and not explored in depth. The focus here is on clarifying the dramatic transition in interior structure that occurs between planets with lower instellation than the runaway threshold and planets with higher instellation, the nature of equilibration of water-rich postrunaway sub-Neptunes, and the use of the runaway threshold as a simply stated screening criterion for the existence of Hycean worlds. Our results have some observational implications for the mass–radius distribution of the population of weakly irradiated sub-Neptunes spanning the runaway greenhouse threshold.

There are two situations in which coexistence of gaseous and condensed phases restricts the state of the system to the p − T phase boundary between gaseous and condensed water. First, within an atmosphere dominated by the gaseous phase, saturation can lead to the formation of droplets or particles of the condensed phase, and so long as condensed phase continues to form (e.g., through radiative cooling), the system will be restricted to the phase boundary. The phase boundary is calculated by solving the Clausius–Clapeyron relation for temperature as a function of the partial pressure of water p_H2O. This will be called T_df(p_H2O). It is the usual meteorological definition of dew or frost point. For the pure-steam atmospheres, which are the central focus of this article, p_H2O is in fact the total pressure and can simply be written as p. Alternately, if the planet exhibits a layer of essentially pure condensed phase (an "ocean"), then the vapor pressure of the gas phase in contact with the ocean surface is also restricted to the phase boundary, though within the ocean, where there is no coexisting gas phase, the p − T trajectory is not constrained. In this case, the mass of the atmosphere increases with the ocean surface temperature as more of the liquid phase is converted to atmosphere and the ocean becomes shallower.

The essence of the runaway greenhouse phenomenon is that, when the infrared radiating level of the atmosphere lies in the portion where temperature is controlled by T_df(p), increasing the ocean surface temperature cannot increase the radiative cooling to space, so the ocean surface continues to heat until some process intervenes to allow the temperature of the radiating level to increase (Pierrehumbert 2010, Chapter 4). In the conventional picture, the ocean surface heats up until the entire ocean has evaporated into the atmosphere. At that point, there is no more condensed reservoir to feed into the atmospheric mass as the surface heats up, and a noncondensing region governed by the noncondensing ("dry") adiabat forms, extending from the surface to the layer where condensation occurs. As the surface heats up, the dry adiabatic region progressively eats into the bottom of the condensing region, allowing the radiating region of the atmosphere to warm sufficiently to restore energy balance (Boukrouche et al. 2021; at very high temperatures, there is also some contribution to the equilibration by the more widely recognized mechanism of short-wave thermal radiation escaping through atmospheric window regions.)

But what happens if essentially your entire planet is made of water, and it is impossible to "run out" of ocean to evaporate? How, then, does the planet equilibrate in the postrunaway stage? The key is that once the ocean surface temperature exceeds the critical point, there is no longer a distinction between the fluid and gaseous phase, so the atmospheric adiabat connects seamlessly to the supercritical water adiabat that extends into the deep interior of the planet (and possibly connecting to high-pressure fluid or ice phases there). Unlike the dew or frost point adiabat, which has no free parameters once the composition is fixed, the supercritical adiabat (like the conventional dry or noncondensing adiabat) has a free parameter, which can be regarded as its entropy or its temperature at a reference pressure level. To equilibrate, the adiabat needs to heat up to the point where the radiation to space can increase. Equilibration then proceeds as per the usual scenario, primarily by thinning the condensing layer, but also with some contribution from the hot, deep atmosphere radiating through short-wave windows. This picture will be made quantitative in the following.

The central concept of how a runaway greenhouse plays out for a sub-Neptune will be illustrated in the context of a highly idealized set of assumptions regarding the planetary structure. Some ways in which deviations from this idealization might alter our results will be discussed in Section 4. We assume the planet to be completely composed of water, all the way to its center, and to be convective all the way to the center, apart from a radiative stratosphere at the top of the atmosphere. These assumptions constrain the T(p) profile below the stratosphere to lie on the adiabat for water substance in the appropriate phases. In the upper portions of the atmosphere, where the noncondensing adiabat would be colder than the dew or frost point temperature T_df(p), the temperature is set to T_df(p), which defines the phase boundary,

The noncondensing adiabat extending into the interior is computed using the AQUA equation-of-state database for water, which provides tabulations of the adiabatic exponent as well as density (Haldemann et al. 2020). The equation of state covers all likely liquid or ice phases that may be found in the interior.

The depth of the condensing layer is determined by the radiation balance of the planet. Although the analysis in (Boukrouche et al. 2021) was done for a terrestrial planet with a rocky surface, the dominant radiation escaping to space comes from the atmosphere and not from the surface, and so the results can be applied to sub-Neptunes. The atmospheres considered were shallow compared to a typical planetary radius, so the "surface gravity" in this study is essentially the same as the gravity at the level of atmosphere from which radiation escapes to space, which is the relevant "surface gravity" to use in the context of a sub-Neptune, which has no distinct surface anywhere near the radiating layer. For surface gravity 10 m s⁻², the limiting outgoing long-wave radiation (OLR) for a saturated pure steam atmosphere, which defines the runaway greenhouse threshold, is approximately 280 W m⁻² (as read off the flat portion of Figure 8 in Boukrouche et al. 2021). For absorbed stellar radiation (ASR) above this value, the deep atmosphere will continue to warm until the noncondensing adiabat extends sufficiently far into the condensing region to allow the radiating level to warm sufficiently that the OLR increases so as to balance the ASR.

The corresponding gray-equivalent radiating temperature to the runaway threshold flux is 265 K based on the dew-point formula. This corresponds to a radiating pressure of 333 Pa. However, because of radiation through fairly transparent window regions, a considerable amount of flux originates somewhat deeper in the atmosphere than this layer. Based on the surface temperature at which the OLR begins to rise in Figure 8 of Boukrouche et al. (2021) and the corresponding profiles shown in Figure 3 of that work, the actual base of the condensing layer in a slightly superrunaway state (i.e., one for which the ASR slightly exceeds the 280 W m⁻² runaway threshold) is closer to 1000 Pa, and this is the value we shall take as our estimate of the base of the condensing layer for the slightly subrunaway and slightly superrunaway cases. The corresponding temperature on the dew-point adiabat is 280 K and thus slightly above freezing. This point is a key parameter since it determines the pressure at which the p − T trajectory intersects the phase boundary and therefore provides the starting point for continuing the adiabat into the planetary interior. This prescription allows a simplified treatment of the atmospheric radiative transfer, avoiding the complexity of the more precise use of a full radiative-convective calculation such as employed in Mousis et al. (2020). The limiting OLR, upon which our analysis is based, is only weakly (approximately logarithmically) dependent on gravity (Pierrehumbert 2010). Further results on the effect of surface gravity on the runaway threshold can be found in Kopparapu et al. (2014).

T_df becomes arbitrarily cold as p → 0, but the temperature decrease is halted by the formation of a radiative-equilibrium stratosphere above some level, which is warmer than condensation threshold. Stratospheres are maintained both by absorption of upwelling thermal radiation from below and by in situ absorption of incoming stellar flux. As is common practice, we will represent the stratosphere as an isothermal layer though optically thick stratospheres are not generally isothermal. For illustrative purposes, it will be assumed that the stratosphere is thin enough that it does not significantly affect the top-of-atmosphere radiation balance. This can be a poor approximation for very highly irradiated planets, but our focus here is on planets that are near or below the runaway threshold. As the ASR increases, the stratospheric temperature will increase both because the OLR needed to balance the ASR increases (implying an increase in the infrared flux heating the stratosphere from below) and because the in situ heating by stellar absorption increases.

When the ASR is just below the runaway threshold, a liquid ocean can be maintained, and the bottom of the condensing layer is at the liquid ocean interface. In this case the saturated condensing layer is in contact with the surface, and given our estimate of the base of the condensing layer, only a thin water vapor atmosphere would remain, with surface pressure 1000 Pa (or less at lower ASR and lower surface temperature). At greater depths, the adiabat continues into the liquid adiabat, which is not constrained to be on a phase boundary. This situation is depicted in Figure 1(a) for a case where ASR is just slightly below the runaway threshold. When the ASR is above the threshold, a liquid ocean cannot be maintained, and instead the interior of the planet must heat up until water becomes supercritical at some point. In this case, the bottom of the condensing layer is at the boundary of a noncondensing (dry adiabatic) region, as in Boukrouche et al. (2021). This adiabat is subcritical at the point of contact but continues seamlessly into the supercritical fluid adiabat. The situation is depicted for ASR slightly above the runaway threshold in Figure 1(b).

Zoom In Zoom Out Reset image size

The situation described above corresponds most closely to a hot-start case, in which the planet initially forms with a hot interior, which drives a deep convective layer as the heat escapes. The calculation then describes the state the planet attains as the interior cools down sufficiently that the the top of atmosphere comes into balance with ASR though while there is still enough interior heat flux to maintain convection. Superrunaway planets will equilibrate in a supercritical fluid state as they cool, whereas subrunaway planets will eventually condense a liquid water ocean in the interior. The hot-start scenario also applies to the period of declining ASR characteristic of the premain sequence stage of low- mass M and K stars.

The cold-start case, in which the planet starts in a subrunaway condition with a cool interior and liquid ocean but then crosses the runaway threshold owing to luminosity increase of its star, is considerably more complex owing to the difficulty of mixing heat downward against a stable buoyancy gradient. In such cases, the assumption of a fully convective interior is clearly inappropriate.

Once T(p) has been determined, the mass–radius relation can be computed in the usual fashion. The equation of state yields the density ρ(p). The hydrostatic relation yields dp/dr given ρ and the local acceleration of gravity g(r). g(r) is, in turn, given by Newton's law of gravitation, in term of the mass M(r) contained within a shell of radius r, while dM/dr is determined by the density ρ(p). This yields a pair of coupled ordinary differential equations that can be integrated from r = 0 given a guess at p₀, the pressure at the center of the planet. The equations are integrated outward until a specified low pressure p_top is reached. This procedure yields a planetary mass M(p₀) and radius r(p₀), so that the mass–radius relation can be mapped out parametrically by varying p₀. The value of p_top should in principle be determined by the pressure at which the limb of the atmosphere becomes opaque to the wavelength band in which the planet is being observed, but to specify it in this way would become dependent on the nature of the observation. In practice, the computed radius is relatively insensitive to the precise value chosen for p_top because pressure decays approximately exponentially in the outer atmosphere, whence the estimated radius depends only logarithmically on p_top. The resulting radius can be taken as a baseline from which apparent radius as a function of wavelength can be computed in the analysis of transit-depth observations. In results presented here, we use p_top = 1000 Pa. All references to "surface gravity" below refer to the gravitational acceleration computed at this pressure, which is close to the gravitational acceleration at the infrared radiating level.

The curves shown below the stratosphere in Figure 1 are the actual adiabats computed from the AQUA equation of state, with corresponding phases indicated.

The pressure that corresponds to the center of the planet depends on the planet's mass and will be discussed shortly. For slightly subrunaway conditions, there is a deep, cool liquid ocean. This is nearly isothermal, even if fully convective, because the adiabatic exponent for liquid water is very small. At higher pressures, there is a transition to ice X, which has a higher exponent, leading to a layer of increasing temperatures. For lower ASR (not shown), the temperature of the condensed surface goes down, and the surface pressure (marking the lower boundary of the condensing atmosphere) goes down correspondingly. It takes very little reduction in ASR for the ocean surface to freeze, leading to an increase in albedo and further reduction in surface temperature and pressure. This ice-albedo feedback is especially important for the case of G stars but less so for the case of M or K stars whose infrared-rich spectrum leads to lower ice albedo (Joshi & Haberle 2012; Shields et al. 2013; von Paris et al. 2013). A computation based strictly on the AQUA adiabat would lead to an ice-I shell transitioning to higher-pressure ice phases in the deeper interior, but in reality, the rigidity and low thermal conductivity of the shell would inhibit convection, so even a slight interior heat flux would lead to a liquid water ocean under a relatively thin ice shell, as is well known from Snowball Earth calculations (e.g., Li & Pierrehumbert 2011). This is also the state of Europa, so the resulting sub-Neptune could be considered a super-Europa.

From a hot start, the subrunaway equilibrium situation depicted in Figure 1(a) would take some time to emerge because the interior would only gradually cool as the excess heat of formation in the supercritical water deep interior radiates away to space. Liquid or ice layers would only start to accumulate once the interior cooled sufficiently for condensate to persist in the deep interior. At intermediate times, the liquid and ice-X deep interior in the equilibrium state would be replaced by hot supercritical water, leading to a corresponding increase in the radius of the planet relative to the equilibrium state. Modeling the thermal evolution of sub-Neptunes is beyond the scope of the present work.

When the ASR is slightly above the runaway threshold, the gas and supercritical adiabat (Figure 1(b)) have a uniformly large adiabatic exponent, leading to a very hot (≈10,000 K) interior, which will reduce the density and inflate the planet. The contrast between Figures 1(a) and 1(b) indicates the dramatic transition that occurs for sub-Neptunes below versus above the threshold. In the superrunaway case, condensation and consequent formation of water clouds occur only in a thin layer near the top of the atmosphere, even if the ASR only slightly exceeds the runaway threshold. This results from the very nature of the postrunaway equilibration and applies for terrestrial planets as well.

At higher ASR (not shown) the gas/supercritical adiabat intersects the condensing adiabat at lower pressures and somewhat lower temperatures so as to allow the greater radiation from the hot noncondensing layer to balance ASR. Since the reduction of temperature with pressure on T_df is slight, this leads to a hotter adiabat and yet hotter supercritical interior temperatures. For sufficiently large ASR, the increase in stratospheric temperature and deepening of the noncondensing layer would eliminate the condensing layer altogether. This is a familiar situation for highly irradiated sub-Neptunes such as GJ1214b, for which water never condenses (Miller-Ricci & Fortney 2010). As a simple estimate, if one estimates the stratospheric temperature by gray-gas skin temperature, GJ1214b would have an upper-atmosphere temperature of 506 K, well above the temperature of the condensing layer in a postrunaway steam atmosphere.

The radius, surface gravity, and central pressures are given for various planetary masses in Table 1 for the slightly subrunaway and slightly superrunaway cases. As expected, the superrunaway cases are very inflated relative to the subrunaway cases though the proportional inflation decreases with increasing mass as the high-pressure interior in the hot case becomes less compressible. The inflated runaway cases are similar to the runaway greenhouse radius inflation discussed by Turbet et al. (2019), except that they involve the entire interior of the planet and are the same kinds of states as the supercritical water states discussed comprehensively in Mousis et al. (2020). The 300 K equilibrium temperature case in Figure 2 of Mousis et al. (2020) is closest to being comparable to our marginally superrunaway states, and our computed superrunaway radii are reasonably consistent with the results there. For example, for 3M_⊕ our radius is 2.47r_⊕ versus 2.7r_⊕ in Mousis et al. (2020); for 9M_⊕ it is 2.9r_⊕ versus 3.4r_⊕. It is to be expected that the cases in Mousis et al. (2020) are somewhat inflated relative to our marginally superrunaway states since 300 K equilibrium temperature corresponds to an implied ASR of 460 W m⁻², versus the 280 W m⁻² assumed in our marginally superrunaway states.

Table 1 also shows the surface gravity for each case. The runaway threshold is weakly dependent on surface gravity (Pierrehumbert 2010) because the hydrostatic relation implies that for low gravity it takes less pressure at the bottom of the condensing layer to make it optically thick though the dependence is reduced by the pressure-dependence of opacity. Low gravity reduces the ASR threshold for runaway, and high gravity increases it. The superrunaway cases below 9M_⊕ all have lower gravity than the nominal 10 m s⁻² we assumed in the calculation, so in fact they exceed the runaway threshold by somewhat more than assumed in the calculation; accordingly, the interiors would be somewhat warmer than computed, and the planet would be somewhat puffier. A fully self-consistent calculation would need to iterate on surface gravity and ASR threshold, but except for very low masses, the effect is not likely to be great. It is interesting that the runaway actually helps to reinforce itself through reduction in surface gravity. The subrunaway cases can also have reduced gravity for the lowest-mass cases, but by 9M_⊕ the gravity is 40% above the nominal value assumed, which would render the planet somewhat more subrunaway than assumed, yielding a somewhat cooler interior. Similarly to the superrunaway case, the gravity feedback reinforces the state, pushing the planet further into the subrunaway regime.

In Figure 2 we plot the adiabats on the AQUA phase diagram, including a case of a hotter adiabat corresponding to ASR somewhat in excess of the runaway threshold. Note that the boundary between the gas and supercritical regions is not actually a phase transition but only indicates a switch in the equation of state employed in the AQUA compilation. AQUA does not distinguish between supercritical and superionic states, but an examination of the phase diagram in Millot et al. (2019) shows that the superrunaway adiabats are too hot and too low pressure to venture into the superionic ice region. This contrasts with the higher-pressure adiabats starting from colder upper-atmosphere conditions, shown for Uranus and Neptune in Millot et al. (2019). The adiabats shown in the figure are computed out to 1000 GPa, but the pressures given in Table 1 determine where the adiabats actually terminate for planets of a given mass. In all cases, the superrunaway planets terminate in a supercritical fluid region.

Zoom In Zoom Out Reset image size

Our most vulnerable assumption is that the entire atmosphere and interior below the condensing region is convective. The low thermal and radiative diffusivity of the deep interior may mean that it takes very little interior heat flux to maintain convection there, but stellar flux attenuates quite strongly with depth in a steam atmosphere, particularly in the case of M stars (see, e.g., Figure 5.14 in Pierrehumbert 2010). Analytic gray-gas solutions (Guillot 2010; Pierrehumbert 2010) show that in the absence of interior heat flux, the atmosphere asymptotes to a deep, hot isothermal region below the level to which appreciable stellar flux penetrates. The temperature of the radiative layer increases like the fourth root of the infrared optical depth measured at the level to which stellar flux penetrates (Pierrehumbert 2010, Equation (4.56)). Insofar as atmospheres (including water vapor) are much more optically thick in the thermal infrared than at shorter wavelengths, the radiative layer is generally very hot though less so for M stars than for hotter stars. Assuming convection resumes in the deep interior, the isothermal radiative layer would set the interior adiabat to a lower (though generally still supercritical) temperature than the adiabat connecting without interruption to the bottom of the condensing layer. Such states would reduce the inflation of runaway sub-Neptune waterworlds, with the magnitude of the effect increasing as the depth of the radiative layer increases.

The same solutions show that convection can be restored by sufficient interior heat flux if the absorption coefficient of the atmosphere increases sufficiently rapidly with pressure. Nearly isothermal radiative layers are familiar from studies of the atmospheric structure of GJ1214b (Miller-Ricci & Fortney 2010) though the hot upper atmosphere of such highly irradiated planets makes it easier to suppress convection and planets closer to the runaway threshold likely have thinner radiative layers, if any. We do not attempt to estimate the amount of interior heat flux needed to yield a fully convective atmosphere because the threshold is highly sensitive to the opacity behavior of water vapor at high temperatures and pressures. For the most part, these spectroscopic properties—particularly the crucial continua—have been validated only for relatively Earthlike conditions.

According to current observations, K2-18b has an H₂-dominated outer atmosphere, with possibly water as a minor constituent. Given the many ways planets can accrete both H₂ and water, this may be a common situation among sub-Neptunes. Adding a noncondensible background gas that has fairly strong infrared opacity generally brings down the OLR for a given ocean surface temperature but does not change the threshold for a runaway, which is largely governed by the pure-steam limit. For the case of addition of noncondensible CO₂, this can be seen in the OLR curves computed in Wordsworth & Pierrehumbert (2013). Because of its strong collisional opacity, H₂ has a similar effect, as was shown by Koll & Cronin (2019) for the case of a water saturated H₂ − H₂O atmosphere assumed to be on the moist adiabat, though in that case there are thermodynamic as well as opacity effects of the H₂ background. Thus, in a subrunaway case, the addition of a noncondensible greenhouse gas to the atmosphere of an ocean world increases surface temperature and moisture content of the lower atmosphere but does not trigger a runaway. In their analysis of habitability of Hycean worlds, Madhusudhan et al. (2021) incorporated the effects of H₂ and H₂O opacity in the determination of ocean surface temperature, but because the H₂O concentration was assumed fixed (rather than increasing with temperature as it would in equilibrium with an oceanic reservoir), the calculation did not factor in the effect of the runaway greenhouse threshold, in contrast to Koll & Cronin (2019), who allowed for the H₂O source arising from equilibration with an oceanic reservoir. In any event, the increased ocean temperature due to additional opacity sources would lead to a somewhat larger planet relative to the equilibrated all-water subrunaway solutions we have exhibited.

The preceding remarks imply that for a superrunaway case, a planet with a moderately thick H₂ envelope over a massive water ocean would run away until the atmosphere became steam dominated. In a hot-start condition with large water inventory and moderate H₂ inventory, the planet would equilibrate in the supercritical postrunaway state we have exhibited, with the H₂ diluted by the water inventory. "Moderate" in this context means that the planet's H₂ inventory is small compared to its H₂O inventory. The reason H₂ would become diluted into the supercritical water interior is that the dissolution of gases in supercritical fluids is no longer constrained by Henry's law of solubility; rather, the gases become completely miscible with the supercritical fluid (Budisa & Schulze-Makuch 2014). We have been unable to find published experiments specifically on the case of H₂ mixing with supercritical water, but experiments with CH₄ solubility in water show that gases can become nearly miscible with high-temperature water even before the critical point is reached (Pruteanu et al. 2017).

It is of interest to explore the consequences of our results for the state of K2-18b. The fact that K2-18b does not have a steam-dominated outer atmosphere means that either cloud effects have put it into a subrunaway state with a deep liquid ocean, or (if it has superrunaway ASR) it has a small water inventory relative to its H₂ inventory. The 9M_⊕ case in Table 1 is close to the estimated mass of K2-18b. The radius of an all-water subrunaway state is 2.5r_⊕, which only slightly exceeds the observed 2.37r_⊕ radius, but the upper bound on the radius estimate is 2.59r_⊕, bracketing our subrunaway state. If the planet can be kept sufficiently cold through high albedo, the addition of a modest H₂ envelope could dilute the thin water vapor atmosphere and cold-trap water at low levels, yielding the observed H₂-dominated atmosphere. A superrunaway waterworld state at 2.9r_⊕ is well in excess of the observed radius. A silicate core could perhaps bring the radius down enough to match observations, but the superrunaway waterworld state is already ruled out by the observed atmospheric composition, based on our inference that a water-dominated superrunaway state would have a water-dominated outer atmosphere. This line of argument lends further weight to existing observational work tending to rule out the existence of a liquid ocean. Specifically, Scheucher et al. (2020) concluded that for K2-18b, equilibration with an interior pure-water reservoir would yield a high-molecular-weight outer atmosphere and hence a transit-depth spectrum incompatible with the observations. Similarly, Blain et al. (2021) and Charnay et al. (2021) concluded that transit-depth observations are most consistent with a low-metallicity, H₂-dominated atmosphere with water present as a minor constituent. The speculative possibility remains that an as-yet unidentified mechanism could yield high albedo clouds that reflect enough instellation to put the planet in a subrunaway state and permit a liquid ocean and condensed interior, but such a cloud deck would need to be deep enough in the atmosphere so as to not suppress the amplitude of the transit-depth spectrum. In an extensive sensitivity study, Charnay et al. (2021) did not find any indication of a cloud regime that would permit the existence of a liquid water ocean.

Sub-Neptunes may have a silicate/iron core. In principle it would be easy to generalize our calculation by terminating the AQUA adiabats at the core interface, using core radius and mass based on established compressibility estimates. For the hot superrunaway cases, though, the calculation is greatly complicated by the fact that the temperature at the silicate boundary would be high enough to create a magma ocean and associated silicate vaporization, unless the supercritical water envelope is very thin (less than about 10 MPa). Adequate treatment of the problem in that case would require consideration of the effects of supercritical water exchange with the silicate melt and the stabilization of the atmosphere against convection by the high molecular weight of silicate vapor relative to water, such as considered for H₂ envelopes by Misener & Schlichting (2022) and Dorn & Lichtenberg (2021). Incorporation of a moderate silicate/iron core would contract the planet somewhat but not alter our fundamental conclusion that, for water-rich sub-Neptunes, crossing the runaway threshold effects a profound change in the interior state of the planet.

A sub-Neptune waterworld could undergo a cold-start runaway if it cooled off sufficiently, early in its life, to form a condensed interior, but then the host star luminosity increased sufficiently to cross the runaway threshold. This scenario is most likely for G or higher-mass main sequence stars, which have a short premain-sequence stage but become significantly more luminous on the main sequence on a timescale of a few billion years or less. For a cold-start runaway, the atmosphere would still heat up to the depth to which significant stellar flux penetrates. For the cold start, there would be little heat flux escaping the interior and hence nothing to drive convection in the interior. This would result in a hot (and possibly supercritical) isothermal upper layer in contact with a cold liquid or ice boundary. As the isothermal layer radiated into the ocean or ice, it would progressively advance toward the center of the planet through evaporating part of the condensed region and heating it to match the isothermal layer temperature. The energy required to do this would cool the lower part of the radiative layer, and the energy needed to support further advance of the front must be supplied by radiative or molecular diffusion through the deep isothermal radiative layer. The lack of opacity data for hot, high-pressure water poses a severe challenge for quantifying the rate of advance of the isothermal layer toward the planet's center, but the timescale for equilibration of the planet's full interior is certain to be longer than in cases where heat transport is dominated by convection. It is an open question how far the isothermal layer can progress during the ages of observed sub-Neptunes, and this topic requires further study, as does the general issue of thermal evolution of sub-Neptune waterworlds.

We have argued that for water-dominated sub-Neptunes, crossing the runaway greenhouse instellation threshold leads to a profound transition in the interior state of the planet, from one dominated by liquid water and high-pressure ice at moderate temperatures to one dominated by supercritical water. Although our calculations have been performed explicitly only for the case of pure-water composition, we have argued that the fundamental conclusion would not be greatly altered by a moderate silicate/iron core or the presence of a minor proportion of H₂.

The transition is manifest in an inflation of the planet of a given mass when the threshold is exceeded, similarly to the runaway greenhouse radius inflation discussed for more limited ranges of depth by Turbet et al. (2019). For a universe of pure-water sub-Neptunes, one would expect a discontinuous population of planets on the mass–radius diagram for planets below versus above the runaway threshold. Any attempt to confirm this prediction against data will run up against the usual degeneracies in interpretation of the mass–radius data: those associated with possible presence of a silicate/iron core and those associated with presence of a hydrogen envelope. To some extent the hydrogen degeneracy can be ameliorated through observations of atmospheric composition, but this only constrains the outer portions of the atmosphere. We have argued that if a planet has a composition in which hydrogen is a minor constituent relative to water, then in a superrunaway state the hydrogen would be diluted into the massive supercritical water layer, leaving a steam-dominated atmosphere that is incompatible with the observed hydrogen-dominated atmospheres of planets such as K2-18b. Further, our results position the runaway threshold of ASR as a simple screening criterion for the existence of Hycean worlds since any candidate planets that exceed the runaway threshold would evolve into a state with a supercritical water interior (as in Mousis et al. 2020) rather than a state with a habitable liquid ocean. Of the 11 Hycean candidates in Table 1 of Madhusudhan et al. (2021), all but one (K2-18b) is above the runaway greenhouse threshold. For planets not too far above the threshold, though, there is the possibility that cloud effects could suppress the runaway and permit a liquid ocean.

It is of interest to see what candidates are available for investigation of systematic differences between sub and superrunaway sub-Neptune states. The precise instellation threshold for runaway depends on stellar spectral class (which affects albedo), the planet's surface gravity (a weak dependence), and cloud effects (a potentially profound dependence), but as a general guideline, the threshold is in the vicinity of Earth's present instellation. Figure 3 shows cumulative histograms of the instellations of sub-Neptunes in the NASA exoplanet archive; sub-Neptunes are defined for this figure as planets with an estimated radius between 1.5 and 3r_⊕, and instellations were computed from stellar effective temperature, stellar radius, and orbit semimajor axis. Of the roughly 1800 sub-Neptunes in the archive, the overwhelming majority are highly irradiated objects with instellation in excess of 30 times Earth's. These would all be clearly in the superrunaway range, and most would be in regimes like that of GJ1214b where water could not condense at any point in the atmosphere. There are, however, 80 confirmed planets with instellations below twice that of Earth, and this population does straddle the runaway threshold. Only two of these (TOI-2257 b and TOI-2285 b) are TESS detections, but it is to be hoped that TESS will add more candidates to this population in the future. At present, only six of the 80 have measured radii and masses, and three of those only have a weak upper bound on mass (Table 2). The case of K2-18b has already been discussed. K2-3d is superrunaway, but its radius is far smaller than what would be expected of a superrunaway waterworld, so it is most likely a primarily rocky super-Earth, though at 80% of Earth's mean density it most likely has a fairly substantial volatile envelope, which could involve either supercritical water or hydrogen; future atmospheric characterization could help discriminate between the two cases. LHS1140b is subrunaway but has a radius much smaller than the subrunaway pure-water case; it is ultradense and most likely a rock-ball with at most a thin atmosphere or ice/ocean shell. Kepler-138 b has a radius in the range that is generally considered to be likely to be a rocky super-Earth, but its mass estimate of a mere .64M_⊕ firmly positions it as a sub-Neptune, so we have included it in the table. Its instellation puts it in the superrunaway category, unless the planet has an unusually high albedo. We did not include such low masses in Table 1 because the low surface gravity of such planets in the superrunaway state compromises our use of a fixed 280 W m⁻² runaway threshold, but the estimated radius for a marginally superrunaway state using the nominal threshold is 2.3r_⊕. The large radius arises because of the weak gravity associated with low mass. In reality this is an underestimate of the radius since with an estimated surface gravity of just 1.2 m s⁻², the runaway threshold would be considerably below 280 W m⁻², so the planet would be hotter and even more inflated. Even the nominal estimate is nearly double the observed radius, so Kepler-138 b cannot be a superrunaway waterworld with a supercritical interior. If the albedo were large enough (and sufficiently dominant over cloud greenhouse effects) to put the planet in a subrunaway state, then the estimated radius falls to 1.25r_⊕, which is similar to the observed value. However, Kepler-138 is a cool M star, so it would be hard to achieve a high albedo given the weak Rayleigh scattering and low ice albedo associated with M star illumination. It is more plausible that the planet consists of a small rocky core surrounded by a deep, H₂-dominated envelope, but future atmospheric characterization of this interesting target would help to resolve the matter.

Zoom In Zoom Out Reset image size

This does not leave us with any good candidates for runaway versus subrunaway waterworlds though the possibility that K2-3d is a silicate/iron world with an extensive (though not dominant) runaway atmosphere bears further investigation. Hopefully future discoveries and mass measurements will add more candidates that could be fodder for analysis of the runaway transition, but there is a long way to go.

Even though the observational implications of the runaway transition for interpretation of mass/radius data are clouded by degeneracies and (at present) by lack of data, the transition in interior state affects numerous key planetary processes, such as generation of magnetic fields in conducting fluids or ices, and formation of interior magma oceans. Hence, the runaway transition indicates a significant transition in the behavior of the population of sub-Neptunes, which will be increasingly important as more low-instellation sub-Neptunes are discovered and characterized.

This work was supported by the European Research Council Advanced grant EXOCONDENSE, #740963.

Facility: Exoplanet Archive - .