Rocky Worlds Limited to ∼1.8 Earth Radii by Atmospheric Escape during a Star's Extreme UV Saturation

In the past decade, thousands of exoplanet candidates and diverse planetary systems have been found (e.g., Hatzes 2016). The variety of characteristics observed among these planets has raised many questions about planetary formation and evolution. Of particular interest is how low-mass (defined here as less than ∼10 M_⊕), rocky planets form and evolve, given their potential to support habitable conditions (e.g., Forget & Leconte 2014). Central to this question is whether a rocky planet will accrete and retain a massive H₂/He protoatmosphere that represents a non-negligible fraction of the total planetary mass, which may produce uninhabitable surface temperatures, or if such an atmosphere will be lost. These H₂/He protoatmospheres can form on even the smallest planets, because once a protoplanet reaches ∼0.1 M_⊕, it can accrete H₂/He directly from the stellar disk (Hayashi et al. 1979; Lammer et al. 2011; Ikoma & Hori 2012).

If low-mass planets form after the protoplanetary nebula dissipates, they will not accrete H₂/He protoatmospheres (e.g., Massol et al. 2016). However, from planetary formation models, 1–10 M_⊕ planets with periods <100 days may have formed with longer periods when the stellar disk was still present, and migrated to their observed locations through interactions with the disk (e.g., Cossou et al. 2014; Raymond & Cossou 2014). The planets considered in this study have short periods (less than 100 days), and thus may have formed before the stellar disk dissipated. Theoretical calculations suggest that such planets may form with H₂/He protoatmospheres of 1–10 wt.% (Ikoma & Hori 2012; Bodenheimer & Lissauer 2014; Inamdar & Schlichting 2015). Indeed, a number of low-mass exoplanets with thick atmospheres have been observed, indicating there is likely no barrier for such protoatmospheres to form on low mass planets (e.g., Lissauer et al. 2013; Masuda 2014; Cubillos et al. 2016).

Observations and subsequent analysis have shown that rocky planets without thick protoatmospheres are only found up to 1.5–2 R_⊕ in size (Marcy et al. 2014; Weiss & Marcy 2014; Dressing et al. 2015; Rogers 2015), with some recent observations indicating that rocky planets can reach ∼1.9 R_⊕ in size (Buchhave et al. 2016; Demory et al. 2016). If most close-orbiting planets indeed form with thick protoatmospheres, then the least-massive planets must have lost their captured H₂ and He (see Lopez & Rice 2016 for a discussion of atmospheric formation versus atmospheric loss). Thermally driven atmospheric loss in the hydrodynamic escape regime is thought to be able to effectively strip the protoatmospheres from planets up to 5–10 M_⊕ (e.g., Lopez et al. 2012; Owen & Jackson 2012; Lopez & Fortney 2013; Owen & Wu 2013, 2016; Jin et al. 2014; Howe & Burrows 2015; Wolfgang & Lopez 2015; Chen & Rogers 2016; Lopez & Rice 2016; Owen & Morton 2016).

The rapid loss of the protoatmosphere via hydrodynamic escape is driven by the XUV emissions (where XUV is X-ray plus extreme UV) from the host star, which heats the upper atmosphere of the planet. For young, Sun-like stars, this XUV flux can be orders of magnitude larger than the modern Sun (Lammer et al. 2014; Johnstone et al. 2015). A saturated XUV flux can last for ∼100 Myr (Ribas et al. 2005; Jackson et al. 2012; Lammer et al. 2012). The XUV-driven hydrodynamic escape of a protoatmosphere will occur largely during this saturation time, after which the XUV flux and XUV-driven hydrodynamic escape decrease exponentially. While Sun-like stars can erode atmospheres from closely orbiting planets (0.1 au for this model), the modern Earth is orbiting at a sufficient distance that its hydrogen-poor atmosphere is not subjected to hydrodynamic escape (see Catling & Kasting 2017 p. 175 for a discussion of the topic). However, evaporation of Earth-like planets may occur at orbital distances similar to the modern Earth, for some low-mass planets, via water vapor photolysis and subsequent hydrogen escape (e.g., Kasting et al. 2015; Luger & Barnes 2015).

Several recent numerical studies on XUV-driven hydrodynamic escape from low-mass planets have results that overlap, but differ from each another when simulating the loss of protoatmospheres. Lopez & Fortney (2013), using the model of Lopez et al. (2012), showed that planets less than ∼4 R_⊕ could easily lose their atmospheres from hydrodynamic escape, and follow-on work by Lopez & Fortney (2014) suggested that 1.75 R_⊕ was a likely upper size limit for rocky bodies. Similarly, Wolfgang & Lopez (2015) applied the model from Lopez et al. (2012) to data from the Kepler mission and found planets with radii above 2 R_⊕ should have atmospheres of at least ∼1 wt.%, while planets below 2 R_⊕ should have atmospheres less than 1 wt.%. A parameter study of Howe & Burrows (2015) for XUV-driven hydrodynamic escape found that a cutoff between rocky and gas-enveloped planets occurs between 2 M_⊕ and 10 M_⊕, depending on the model orbital distance, which corresponds to a cutoff of 1.26 R_⊕ to 2.16 R_⊕ for an Earth-like density of 5.5 g cm⁻³.

Other studies have discussed an apparent "valley" or "radius gap" in the distribution of exoplanet sizes. Owen & Wu (2013) showed that XUV-driven hydrodynamic escape from low-mass planets results in rocky planets less than 1.5 R_⊕, gas-enveloped planets above 2.5 R_⊕, and a lack of planets with intermediate radii. Jin et al. (2014) studied several hydrodynamic escape models and found a valley between 1 and 2.5 R_⊕ depending on the orbital distance chosen for their model, below which planets are rocky, and above which planets have thick atmospheres, typically of at least a few wt.%. Similarly, the XUV-driven hydrodynamic escape model of Chen & Rogers (2016) found a valley in the range of 1–2 R_⊕.

In this study, using an XUV-driven hydrodynamic escape model in which only the atmospheric mass changes over time, we look at the likelihood that XUV-driven hydrodynamic escape can reproduce the observed ${1.62}_{-0.08}^{+0.67}\,{R}_{\oplus }$ cutoff from Rogers (2015), and seek to examine the dominant factors that lie behind the cutoff. We do so by running our model with parameter ranges that describe the most escape-vulnerable planets studied by Rogers (2015). From these escape-vulnerable planets, our model provides an upper limit on the atmospheric loss rate, and thus the radii of planets that can lose their entire protoatmospheres and become rocky.

We calculate the cutoff between rocky and gas-enveloped planets with 10,000 different model parameter combinations. We consider only Sun-like stars in this work, because the planets used in the study of Rogers (2015) all orbited stars with effective temperatures between 4700 and 6300 K. The atmospheric loss from XUV-driven hydrodynamic escape around M dwarfs has been considered elsewhere (e.g., Tian 2009; Luger et al. 2015), and we do not address such systems here.

During hydrodynamic escape, a high-altitude portion of an atmosphere is heated by XUV flux and flows hydrodynamically outward (Mordasini et al. 2012; Johnstone et al. 2015). However, for thick protoatmospheres that likely represent at least a few wt.% of a planet (Ikoma & Hori 2012; Bodenheimer & Lissauer 2014; Inamdar & Schlichting 2015), the bulk of the lower atmosphere will approximately remain in hydrostatic equilibrium. As such, we assume for this model that, above the XUV absorption level, R_XUV, where the optical depth for the XUV is near unity, the atmosphere is in the hydrodynamic regime, and below R_XUV, the atmosphere is in approximate hydrostatic equilibrium.

The energy-limited rate of XUV-driven hydrodynamic escape from a planet can be approximated by a first-order equation, as follows

$$\begin{eqnarray}&&\displaystyle \frac{{dM}}{{dt}}=\displaystyle \frac{\eta \pi {F}_{\mathrm{XUV}}{R}_{\mathrm{XUV}}^{3}}{{{GM}}_{{\rm{p}}}},\end{eqnarray} \tag{ 1 }$$

where dM/dt is the rate of hydrodynamic escape in kg s⁻¹ (Watson et al. 1981). The parameter η in Equation (1) is an efficiency factor that is typically taken to be 0.1 < η < 0.6 (e.g., Lammer et al. 2013; Owen & Wu 2013; Koskinen et al. 2014; Bolmont et al. 2017). The XUV flux incident on the planet in W m⁻² is given by F_XUV, R_XUV is the radial distance from the planetary center at which the optical depth for broadband XUV radiation is unity, G is the gravitational constant, and M_p is the mass of the planet. For the range of F_XUV values considered in this study, the rate of hydrodynamic escape may border on the recombination-limited regime described by Murray-Clay et al. (2009). In the recombination-limited regime, for large XUV fluxes (greater than ∼10 W m⁻²), the protoatmospheres could lose energy via recombination of ionized gas slowing the hydrodynamic loss rate. In this recombination-limited regime ${dM}/{dt}\propto {F}_{\mathrm{XUV}}^{1/2}$, while dM/dt ∝ F_XUV in the energy-limited regime. However, on small planets, the rate of hydrodynamic escape is dominated by the R_XUV term in Equation (1), so the difference between the recombination-limited and energy-limited regimes is likely small. As such, we consider only the energy-limited case in this study. To find dM/dt, we need only determine the XUV flux and R_XUV.

The XUV flux from young FGK stars is largest for the first 100 ± 20 Myr after formation (Ribas et al. 2005; Jackson et al. 2012; Lammer et al. 2012). Emissions of XUV are saturated during that time and remain approximately constant. Afterward, the XUV flux diminishes exponentially and the hydrodynamic loss rate of a planetary atmosphere drops with it. During the saturated regime, for a Sun-like star, the XUV flux can reach ∼0.1% of the bolometric luminosity (Jackson et al. 2012; Lammer et al. 2014). Given the uncertainty of stellar evolution, a Sun-like star could generate 43–172 W m⁻² in the XUV at 0.1 au following Pizzolato et al. (2003). For comparison, at 1 au, the present Earth receives an XUV flux of only ∼5 mW m⁻² (Lammer et al. 2014).

Around Sun-like stars, the protoatmospheres of some low-mass planets orbiting interior to 0.1 au will likely extend beyond the planet's Roche lobe and be rapidly lost (Owen & Wu 2013; Ginzburg & Sari 2017). Not only would the atmosphere be rapidly removed by the high flux and gravitational effects of the host star interior to 0.1 au, but the rocky core could begin to evaporate as well (Perez-Becker & Chiang 2013). The planets used in the Rogers (2015) work fall mostly outside this limit, with only four planets, Kepler-10b, Kepler-21b, Kepler-98b, and Kepler-407b, receiving a flux greater than the 0.1 au equivalent around a Sun-like star. Kepler-10b, Kepler-21b, and Kepler-407b are less than ∼1.6 R_⊕ in radii. In contrast, Kepler-98b has a radius of 1.99 ± 0.22 R_⊕, with a density of 2.18 ± 1.21 g cm⁻³, so it is not a rocky planet (Marcy et al. 2014). Of the three rocky planets, none exceed the radius limit described by Rogers (2015). Thus, we use an orbital distance of 0.1 au as the inner bound for rocky planets in our model.

We assume that absorption of XUV occurs downward through an upper atmosphere, and is fully absorbed by the base of a thermosphere. The problem of finding R_XUV then becomes a matter of finding the radial distance to the base of the thermosphere. For a neutral H₂ atmospheric column, the broadband XUV flux is typically absorbed within a column of density 10²⁶ m⁻² (e.g., Glassgold et al. 2004; Cecchi-Pestellini et al. 2006; Ercolano et al. 2009; Owen & Jackson 2012). For an Earth-mass planet, with gravity assumed constant, this corresponds to a pressure at the base of the thermosphere of p_XUV = 3.3 Pa and p_XUV = 7.1 Pa for a planet of 10 M_⊕. However, on the modern Earth, the base of the thermosphere can occur at pressures as low as p_XUV = 0.1 Pa (Catling & Kasting 2017, p. 4) so we will consider a range of pressures from 0.1 ≤ p_XUV ≤ 10 Pa. Once p_XUV is known, it remains fairly constant over a planet's lifetime, even if the surface pressure changes by orders of magnitude (Erkaev et al. 2013). Thus, as rocky planets lose their substantial protoatmospheres, p_XUV does not change, but rather moves closer to the planetary surface as surface pressure drops.

The protoatmosphere of a young planet will be in approximate hydrostatic equilibrium from the surface to the base of the thermosphere. The radial distance to the base of the thermosphere, R_XUV, can then be found via the hydrostatic equation, which can be written as

$$\begin{eqnarray}&&-{R}_{{\rm{g}}}T\displaystyle \frac{1}{p}{dp}={g}_{{\rm{s}}}{\left(\displaystyle \frac{{R}_{{\rm{s}}}}{R}\right)}^{2}{dR}\end{eqnarray} \tag{ 2 }$$

for specific gas constant R_g, pressure p, isothermal temperature T, surface radius of the rocky core R_s, and surface gravity g_s. In Equation (2), we have approximated the gravity term of the hydrostatic equation as g(R) = g_s [R_s/R]² (see Catling & Kasting (2017), p. 11 for a discussion of the topic). In this study, we assume all rocky planets form with an Earth-like structure and relate mass to radius via

$$\begin{eqnarray}&&{R}_{{\rm{s}}}=1.3\,{M}_{{\rm{p}}}^{0.27},\end{eqnarray} \tag{ 3 }$$

which provides a good approximation for the rocky cores of planets in our model. Equation (3) is derived from the relation ${R}_{{\rm{s}}}\propto {M}_{{\rm{p}}}^{0.27}$ given by Zeng et al. (2016). To make this relationship hold when using SI units for the mass (5.9742 × 10²⁴ kg) and radius (6.371 × 10⁶ m) of the Earth, a scalar value of 1.3 must be used, as seen in Equation (3). Integrating Equation (2) from the surface to p_XUV and solving for R = R_XUV we find

$$\begin{eqnarray}&&{R}_{\mathrm{XUV}}=\displaystyle \frac{{R}_{{\rm{s}}}^{2}}{H\mathrm{ln}({p}_{\mathrm{XUV}}/{p}_{{\rm{s}}})+{R}_{{\rm{s}}}},\end{eqnarray} \tag{ 4 }$$

with scale height H = R_gT/g_s. Surface pressure, p_s, can be defined as

$$\begin{eqnarray}&&{p}_{{\rm{s}}}(M)=\displaystyle \frac{{g}_{{\rm{s}}}(\alpha {M}_{{\rm{p}}}-M)}{4\pi {R}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 5 }$$

where α is the initial protoatmospheric mass fraction and M is the integrated mass loss from Equation (1). We assume an initial α for our model, leaving only the scale height, H, unknown.

The atmospheric scale height depends on the temperature profile and atmospheric composition. We are interested in the upper limit on rocky planet radii, so the atmospheric composition was assumed to be pure H₂, with R_g = 4157 J kg⁻¹ K⁻¹as an upper limit. This provides an upper bound on H, and thus on the loss rate and radii limit. The atmospheres of Uranus and Neptune have specific gas constants of R_g ≈ 3600 J kg⁻¹ K⁻¹ (Lunine 1993) and may represent a composition similar to the protoatmospheres we are modeling, so we consider the range 3600 ≤ R_g ≤ 4157 J kg⁻¹ K⁻¹. To calculate the scale height with R_g, we use an isothermal atmospheric temperature.

For the protoatmospheres in this study, which represent up to a few wt.% of the total planetary mass, ∼3000 K is a reasonable upper limit on surface temperature, and thus the isothermal upper atmospheric temperature, based on temperature profiles from more complex calculations (Mordasini et al. 2012; Jin et al. 2014). In addition, on larger planets with H₂ dominated atmospheres, cooling from gas expansion and Lyα radiation in the upper atmosphere likely results in temperatures of ∼3000 K in the thermosphere (Murray-Clay et al. 2009). We expect similar processes to occur in the protoatmospheres of low-mass planets. With 3000 K being a reasonable upper limit at both the planetary surface and the base of the thermosphere, we set the upper limit for the isothermal atmospheric temperature in our model to 3000 K. For a lower bound, we set the isothermal atmospheric temperature equal to the effective temperature at 0.1 au. For a Bond albedo of 0, this distance corresponds to an effective temperature of 880 K around a Sun-like star.

The atmospheric loss rate is easily calculated from Equation (1). First, an orbital distance is chosen (0.1 au in this model), as well as a planetary mass. To calculate the atmospheric loss rate, and thus the cutoff between rocky and gas-enveloped planets, seven additional model parameters must be specified. These parameters are: isothermal atmospheric temperature, T; XUV flux at the given orbital distance, F_XUV; escape efficiency, η; initial atmospheric mass fraction, α; pressure at the base of the thermosphere, p_XUV; specific gas constant of the atmosphere, R_g; and XUV saturation time, τ. We consider 10,000 combinations of these seven parameters, which were selected randomly from a uniform distribution from the values in Table 1. The result for the radius cutoff and its uncertainty is insensitive to the exact number of parameter combinations, whether 10,000, 20,000, or 5000, based on sensitivity tests. For each parameter combination, we calculate the cutoff between rocky and gas-enveloped planets using a time step of 10,000 years, which runs quickly and is as accurate as smaller step sizes. Figures 1(a) and (b) show an example of a model run where the cutoff occurs in the 1.2–1.6 R_⊕ range for the given parameters. Supplemental material provides an animated version of the model shown in Figure 1(b), where one can easily see the how gas-enveloped rocky cores of low-mass planets evolve into dense rocky planets that become bereft of their primordial atmospheres over time.

We examined the protoatmospheric loss from planets between 0.5 and 10 M_⊕. The results can be seen in Figure 2, which shows the counts of the calculated cutoff radii, R_s, between rocky and gas-enveloped planets for the 10,000 random parameter combinations. The red dot shows the observed cutoff radius for rocky planets from Rogers (2015) with a 2σ uncertainty. The black dot and error bar show the model mean and 2σ uncertainty. The distribution in Figure 2 aligns well with the observed rocky planet limit with both the mean and mode falling within the 95% confidence interval of the observed rocky planet limit. While the mean of our model is 1.76 ± 0.38 (2σ) R_⊕, the mode falls closer to ∼1.9 R_⊕, which is where the largest rocky planets have been found. Our model predicts that, beyond ∼1.9 R_⊕, there is a fairly sharp drop-off in the likelihood that hydrodynamic escape can erode a planet. This agrees with recent observations that the largest rocky planets are found up to ∼1.9 R_⊕ (Buchhave et al. 2016; Demory et al. 2016). The largest cutoffs predicted in our model are due to parameter combinations with high isothermal atmospheric temperatures and large specific gas constants. In addition to planetary mass, these two parameters control the R_XUV term in Equation (1), which dominates the loss rate, as discussed below.

Zoom In Zoom Out Reset image size

Upon close examination, we find that a key aspect of the cutoff between gas-enveloped and rocky planets is that the atmosphere remains very distended, up to several planetary radii in size, and available for XUV absorption on low-mass bodies even as atmospheric mass is lost. That protoatmospheres remain puffy, even at low mass, is a result of the logarithmic term in Equation (4) generating large values for R_XUV until the atmosphere is completely removed. This is seen in Figure 3, where we show R_XUV of a 2 M_⊕ planet over time. Even when less than 20% of the original protoatmosphere remains at 20.3 Myr, the radius R_XUV is roughly twice the radius of the core. The large radius is caused by the lower gravity on low-mass planets, coupled to high temperatures (up to several thousand Kelvin) and light atmospheric compositions of H₂/He that lead to substantial scale heights and R_XUV values that increase rapidly with decreasing mass. In Equation (1), we see that the R_XUV term is cubed, and it is the only term that changes greatly as mass is lost (for a 3 wt.% atmosphere M_p will only change by at most 3%). Thus, any change in the loss rate will be dominated by the ${R}_{\mathrm{XUV}}^{3}$ term.

Zoom In Zoom Out Reset image size

In Figure 4, we see how, for planets with mass <2.5 M_⊕, ${R}_{\mathrm{XUV}}^{3}$ is orders of magnitude larger than R_⊕ while the protoatmosphere remains, but for planets with masses of 6–7 M_⊕, ${R}_{\mathrm{XUV}}^{3}$ levels off to ∼17 R_⊕. Not only do low-mass planets have much larger loss rates due to this exponential increase in R_XUV, they also have less overall atmospheric mass to lose. The strong nonlinearity of the hydrostatic equation shows us that there will exist a critical planet size below which R_XUV increases rapidly, leading to substantial hydrodynamic escape. This nonlinear dependence on R_XUV has been noted in previous work (e.g., Lopez & Fortney 2013; Chen & Rogers 2016), and on average, in our model, it results in planets with cores larger than 1.76 ± 0.38 (2σ) R_⊕ retaining a significant portion of their protoatmospheres. The largest planets that can lose their entire protoatmosphere are thus planets with rocky cores less ∼1.8 R_⊕.

Zoom In Zoom Out Reset image size

In addition to the predicted cutoff at 1.8 R_⊕, our model shows that there should be a lack of planets with radii immediately larger and smaller than the cutoff radius. This is seen in Figure 1, where the abrupt jump from rocky to gas-enveloped planets may result in a void where planets are unlikely to exist. This agrees with a number of previous studies (see Introduction) and the recent work by Fulton et al. (2017), which found that such a deficit is indeed present in the current exoplanet data.

The exoplanet data from Fulton et al. (2017) is shown in Figure 5, with our model predictions. Our model indicates that, for radii below ∼2 R_⊕, planets are less likely to be gas-enveloped (the blue region in Figure 5), and for radii above ∼1.5 R_⊕, planets are less likely to be rocky (the red region in Figure 5). The paucity of 1.5–2 R_⊕ planets predicted by our model is seen in the Fulton et al. (2017) data. The missing planets fall into the evaporation valley, or radius gap described by previous XUV-driven hydrodynamic escape studies (see Introduction) and show that XUV-driven hydrodynamic escape is able to reproduce the major characteristics of the observed, low-mass exoplanet population.

Zoom In Zoom Out Reset image size

The transition from rocky to gas-enveloped planets occurs where XUV-driven hydrodynamic escape predicts such a transition should occur (see Figure 5). That the transition has been predicted across numerous studies with models of varying complexity (see Introduction) and agrees with current exoplanet data leads us to conclude that hydrodynamic escape is plausibly the cause of the observed limit in rocky planet radii. The closely orbiting exoplanets (periods less than ∼100 days) modeled in this study comprise the majority of known, low-mass exoplanets (e.g., Batalha 2014). That XUV-driven hydrodynamic atmospheric escape is important for these planets is not surprising, given the large XUV fluxes present at such short orbital periods.

As additional, longer-period, rocky planets are discovered, the average cutoff in rocky planet size may decrease with increasing orbital period, as noted by Lopez & Rice (2016). However, at large orbital distances where the XUV flux is small and XUV-driven hydrodynamic escape becomes negligible, other processes may limit the size of rocky planets. Indeed, Zeng et al. (2017) found a bimodal distribution in the current exoplanet data similar to Fulton et al. (2017), but note that it could be explained by formation scenarios rather than evolutionary ones (i.e., XUV-driven hydrodynamic escape). The radius limit for closely orbiting rocky planets appears to be set at ∼1.8 R_⊕ by XUV-driven hydrodynamic escape, but to address the limit in rocky planet size for longer period planets, additional studies on rocky planet formation should be conducted.

We would like to thank Kevin Zahnle and Joshua Krissansen-Totton for their advice and guidance on this paper. D.C.C. was supported by NASA Planetary Atmospheres grant NNX14AJ45G. O.R.L. was supported by the NASA Pathways Program.