Mass Loss by Atmospheric Escape from Extremely Close-in Planets

Atmospheric escape can have significant consequences for the evolution of planetary atmospheres as well as the structure and bulk density of different planets. For example, it is widely accepted that the early atmospheres of Venus, Earth, and Mars were shaped by escape (e.g., Lammer et al. 2008). Exoplanets offer an extended parameter space that we can use to test models of atmospheric escape. The population of known exoplanets includes features such as the radius valley between rocky planets and planets with a more substantial gaseous envelope (Fulton et al. 2017), and the hot Neptune/super-Earth deserts (e.g., Mazeh et al. 2016; Berger et al. 2018) that are often linked to the loss of gaseous envelopes over time (e.g., Owen & Wu 2017; Gupta & Schlichting 2020). Many transiting exoplanets also exhibit signatures of atmospheric escape in their transmission spectra (e.g., Vidal-Madjar et al. 2003; Linsky et al. 2010; Lecavelier des Etangs et al. 2012; Ben-Jaffel & Ballester 2013; Ehrenreich et al. 2015; Bourrier2018; Spake et al. 2018; Sing et al. 2019) that have inspired new developments in the theory of atmospheric escape.

The implications of atmospheric escape on the envelopes of close-in planets have been debated since the first hot Jupiter, 51 Pegasi b, was discovered (Mayor & Queloz 1995). Since then, observations of transiting planet atmospheres have constrained the escape mechanism and mass-loss rates in theoretical models and revealed qualitative differences between different types of planets and host stars. For example, transit observations in the far-ultraviolet (FUV) obtained by the Hubble Space Telescope probe hydrogen in the escaping upper atmosphere, along with some heavier elements such as carbon, oxygen, and silicon if the latter are present at high altitudes with sufficient abundances. These observations show that many close-in planets lose mass to hydrodynamic escape. They also indicate that Neptune-mass planets such as GJ436b and GJ3470b have larger upper atmosphere transit depths (Ehrenreich et al. 2015; Bourrier 2018) than hot Jupiters such as HD209458b or HD189733b (Vidal-Madjar et al. 2003; Lecavelier des Etangs et al. 2012).

Generally, escape is powered by ionization and heating of the planetary upper atmosphere by stellar X-ray and UV (XUV) radiation. Most models of this process agree that typical hot Jupiters do not lose a substantial fraction of their envelopes to mass loss, even when their atmospheres escape hydrodynamically. Detailed models of the upper atmosphere that have been used to explain transit observations generally predict mass-loss rates of the order of 10⁷–10⁸ kg s⁻¹ for "mainstream" hot Jupiters such as HD209458b or HD189733b (e.g., Yelle 2004, 2006; García Muñoz 2007; Koskinen et al. 2013a, 2013b; Chadney et al. 2017). With these rates, cumulative mass loss from many hot Jupiters amounts to at most a few percent of the current planet mass, depending on the XUV evolution of the host star. This conclusion, however, does not necessarily apply to all hot Jupiters.

Extremely close-in planets such as WASP-12b and WASP-121b have near-UV (NUV) transit depths that substantially exceed the visual transit depth and therefore probe the escaping upper atmosphere (Fossati et al. 2010; Sing et al. 2019). Since hydrogen and helium do not absorb in the NUV, the transit depths are due to heavier elements and metals that escape the atmosphere in large quantities. As might be expected for such planets, this implies exceptionally high mass-loss rates. These planets orbit so close to their host stars, with an orbital period of just ∼1 day, that their envelopes are substantially distorted by the stellar tide. Therefore, they fill a significant fraction of their Roche lobe i.e., the smallest equipotential surface that can exist around the planet in the combined gravity field of the planet and the star. Under these circumstances, the middle atmosphere extends to the Roche lobe, and the resulting high mass-loss rates can have cataclysmic consequences (Lai et al. 2010; Li et al. 2010). In this case, mass loss is powered by the stellar bolometric luminosity and depends roughly on the effective temperature of the planet instead of the stellar XUV radiation that heats the upper atmosphere.

Observations of transiting warm Neptunes (10–25 M_⊕), at Lyα and the He i 1083 nm line, imply mass-loss rates of about 10⁶–10⁹ kg s⁻¹, i.e., of the same order of magnitude as for the known hot Jupiters (e.g., Ehrenreich et al. 2015; Loyd et al. 2017a; Bourrier 2018; Mansfield et al. 2018). While mass-loss rates ostensibly retrieved from observations are actually model-dependent, in the sense that they depend on what physics is included in the model, and on the assumed planet and host star properties, these values are consistent with the idea that the atmospheres of hot Neptunes and lower-mass planets in general are even more likely to undergo hydrodynamic escape than the atmospheres of hot Jupiters (Koskinen et al. 2014). If this is the case, the mass-loss rate is "energy-limited" (e.g., Watson et al. 1981; Erkaev et al. 2007) roughly to within an order of magnitude, and the escape rate depends inversely on the bulk density of the planet, which is similar for hot Neptunes and Jupiters. The relative impact of mass loss on hot Neptunes, however, is larger, especially if one accounts for their evolution over time and the higher stellar XUV fluxes of young stars (e.g., Lopez et al. 2012). Hot Neptunes and lower-mass planets are also more susceptible to mass loss by Roche lobe overflow than hot Jupiters. Lower-mass planets may therefore lose a substantial fraction of their gaseous envelope over time, explaining the occurrence of the hot Neptune/super-Earth desert in the population of known planets (e.g., Owen & Wu 2017).

In this work, we focus on close-in planets that are likely to have the highest mass-loss rates. In particular, we explore the transition from stellar XUV-driven escape of the upper atmosphere to rapid Roche lobe overflow of the lower and middle atmosphere. This transition relates to a fundamental stability limit for gaseous envelopes. We present a relatively simple theory for it, highlight the underlying physics, and show how the equilibrium shape of a planet can be used to assess the stability of its envelope. While Roche lobe overflow is a well-known phenomenon in binary star, exoplanet, and planetary systems in general, the results here point to some new insights to the nature of its application to close-in exoplanets. In binary star systems, for example, Roche lobe overflow occurs when one of the stars fills its Roche lobe and begins to lose mass to the other star. The same is expected to occur in close-in exoplanet systems where the material lost from the planet forms a torus around the star (Fossati et al. 2013), analogous to, say, the Io plasma torus around Jupiter. Near the classical Roche limit, the surface of the planet reaches the Roche lobe, and the planet loses its envelope to cataclysmic mass loss. The visible surface need not fill the Roche lobe, however, for the planet to lose its envelope over a relatively short time (e.g., Jackson et al. 2017). Here, we quantify the limit for significant envelope loss in terms of the tidal perturbation of the planet and extent of its atmosphere, and properly relate it to stellar XUV-driven escape.

We compare our predictions with detailed simulations of atmospheric escape and validate them by comparison with the observed properties of known exoplanets. Since Neptune-mass planets undergo rapid mass loss more readily than hot Jupiters, we focus the detailed simulations on a Uranus-like planet orbiting close-in to a Sun-like star. An additional purpose of these simulations is to identify qualitative differences between planets similar to Uranus and Neptune and previously published simulations of hot Jupiters. This helps to inform approximations of mass loss used in other studies and consistently explain some of the general differences between transit observations of these types of planets. Here, we model a Uranus-like planet because its mass, radius, and overall structure are consistent with and represent a real planet, instead of depending on uncertain models of interior structure and composition. We chose Uranus instead of Neptune because its surface gravitational potential is slightly lower than that of Neptune, resulting in a somewhat more extended atmosphere that might be more typical of close-in planets.

Below, we provide a relatively concise summary of our methods and underlying theory. In Section 2.1, we explain how we include the Roche potential in our models and calculate the shape of the equipotential surfaces for different planets. In Section 2.2, we present a straightforward method for estimating the extent of the lower and middle atmosphere that affects the mass-loss rates. In Sections 2.3 and 2.4, we provide a simple approximation for calculating mass-loss rates due to Roche lobe overflow and review the equations for energy-limited escape that are often used to represent stellar XUV-driven escape (see Owen et al. 2020, for review). In Section 2.5 and Appendix B, we describe a numerical, multispecies model of hydrodynamic escape that we use to test and validate our simple equations of atmospheric escape, in addition to comparing our predictions with the properties of known exoplanets.

2.1. Equilibrium Shape

At close-in orbits, pressure levels in the gaseous envelopes of planets are significantly perturbed by the stellar gravity field. As the orbital distance decreases, the outer layers of the atmosphere and eventually the visible surface of the planet approach the Roche lobe. In hydrostatic equilibrium, isobars coincide with surfaces of constant gravitational potential, provided that horizontal variations in temperature can be ignored. The relevant gravitational potential here is the Roche potential (e.g., Pringle & Wade 1985):

$$\begin{eqnarray}\begin{array}{rcl}U(x,y,z) & = & -\displaystyle \frac{{{GM}}_{p}}{{\left({x}^{2}+{y}^{2}+{z}^{2}\right)}^{1/2}}-\displaystyle \frac{{{GM}}_{s}}{{\left[{\left(x-a\right)}^{2}+{y}^{2}+{z}^{2}\right]}^{1/2}}\\ & & -\displaystyle \frac{1}{2}{{\rm{\Omega }}}^{2}\left[{\left(x-\mu a\right)}^{2}+{y}^{2}\right],\end{array}\end{eqnarray} \tag{ 1 }$$

where z is the distance along a coordinate axis pointing to the north pole from the center of the planet, x axis points toward the star, y axis completes a right-handed coordinate system, G is the gravitational constant, M_p is the mass of the planet, a is the orbital distance, M_s is the mass of the star, Ω is the orbital angular rotation rate, given by Kepler's third law:

$$\begin{eqnarray}&&{{\rm{\Omega }}}^{2}=\displaystyle \frac{G({M}_{s}+{M}_{p})}{{a}^{3}},\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{M}_{s}}{{M}_{s}+{M}_{p}}\end{eqnarray} \tag{ 3 }$$

is the mass ratio. In this study, we assume that planets orbit their host stars in a circular orbit within 0.1 au and are tidally locked (e.g., Guillot et al. 1996). In spherical polar coordinates centered at the planet, the Roche potential is

$$\begin{eqnarray}\begin{array}{rcl}U(r,\theta ,\phi ) & = & -\displaystyle \frac{{{GM}}_{p}}{r}-\displaystyle \frac{{{GM}}_{s}}{\left[{\left(a-r\zeta \right)}^{2}+{r}^{2}{\xi }^{2}+{r}^{2}{\nu }^{2}\right]}\\ & & -\displaystyle \frac{1}{2}{{\rm{\Omega }}}^{2}\left[{\left(r\zeta -\mu a\right)}^{2}+{r}^{2}{\xi }^{2}\right],\end{array}\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&\zeta =\sin \theta \cos \phi ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\xi =\sin \theta \sin \phi ,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\nu =\cos \theta ,\end{eqnarray} \tag{ 7 }$$

where θ and ϕ are the co-latitude and longitude, respectively. We solve for equipotential surfaces numerically by using an iterative method (see Appendix A).

As we will demonstrate, the equilibrium shape of the pressure levels based on the Roche potential can be used to predict the stability of gaseous envelopes. In Section 3, we describe the equilibrium shape by calculating the ratio of the substellar radius to the polar radius (x/z) for different planets. This ratio is obtained by equating

$$\begin{eqnarray}&&U(x,0,0)=U(0,0,z)\end{eqnarray} \tag{ 8 }$$

and solving for x/z. The L1 point r_L1 = x_R is the substellar radius of the Roche lobe equipotential surface, for which x_R /z_R ∼ 1.5. We note that substellar gravity is given by

$$\begin{eqnarray}\begin{array}{rcl}{g}_{s}(r) & = & -\displaystyle \frac{\partial U(r,\pi /2,0)}{\partial r}=-\displaystyle \frac{{{GM}}_{p}}{{r}^{2}}\\ & & +\displaystyle \frac{{{GM}}_{s}}{{\left(a-r\right)}^{2}}+{{\rm{\Omega }}}^{2}(r-\mu a).\end{array}\end{eqnarray} \tag{ 9 }$$

We obtain the radius at the L1 point by setting the substellar gravity to zero and solving the resulting equation for r_L1 by using Newton's method.

2.2. Lower and Middle Atmosphere

In addition to the gravity field, the extent of the lower and middle atmosphere depends on the temperature and composition through the mean molecular weight. State-of-the-art models for the composition and thermal structure of the lower and middle atmosphere include self-consistent treatments of photochemistry/thermochemistry and radiative transfer in the presence of clouds and hazes (e.g., Lavvas & Arfaux 2021). We do not aim to reproduce such complexity here. Instead, we develop a simple approximation to estimate the extent of the lower and middle atmosphere for a H/He envelope to highlight its effect on mass loss. This approximation is roughly valid for a solar-metallicity gas where the mean molecular weight is largely set by the relative abundances of H₂, He, and H. While higher-metallicity envelopes are possible, this assumption is consistent with most current studies of mass loss and exoplanet populations (e.g., Owen & Wu 2017).

We calculate the extent of the atmosphere, or pressure level heights, by numerically integrating the hypsometric equation:

$$\begin{eqnarray}&&{\int }_{{h}_{0}}^{h}\displaystyle \frac{m(h^{\prime} )g(h^{\prime} ){\rm{d}}h^{\prime} }{{kT}(h^{\prime} )}=-\mathrm{ln}\left[\displaystyle \frac{p(h)}{p({h}_{0})}\right],\end{eqnarray} \tag{ 10 }$$

where h is the geopotential height, dU = g(h)dh, h₀ is a reference height at the lower boundary, g(h) is the gravity (perpendicular to equipotential surfaces), m(h) is the mean molecular weight, k is the Boltzmann constant, and T(h) is the temperature. We assume that the temperature is constant throughout the lower and middle atmosphere and equal to the effective temperature of the planet. In order to calculate the effective temperature, we use a generic Bond albedo A_B = 0.3 and assume uniform redistribution of energy around the planet.

We assume that the mean molecular weight of the atmosphere depends only on the relative abundances of H₂, He, and H. For a solar-metallicity gas, the volume mixing ratios of these species can be approximated as (Visscher et al. 2006):

$$\begin{eqnarray}\begin{array}{rcl}{q}_{{{\rm{H}}}_{2}} & \approx & \displaystyle \frac{1.9845+{10}^{u}-\sqrt{{10}^{u}(3.9690+{10}^{u})}}{2.3670}\\ u & = & -\displaystyle \frac{23672}{T}-{\mathrm{log}}_{10}p+6.2645,\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{q}_{\mathrm{He}}=\displaystyle \frac{{f}_{\mathrm{He}}(1-{q}_{{{\rm{H}}}_{2}})+2{f}_{\mathrm{He}}{q}_{{{\rm{H}}}_{2}}}{{f}_{\mathrm{He}}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{q}_{{\rm{H}}}=1-{q}_{{{\rm{H}}}_{2}}-{q}_{\mathrm{He}},\end{eqnarray} \tag{ 13 }$$

where p is pressure in bar, and f_He = 7.93 × 10⁻² is the abundance of He relative to H (Lodders 2003). In order to test the validity of these approximations, we use the NASA Chemical Equilibrium Applications (CEA) code (Gordon & McBride 1994) to calculate the composition at temperatures of 1000, 1500, 2000, and 2500 K. Figure 1 compares the mean molecular weight based on the model with that based on the above equations. The lack of heavy elements leads to a slightly lower mean molecular weight in our simple model but otherwise, the agreement is relatively good. We note that our approach here ignores photochemistry. Since photolysis leads to more efficient dissociation of H₂ and other molecules, our results underestimate the extent of the atmosphere for cooler planets with effective temperatures of ∼1000–2000 K. At higher temperatures, the composition is expected to approach thermochemical equilibrium (Moses et al. 2011).

Zoom In Zoom Out Reset image size

2.3. Roche Lobe Overflow

Roche lobe overflow occurs when the atmosphere of the planet extends to the Roche lobe (e.g., Lecavelier des Etangs et al. 2004; Li et al. 2010). Otherwise, the outer limit of the atmosphere is the exobase where the mean free path between collisions exceeds the pressure scale height (p ≈ 10⁻¹² bar). For our purposes, Roche lobe overflow occurs when the L1 point is located below the exobase. If the temperature and mean molecular weight do not change significantly with latitude and longitude, the Roche lobe is a surface of constant pressure and density, and once the L1 point is below the exobase, the entire Roche lobe is below the exobase.

Roche lobe overflow can support substantial mass loss even in the absence of stellar XUV heating (e.g., Li et al. 2010). In general, we can estimate the mass-loss rate based on the underlying atmospheric structure by adapting Jeans escape for "geometrical blow-off" (Lecavelier des Etangs et al. 2004):

$$\begin{eqnarray}&&\dot{M}\approx 4\pi {x}_{c}{z}_{c}\displaystyle \frac{{\rho }_{c}}{2\sqrt{\pi }}\sqrt{\displaystyle \frac{2{{kT}}_{c}}{{m}_{c}}}\left({{\rm{\Gamma }}}_{c}+1\right)\exp (-{{\rm{\Gamma }}}_{c}),\end{eqnarray} \tag{ 14 }$$

where subscript c refers to the exobase, ρ_c is mass density, x_c and z_c are the substellar and polar radii, respectively, of the corresponding equipotential surface, and Γ_c is the thermal escape parameter:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{c}=\displaystyle \frac{{m}_{c}{\rm{\Delta }}U}{{{kT}}_{c}}.\end{eqnarray} \tag{ 15 }$$

Here, ΔU is the potential difference from the exobase to the Roche lobe. This potential difference approaches zero when the Roche lobe approaches the exobase. When the Roche lobe is at the exobase or below, the mass-loss rate is

$$\begin{eqnarray}&&{\dot{M}}_{R}\approx 4\pi {x}_{R}{z}_{R}\displaystyle \frac{{\rho }_{R}}{2\sqrt{\pi }}\sqrt{\displaystyle \frac{2{{kT}}_{R}}{{m}_{R}}}.\end{eqnarray} \tag{ 16 }$$

Effectively, this means that the atmosphere escapes through the Roche lobe close to the speed of sound. In Section 3.2, we use a detailed numerical model of hydrodynamic escape to show that this approximation holds for flow through the L1 point and argue that, in principle, Equation (16) provides a reasonably good estimate of the global mass-loss rate due to Roche lobe overflow.

If the nightside is substantially colder than the dayside, the density at the Roche lobe boundary could be significantly lower than on the dayside. In the unlikely case of no escape at all from the nightside, the mass-loss rate would be reduced by a factor of 2. A more conservative assumption is that Roche lobe overflow is limited to a nozzle around the L1 point. Lai et al. (2010) estimate the radius of this nozzle as:

$$\begin{eqnarray}&&{\rm{\Delta }}{R}_{L}\approx \displaystyle \frac{{c}_{s}}{{\rm{\Omega }}},\end{eqnarray} \tag{ 17 }$$

where c_s is the sound speed, and give the mass-loss rate as:

$$\begin{eqnarray}&&{\dot{M}}_{L}\approx \pi {\rho }_{R}{c}_{s}{\left({\rm{\Delta }}{R}_{L}\right)}^{2}.\end{eqnarray} \tag{ 18 }$$

For the planets considered in this paper, the mass-loss rates based on this approximation are lower by factors of 2–10 compared to Equation (16) once the atmosphere enters the Roche lobe overflow regime. We conclude that Equation (16) is accurate to a factor of a few even if conditions around the planet are significantly different.

2.4. Energy-limited Escape

Energy-limited escape is often assumed in the literature to calculate mass loss due to stellar XUV-driven atmospheric escape or photoevaporation ⁵ . For this reason, we also compare our results with this approach in Section 3. The standard formula for mass loss by energy-limited photoevaporation (Watson et al. 1981) can be written as:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{XUV}}=\displaystyle \frac{\epsilon {L}_{\mathrm{XUV}}}{4{\rm{\Delta }}U}{\left(\displaystyle \frac{{r}_{e}}{a}\right)}^{2},\end{eqnarray} \tag{ 19 }$$

where is the mass loss (or heating) efficiency, L_XUV is the XUV luminosity of the star, r_e is the effective radius in the planet's atmosphere, where XUV radiation is absorbed and ΔU is the potential difference that escaping particles must overcome. This formula assumes uniform redistribution of energy around the planet. Typically, photoionization of H and He are assumed to heat the upper atmosphere, and L_XUV includes radiation at energies higher than the ionization threshold of H at 13.6 eV (λ ≤ 91.1 nm). The apparent simplicity of the formula, particularly the linear dependency on the stellar XUV luminosity, is deceptive. Much of the related atmospheric physics is buried in the efficiency that depends on stellar luminosity and atmospheric composition (e.g., Murray-Clay et al. 2009; Koskinen et al. 2014). The lack of knowledge about atmospheric structure also makes it difficult to properly define the potential difference ΔU.

The simplest form of the equation, often adapted in studies of exoplanet evolution, assumes an isolated object for which ΔU = GM_p /R_p , where R_p is the radius of the planet that would be probed by visual transit observations. In this case, we have (e.g., Owen & Wu 2017):

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{XUV}}=\displaystyle \frac{\epsilon \pi {R}_{p}^{3}{L}_{\mathrm{XUV}}}{4\pi {a}^{2}{{GM}}_{p}}.\end{eqnarray} \tag{ 20 }$$

Here, energy-limited escape is inversely proportional to the bulk density of the planet. This explains why it returns similar mass-loss rates for a wide range of different planets, ranging from Neptune-mass planets to Jovian planets, that tend to have similar bulk densities. Given that bulk densities vary by less than an order of magnitude in general, energy-limited mass-loss rates for rocky planets are also similar in magnitude to the giant planets.

Another commonly used form of energy-limited escape includes a correction factor to approximate Roche lobe overflow. This was proposed by Erkaev et al. (2007) who wrote the potential difference in Equation (19) as

$$\begin{eqnarray}&&{\rm{\Delta }}U\approx \displaystyle \frac{{{GM}}_{p}}{{R}_{p}}\left(1-\displaystyle \frac{3}{2\eta }+\displaystyle \frac{1}{2{\eta }^{3}}\right)=\displaystyle \frac{{{GM}}_{p}}{{R}_{p}}K(\eta ),\end{eqnarray} \tag{ 21 }$$

where η = x_R /R_p and the L1 point is at

$$\begin{eqnarray}&&{x}_{R}\approx {\left(\displaystyle \frac{\delta }{3}\right)}^{1/3}\left[1-\displaystyle \frac{1}{3}{\left(\displaystyle \frac{\delta }{3}\right)}^{1/3}\right]a,\end{eqnarray} \tag{ 22 }$$

where δ = M_p /M_s . This approach assumes that both the surface of the planet at r = R_p and the Roche lobe are surfaces of constant density and that the atmosphere escapes uniformly through the Roche lobe. In line with Equation (20), it does not account for the extent of the atmosphere.

2.5. Upper Atmosphere Escape Model

3.1. A Fundamental Stability Limit

Once the atmosphere of a planet reaches the Roche lobe, the mass-loss rate rapidly increases. Based on this, we define a stability limit for gaseous envelopes by using the simple theory from Sections 2.2 and 2.3, and relate the limit to the equilibrium shape of a planet near the surface pressure level. Here, we ignore the heating of the upper atmosphere and assume instead that the effective temperature of the planet applies throughout the atmosphere up to the exobase or the Roche lobe. This produces a lower limit on mass loss, and conveniently the results depend only on M_p , R_p , M_s , a, the stellar bolometric luminosity L_s , and the Bond albedo of the planet that we set to A_B = 0.3. We focus on Sun-like host stars, so that M_s = M_⊙ and L_s = L_⊙. In order to highlight the underlying physics that determine the stability limit, we consider a tidally locked Jupiter-like planet orbiting a Sun-like star (M_p = 1898.19 × 10²⁴ kg, R_p = 71,492 km) to represent hot Jupiters and a tidally locked Uranus-like planet orbiting a Sun-like star (M_p = 86.813 × 10²⁴ kg, R_p = 25,559 km) to represent hot Neptunes. We assume that the polar radii of the planets at 1 bar are equal to the radii of Jupiter and Uranus.

The shape of the 1 bar surface changes with decreasing orbital distance as the planets become increasingly distorted by the stellar tide. Figure 2 shows the ratio of the substellar radius to the polar radius, X/Z, for Uranus-, Saturn-, and Jupiter-like planets as a function of the orbital distance. Here, the capital letters in the X/Z ratio refer to the observable surface of the planets. We use lower case letters for the perturbation of the pressure levels generally. The results are similar for Jupiter- and Uranus-like planets because they have similar bulk densities. The ratio is close to unity at orbital distances of a ≳ 0.03 au (P ≳ 1.9 days). At orbital distances of a ≲ 0.02 au (P ≲ 1.03 days), the ratio increases rapidly with decreasing orbital distance. Saturn has a lower bulk density, and therefore the ratio begins to increase at somewhat wider orbits. The results are consistent with the classical Roche limit for tidal destruction of a fluid planet:

$$\begin{eqnarray}&&{R}_{L}\approx 2.455{R}_{s}{\left(\displaystyle \frac{{\rho }_{s}}{{\rho }_{p}}\right)}^{\tfrac{1}{3}},\end{eqnarray} \tag{ 23 }$$

where ρ_s and ρ_p are the bulk densities of the star and the planet, respectively. Hereafter, we focus on the Jupiter- and Uranus-like planets as an example. In Section 4, we generalize our approach to different planets. For both Jupiter and Uranus, the Roche limit is R_L ≈ 0.0118 au (P = 0.47 days). In Figure 2, the X/Z ratio exceeds 1.2 at a ≲ 0.0124 au (P ≲ 0.5 days), as the planets approach their Roche limits. We note that X/Z ≈ 1.5 once the surface of the planet fills the Roche lobe.

Zoom In Zoom Out Reset image size

As tidal distortion increases, the upper atmosphere reaches the Roche lobe first. To illustrate this, we use the approach from Section 2.2 to calculate the extent of the atmosphere for the planets. We remind the reader that we ignore the heating of the upper atmosphere by stellar XUV radiation and photochemistry here. Stellar XUV heating would produce a hotter upper atmosphere and photochemistry would lead to more efficient ionization and dissociation of molecules, both of which increase the extent of the atmosphere and mass-loss rate. Thus, the results here represent lower limits on both.

We integrate Equation (10) numerically at the substellar point by using a pressure grid that ranges from p₀ = 1 bar to the exobase at p = 10⁻¹² bar with a pressure level spacing of 0.01 scale heights. The top of the atmosphere is at p = 10⁻¹² bar if the L1 point is above this level. If not, the top of the atmosphere is the L1 point. Figure 3 shows the radius of the L1 point for Uranus- and Jupiter-like planets orbiting a Sun-like star. It also shows the radius at the L1 point for a Uranus-like planet based on the approximate expression from Erkaev et al. (2007) (Equation (22)). Figure 4 shows the corresponding pressure at the top of the atmosphere. By our definition, once the Roche lobe reaches down to the exobase, the atmosphere undergoes Roche lobe overflow and the pressure at the top of the atmosphere (L1 point) exceeds 10⁻¹² bar (see Section 2.3). A Uranus-like planet has lower gravity and a more extended atmosphere than a Jupiter-like planet. Thus, a Uranus-like planet undergoes Roche lobe overflow at a ≲ 0.033 au (P ≲ 2.2 days), while a Jupiter-like planet undergoes Roche lobe overflow at a ≲ 0.014 au (P ≲ 0.6 day).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Many studies of exoplanet evolution rely on simple energy-limited escape to represent stellar XUV-driven photoevaporation (e.g., Owen et al. 2020). Previous work has shown that energy-limited escape can substantially overestimate mass loss unless the atmosphere undergoes hydrodynamic escape and adiabatic cooling due to escape mostly balances heating (Koskinen et al. 2014). Here, we focus on the opposite end of escape regimes. Once the middle atmosphere reaches the Roche lobe, the atmosphere effectively enters a "low gravity" regime and is free to escape (see also Jackson et al. 2017; Kubyshkina et al. 2018). This regime cannot be captured by adding a simple correction to energy-limited escape that represents the potential difference between the surface of the planet and the L1 point (Erkaev et al. 2007; Equation (21) here). Even with this correction, energy-limited escape substantially underestimates the mass-loss rate due to Roche lobe overflow.

We calculate the corresponding mass-loss rates due to Roche lobe overflow by using Equations (14)–(16). As we note above, these are lower limits on mass loss since we ignore the heating of the upper atmosphere entirely. Figure 5 shows the mass-loss rate at different orbital distances in units of M_p Gyr⁻¹ while Figure 6 compares it with the mass-loss rate based on energy-limited photoevaporation (Equations (19) and (21) with r_e = R_p ). For energy-limited escape, we use the mean solar XUV flux of 1.6 W m⁻² at 0.05 au integrated over the wavelength range of 0.1–91.1 nm (Ribas et al. 2005), with an efficiency factor of = 1. As we note earlier, energy-limited escape rates for Uranus- and Jupiter-like planets are similar because the bulk densities of the planets are similar.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As indicated by Figures 4 and 6, stellar XUV heating of the upper atmosphere is the dominant driver of atmospheric escape at orbital distances where the Roche lobe is above the exobase. At short orbital distances where the Roche lobe is below the exobase, however, Roche lobe overflow dominates over photoevaporation. Under these circumstances, energy-limited escape significantly underestimates the mass-loss rate because it does not account for the extent of the atmosphere. In this case, escape is effectively powered by the bolometric luminosity of the star, which is more than 5 orders of magnitude higher than the XUV luminosity. Tidal heating of the planet can provide an additional energy source (Li et al. 2010) but is not necessarily required.

For a Uranus-like planet orbiting a Sun-like star, the cumulative mass lost in 1 Gyr exceeds the mass of the planet at a ≲ 0.024 au (P ≲ 1.36 days). At this limiting distance, the pressure at the Roche lobe is p_R = 3.3 × 10⁻⁹ bar and X/Z = 1.017. The same limits for complete mass loss in only 100 Myr are a ≲ 0.022 au (P ≲ 1.2 days), p_R = 4 × 10⁻⁸ bar, and X/Z = 1.024. For a Jupiter-like planet, the limits for complete mass loss in 1 Gyr are a ≲ 0.0124 au (P ≲ 0.5 day), p_R = 4.6 × 10⁻⁸ bar, and X/Z = 1.19. The same limits for complete mass loss in only 100 Myr are a ≲ 0.0122 au (P ≲ 0.5 day), p_R = 4.8 × 10⁻⁷ bar bar, and X/Z = 1.21. In other words, Jovian planets experience tidal destruction once the surface X/Z ratio exceeds 1.2 near the classical Roche limit. The limiting X/Z ratio for sub-Jovian planets, however, is smaller, and they experience envelope loss outside of the Roche limit. Thus, the visible surface of a planet need not fill the Roche lobe for the planet to lose its envelope over a relatively short time. Our results show that rapid mass loss occurs once the 10⁻⁸–10⁻⁷ bar level in the upper atmosphere reaches the Roche lobe. We note that the results here do not account for the change in mass loss once the envelope is substantially depleted and instead assume that the mass-loss rate is constant over time (see Section 4 for related discussion).

3.2. Numerical Modeling

In order to test the validity of the simplifying assumptions above, we use a multispecies model of hydrodynamic escape (see Section 2.5 and Appendix B) to simulate the transition from stellar XUV-driven escape to Roche lobe overflow in more detail. Obviously, this model accounts for the heating of the upper atmosphere by stellar XUV radiation, unlike the results presented above. We focus on a Uranus-like planet orbiting a Sun-like star. We calibrate the model at 0.05 au where Roche lobe overflow is not expected to be significant. This provides a direct comparison to models of the prototypical hot Jupiter HD209458b that has a similar orbital distance around a Sun-like star (Koskinen et al. 2013a, 2013b). Our reference Model A does not include Roche lobe overflow. In Model B, we replace spherically symmetric gravity with the substellar gravity based on the Roche potential (Equation (9)). This corresponds to escape through the L1 point. In both cases, we use the same mean solar XUV spectrum as Koskinen et al. (2013a, 2013b) and assume uniform redistribution of energy around the planet, in line with energy-limited escape.

The radial mass flux constant predicted by Model A at a = 0.05 au is F_M = ρ vr² = 1.5 × 10⁶ kg s⁻¹ sr⁻¹. Multiplying this by a solid angle factor 4π gives a global mass-loss rate of $\dot{M}=$ 1.9 × 10⁷ kg s⁻¹. Owen et al. (2020) compared 3D simulations of the upper atmosphere of the warm Neptune GJ436b (Shaikhislamov et al. 2018) with 1D globally averaged simulations of the same planet by our model (Loyd et al. 2017a) and found good agreement on the mass-loss rates. This means that our 1D model can be used to calculate global mass-loss rates and that atmospheric extent and the details of the energy balance in the upper atmosphere are likely more important for predicting accurate mass-loss rates than adding multiple dimensions to the model (see Section 4 for related discussion).

Figure 7 shows the temperature and outflow velocity profiles for Model A at a = 0.05 au. The lower boundary of the model is at r₀ = 1.34R_p (p₀ = 10⁻⁶ bar, T₀ = 1, 140 K) where the boundary conditions are consistent with our models of the lower and middle atmosphere (see Section 2.2). Above the lower boundary, the temperature first slightly decreases with the radius, and then increases rapidly to about 4400 K around r ≈ 3 R_p . At higher radii, the temperature increases slowly with the radius to about 5540 K at the upper boundary of the model. The outflow velocity in the model is significant at r ≳ 2 R_p , and it increases slowly with the radius to about 9.8 km s⁻¹ at the upper boundary. The adiabatic sound speed ranges from about 2.5 km s⁻¹ at the lower boundary to about 8 km s⁻¹ at the upper boundary. The corresponding range of values for the isothermal sound speed is from about 2–6 km s⁻¹. The adiabatic and isothermal sonic points are at 6.7 and 4.7 R_p , respectively, while the L1 point is at 7.1 R_p . The exobase is above the upper boundary of the model, and therefore the atmosphere escapes hydrodynamically.

Zoom In Zoom Out Reset image size

The composition is shown by Figure 8. The volume mixing ratio of H₂ at the lower boundary is about q₀ = 0.84 while H is a minor species with q₀ = 0.026. H becomes the dominant species at around r ≈ 2 R_p , but the density of H₂ remains significant at all altitudes. The dominant ion is H⁺, with a density that is at least 2–3 orders of magnitude higher than the densities of other ions. Charge balance means that the electron density is equal to the proton density. The density of neutral H is higher than the proton density at all radii.

Zoom In Zoom Out Reset image size

Simple ionization ${H}_{2}+h\nu \to {H}_{2}^{+}+e$ is the dominant photoionization channel for H₂, and the resulting ${{\rm{H}}}_{2}^{+}$ reacts quickly with H₂ to form ${{\rm{H}}}_{3}^{+}$ (see Appendix B). Dissociative recombination of ${{\rm{H}}}_{3}^{+}$ is a rapid process that keeps the density of ${{\rm{H}}}_{3}^{+}$ relatively small. Despite its small abundance, infrared emissions by ${{\rm{H}}}_{3}^{+}$ are nevertheless significant. These emissions are regularly detected on solar system giant planets where they are used to probe the ionosphere with remote observations. On Jupiter, they also act as an important coolant (Miller et al. 2013). The density of the other simple molecular ion, HeH⁺, created by the reaction He⁺ + H₂ → HeH⁺ + H, is also significant in our model. Rovibrational emissions by HeH⁺ have been detected recently from planetary nebulae, constituting the first detection of the ion in astrophysical environments (Güsten et al. 2019; Neufeld et al. 2020).

The stellar XUV heating peak is at a radius of 1.66 R_p (p = 2.1 × 10⁻⁹ bar), as demonstrated by Figure 9, which shows the prominent heating and cooling terms. The peak volume heating rate is 1.8 × 10⁻⁸ W m⁻³. At r ≲ 2.3 R_p , heating is mostly balanced by ${{\rm{H}}}_{3}^{+}$ cooling. At higher radii, the outflow velocity is significant, and adiabatic cooling due to the expansion of the atmosphere primarily balances the stellar heating rate, with a smaller contribution from ${{\rm{H}}}_{3}^{+}$ cooling. Advection makes a small contribution to cooling and heat conduction is largely negligible. We note that the rate for heat conduction in Figure 9 is diagnostic only. In practice, the model uses a semi-implicit method to update the temperature profile for the effect of conduction, and the diagnostic rate is calculated ex post facto with finite difference methods that produce numerical noise at small values.

Zoom In Zoom Out Reset image size

We now compare our results with energy-limited escape to find values of and ΔU in Equation (19) that agree with our detailed simulations. With a stellar XUV luminosity of L_XUV = 1.12 × 10²¹ W, the flux incident on Model A at a = 0.05 au is F_XUV = 1.6 W m⁻². Thus, the globally averaged flux incident on the planet is about 0.4 W m⁻². The effective flux that heats the upper atmosphere can be obtained from our model as

$$\begin{eqnarray}&&\displaystyle \frac{\epsilon {F}_{\mathrm{XUV}}}{4}={\int }_{{r}_{0}}^{\infty }\left[H(r)-C(r)\right]{\rm{d}}r=0.11\ {\rm{W}}\,{{\rm{m}}}^{-2},\end{eqnarray} \tag{ 24 }$$

where H(r) and C(r) are the radiative volume heating and cooling rates, respectively, from the model. Therefore, the model suggests a mass loss efficiency = 0.28 for energy-limited escape. With this efficiency, Equation (20) gives a mass-loss rate of 4 × 10⁶ kg s⁻¹, which is 4.8 times lower than the actual mass-loss rate predicted by Model A. The model shows, however, that the heating peak is at r_e = 1.66R_p . Replacing R_p in Equation (20) with the radius r_e at the heating peak yields an energy-limited escape rate of 1.8 × 10⁷ kg s⁻¹ that agrees roughly with Model A.

In this case, neglecting the extent of the atmosphere underestimates the energy-limited mass-loss rate by a factor of 5. This factor is larger than similar factors derived for typical hot Jupiters because of the larger extent of the envelopes of Neptune-mass planets. Owen & Wu (2017) suggested that the effect of the "expansion radius" can be folded into the efficiency factor. This is not an ideal solution for low-mass planets with H/He envelopes because our results imply efficiencies close to unity or even higher that are unlikely to be adopted by modelers. In fact, the escape rate of 1.8 × 10⁷ kg s⁻¹ is comparable to the energy-limited escape rate in Figure 6 that we obtained by using Equations (19) and (21) with 100% efficiency and r_e = R_p .

Alternatively, the base radius for escape can be estimated from hydrostatic equilibrium as (e.g., Lopez 2017)

$$\begin{eqnarray}&&{R}_{b}\approx {R}_{p}+{H}_{1}\mathrm{ln}\left(\displaystyle \frac{{p}_{1}}{{p}_{b}}\right),\end{eqnarray} \tag{ 25 }$$

where H₁ is the scale height based on m₁ = 2.3 amu and the effective temperature T₁ of the planet, p₁ = 1 bar is the surface pressure, and p_b = 1 nbar is the base pressure for escape. At a = 0.05 au, we obtain H = 470 km at the 1 bar level and R_b ≈ 1.38 R_p at the 1 nbar level. A better approximation can be obtained by including the dependency of gravity on radius in the above equation:

$$\begin{eqnarray}&&{\rm{\Gamma }}({R}_{b})\approx {{\rm{\Gamma }}}_{1}+\mathrm{ln}\left(\displaystyle \frac{{p}_{b}}{{p}_{1}}\right),\end{eqnarray} \tag{ 26 }$$

where

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{{GM}}_{p}{m}_{1}}{{{kT}}_{1}r}\end{eqnarray} \tag{ 27 }$$

is the thermal escape parameter. We obtain Γ₀ = 54.6 and Γ = 33.9, where the latter corresponds to R_b ≈ 1.6 R_p . This is a straightforward way of roughly reproducing the mass-loss rate predicted by our detailed model with = 0.3. It is reasonably accurate even though it does not account for all the complex physics included in our model.

Model B represents escape through the L1 point. The radial mass flux based on the model is F_M = 1.7 × 10⁶ kg s⁻¹ sr⁻¹. This is higher by a factor of 1.13 than the radial mass flux predicted by Model A. The outflow velocity in Model B is somewhat higher than in Model A, reaching a velocity of 12 km s⁻¹ at the upper boundary. Due to the increased mass-loss rate and the strength of adiabatic cooling, the temperature of 4740 K at the upper boundary is lower than in Model A. Otherwise, Models A and B are qualitatively similar. The adiabatic and isothermal sonic points in Model B are at 5.2 R_p and 4.2 R_p , respectively. The sonic points are below the L1 point of 7.1 R_p , and the equipotential surfaces that coincide with them are approximately spherical. Therefore, the atmosphere escapes through the sonic surface at a global mass-loss rate of 4π F_M = 2.1 × 10⁷ kg s⁻¹.

3.2.1. Comparison with Hot Jupiter Simulations

In order to highlight the main qualitative differences between escape models for hot Jupiters and our model for a Uranus-like planet, we briefly describe the key differences between the model above and the simulations of the prototype hot Jupiter HD209458b that also orbits a Sun-like star at a distance of about 0.05 au. We predict temperatures of 4000–5000 K for the hot Uranus that are significantly cooler than the peak temperature of 8000–12000 K predicted for HD209458b (Koskinen et al. 2013a, 2013b). In addition, the temperature profile for hot Uranus lacks a well-defined peak that appears in hot Jupiter simulations. Both of these differences are explained by the lower gravity of the hot Uranus model that makes it more susceptible to the effects of atmospheric escape that acts as the primary cooling mechanism above the heating peak.

Models of HD209458b predict significant dissociation of H₂ near the bottom of the thermosphere that removes practically all H₂ molecules and related molecular ions from the outflow (Yelle 2004; Koskinen et al. 2013a). In contrast to this, the density of H₂ is significant at all radii of the hot Uranus model. This leads to the appearance of the molecular ions ${{\rm{H}}}_{3}^{+}$ and HeH⁺ at high altitudes, neither of which are significantly produced in hot Jupiter simulations. Another key difference between the composition of the hot Uranus model and the hot Jupiter models is that the density of neutral H is higher than the proton density at all radii. The prevalence of neutral H and H₂ at high radii in the hot Uranus model is due to a combination of a relatively low temperature and high escape flux that replenishes H₂ and H at high altitudes where they would otherwise be predominantly dissociated and ionized (see Appendix B).

3.2.2. Transition to Roche Lobe Overflow

Here, we present results from Model B at orbital distances of a = 0.06, 0.04, 0.03, and 0.02 au. Figure 6 compares the global mass-loss rates predicted by the model for each orbital distance with the minimum mass-loss rates based on the simple treatments of Roche lobe overflow (see Section 3.1) and energy-limited photoevaporation obtained for = 1 with the Roche potential correction factor. For each model, the lower boundary conditions are consistent with the structure and composition of the lower and middle atmosphere (see Section 2.2). We note that, as before, Model B applies at the substellar point. For each model, the tidal perturbation of the lower and middle atmosphere is properly reflected by the lower boundary radius that is consistent with the x/z ratio based on the Roche potential at p₀ = 10⁻⁶ bar.

At orbital distances of a = 0.04–0.06 au, the mass-loss rate is close to the energy-limited escape rate. As expected, the dependency of the mass-loss rate on the orbital distance is close to ∼1 a⁻² at these orbital distances. A slight deviation from this trend is caused by the effect of the stellar tide that enhances the mass-loss rate as the orbital distance decreases. At a = 0.03 au, the mass-loss rate is higher by a factor of 2.6 than at a = 0.04 au. This factor is larger than the expected factor of 1.8 increase based on a 1 a⁻² dependency, because Roche lobe overflow dominates over photoevaporation at a ≲ 0.03 au. At a = 0.02 au, the global mass-loss rate predicted by Model B is 2.8 × 10¹¹ kg s⁻¹. This value is twice as large as the minimum mass-loss rate of 1.4 × 10¹¹ kg s⁻¹ due to Roche lobe overflow that we obtain in Section 3.1, while it is about 3 orders of magnitude higher than the energy-limited escape rate. We note that 2D simulations with a uniform redistribution of energy by Guo (2013) indicate that globally averaged 1D hydrodynamic models with substellar gravity may overestimate Roche lobe overflow by a factor of 2. If this is true, our numerical model predicts a mass-loss rate that agrees well with the simple theory that underlies Equation (16). This shows that (i) Equation (16) and the related theory in Section 2.3 provide a reasonably good approximation of mass loss due to Roche lobe overflow and (ii) photoevaporation does not significantly add to the mass-loss rate once Roche lobe overflow sets in.

Figure 10 shows the temperature and outflow profiles for Model B at a = 0.02 au. The lower boundary radius at the substellar point is r₀ = 2.0 R_p (p₀ = 10⁻⁶ bar). Due to strong adiabatic cooling and advection based on the large escape rate, the temperature decreases with the radius above the lower boundary to T = 650 K at the L1 point, which is located at x_R = 2.8 R_p . The adiabatic sonic point is at r ≈ 2.7 R_p , and therefore our numerical model supports the common approximation that the atmosphere escapes through the L1 point at the speed of sound (e.g., Li et al. 2010; Lai et al. 2010).

Zoom In Zoom Out Reset image size

3.3. The Shape of Known Transiting Planets

We now explore the population of known transiting planets to make sure that their implied shape complies with the new limit on the X/Z ratio and identify the most extreme planets at risk of significant mass loss. Here, we focus only on transiting planets for which both the mass and radius have been measured. We include planets with orbital periods of less than 10 days (a ≈ 0.09 au around a Sun-like host star) and masses higher than 10 M_⊕. For simplicity, we assume that the planets are tidally locked and have circular orbits. We obtain planet properties from the Extrasolar Planets Encyclopedia (www.exoplanet.eu). We screen for confirmed planets with a primary transit detection, measured planet mass, and known stellar mass and luminosity.

We solve Equation (8) numerically to calculate the implied ratio of the substellar radius to the polar radius X/Z at the visible surface of the planets (i.e., the level probed by visual broadband transit). In all cases, we use Kepler's third law to calculate the orbital period and Ω based on the observed M_p , M_s , and a instead of using the published orbital periods (see Section 3.4). The results are shown in Figure 11. Planets WASP-12b, WASP-19b, and WASP-121b have the highest X/Z ratios of 1.19, 1.17, and 1.15, respectively. Two of these planets, WASP-12b and WASP-121b, show evidence of the escape of multiple metals from their atmospheres, indicating large mass-loss rates (Fossati et al. 2010; Sing et al. 2019). The mass-loss rates of these planets exceed the stellar XUV-driven mass-loss rate.

Zoom In Zoom Out Reset image size

As expected, there are no transiting planets in this sample with implied X/Z ≳ 1.2, corresponding to a surface distortion that would lead to overflow of the atmosphere and rapid mass loss for Jovian planets. Furthermore, there are no transiting planets with masses between 10 M_⊕ and 20 M_⊕, similar to Uranus and Neptune, with X/Z ratios higher than 1.06. The highest ratio of 1.06 in this mass range applies to planet K2-266, b while the rest of the hot Neptunes have X/Z ratios lower than 1.01. The ratio for K2-266 b is relatively high, but it is also not particularly reliable. The uncertainty on the planet's mass is relatively large, and it could be twice as large as the listed mean value of 11.3 M_⊕ (Rodriguez et al. 2018). K2-266 b still has the highest X/Z ratio if we increase the upper mass limit of the planets in the sample from 20 to 100 M_⊕ to include Saturn-like planets. This supports the results derived from the simple theory in Section 3.1. Surviving Neptune- and Saturn-like planets have X/Z ratios close to unity while all planets with higher X/Z ratios are hot Jupiters.

3.4. A Note on the Planet Parameters

As we note before, many studies of exoplanet evolution rely on energy-limited escape to represent photoevaporation. This approach underestimates mass-loss rates due to Roche lobe overflow. Therefore, we explore the dependency of our results on an extended set of planet properties and relate them to the population of planets detected by Kepler. Figure 12 shows the periods and radii of Kepler planets orbiting main-sequence stars, following a revision of their properties based on stellar properties from the Gaia Data Release 2 and DR25 Kepler Stellar Properties Catalogue (Berger et al. 2020). The dashed line shows the limit where the minimum mass-loss rate due to Roche lobe overflow exceeds energy-limited photoevaporation. In order to calculate this limit, we create a grid of planet masses from 10 M_⊕ to 960 M_⊕ with a spacing of 5 M_⊕ and use the empirical mass–radius (M–R) relationship of Chen & Kipping (2017) to obtain radii for the planets. This relationship is mostly based on transiting planets that typically orbit close to their host stars. For a clear atmosphere, the visual transit probes the 0.01–0.1 bar level in the atmosphere, and we assume that the radius applies at the 0.1 bar level.

Zoom In Zoom Out Reset image size

The Kepler population includes ultrashort period, or USP, planets with periods of P ≤ 1 day. These planets generally have radii of less than about 2 R_⊕, and they are unlikely to have significant gaseous envelopes. Sanchis-Ojeda et al. (2014) suggested that larger USP planets lose their gaseous envelopes to photoevaporation. Our results confirm that larger USP planets can lose their envelopes to rapid Roche lobe overflow. The cutoff for catastrophic mass loss is at planet radii of 13–14 R_⊕, which correspond to planet masses comparable to or slightly lower than the mass of Jupiter. The population of smaller USP planets can therefore include remnant cores of both subgiant and giant planets. Curiously, there are two USP candidate planets with radii larger than 2 R_⊕. These planets orbit stars with effective temperatures of 5744 and 5889 K and radii of 1.34 and 1.25 R_⊙, respectively. Sanchis-Ojeda et al. (2014) argued that such planets may retain some of their gaseous envelopes if they have large enough core masses and densities. This would make them significant outliers in the M–R distribution, however, and thus the nature of these mysterious planets remains unclear.

Berger et al. (2018) refine the limits of a hot super-Earth and Neptune "desert", a dearth of planets with R_p = 2–4 R_⊕ and host star bolometric fluxes of F ≳ 650 F_⊙. The flux limit corresponds to an orbital distance of a ≈ 0.04 au (P ≈ 3 days ) for a planet orbiting a Sun-like star. Based on the M–R relationship, the lowest mass planet in our grid has a radius of R_p ≈ 3 R_⊕. Extrapolating to slightly lower masses indicates that Roche lobe overflow can play a significant role in shaping the hot Neptune desert. This desert is a part of a larger sub-Jovian desert that extends to P ≈ 4 days, including planets with radii of R_p ≈ 4–10 R_⊕ (e.g., Mazeh et al. 2016). Our results indicate that a part of this desert can also be cleared by Roche lobe overflow. We note that these results are likely to be conservative because gaseous planets have higher internal energies, larger radii, and thus lower bulk densities when they form, making them more sensitive to mass loss in their youth (e.g., Kurokawa & Nakamoto 2014).

Beyond the limit for significant Roche lobe overflow, escape is powered mostly by the stellar XUV radiation, and mass loss is roughly energy-limited as long as escape is hydrodynamic and the extent of the atmosphere is correctly represented. We cannot speculate on the consequences of escape on the planets in this regime as the results depend significantly on the assumed interior structure and envelope composition (e.g., Owen & Wu 2017), the modeling of which is beyond the scope of this work. For the same reason, we have not explored the consequences of Roche lobe overflow with detailed models of planetary or orbital evolution in this work.

Formally, energy-limited escape can be defined as mass loss that is directly proportional to the net energy deposited as heat in the upper atmosphere. We note that even when escape is energy-limited, the heating efficiency of the upper atmosphere remains a significant uncertainty in the models. Not only does the heating efficiency depend on the detailed composition of the upper atmosphere, but it also depends on the incident stellar XUV flux (e.g., Murray-Clay et al. 2009). As if this was not enough, the heating efficiency is also likely to depend on the spectral energy distribution of the XUV flux. For example, it is well known that young Sun-like stars are much more active than the current Sun and emit more of their XUV flux in X-rays that penetrate deeper to the lower thermosphere and middle atmosphere. Thus the heating and cooling profiles in the atmospheres of young planets can differ significantly from planets that orbit current Sun-like stars, and the heating efficiency is likely to change significantly over time.

The key limitations of our models are that they are 1D and only apply to hydrogen and helium envelopes. They also do not include the possibility of a significant planetary magnetic field or the interaction of the planetary upper atmosphere and exosphere with the stellar wind. The degree to which the differences in temperature and composition between the day and nightsides and subsequent horizontal dynamics affect the mass-loss rate is not immediately obvious. One-dimensional models are typically able to include more detailed physics and remain tractable, while multidimensional models are often based on simplified radiative transfer and chemistry to provide a more realistic representation of the redistribution of energy and different species around the planet. This poses a key problem for objective comparisons between 1D and multidimensional models. It is often difficult to work out if the differences between models are due to horizontal dynamics or some other assumptions about chemistry, heating and cooling, lower atmosphere structure, or escape mechanism. In this context, recent work by Guo (2013) and Owen et al. (2020) who compared different models with similar physics provide useful insights.

Owen et al. (2020) compared predictions by our model for the warm Neptune GJ436b (Loyd et al. 2017b) with similar but 3D simulations of Shaikhislamov et al. (2018). There is a relatively good agreement on globally averaged atmospheric structure between the models, and most importantly the mass-loss rates agree to within 10%. We suspect that this is at least in part due to the fact that in extended upper atmospheres, stellar XUV radiation passes the terminator to the nightside, leaving only a small region behind the planet in complete shadow (e.g., see the lower right panel of Figure 2 in Koskinen et al. 2007). Guo (2013), on the other hand, compared results from 1D and 2D versions of his model to explore Roche lobe effects and temperature differences between the day and nightside on a planet like HD209458b. In this case, the symmetry axis of the 2D model was aligned along the star–planet line from the substellar point to the antistellar point. The simulations show that a globally averaged 1D model with substellar gravity overestimates the mass-loss rate by a factor of 2 compared to a 2D model with uniform heating around the planet. If one also includes the change in temperature between the day and nightsides in the 2D model, the mass-loss rate is about 7 times smaller than the mass-loss rate predicted by the globally averaged 1D model with substellar gravity. Note, however, that these differences could be overestimated because dynamics in the 2D model is still restricted, and the difference between the 1D and 3D models compared by Owen et al. (2020) is smaller.

While many models with magnetic fields and stellar wind interaction have been developed for transiting planets (e.g., Ip et al. 2004; Erkaev et al. 2005; Lipatov et al. 2005; Preusse et al. 2007; Trammell et al. 2011; Owen & Adams 2014; Strugarek et al. 2014; Matsakos et al. 2015), few if any models of exoplanet evolution include them. This is understandable, given that observational constraints on the nature and evolution of stellar winds and planetary magnetic fields outside of the solar system are weak or nonexistent. Their effect on the planetary upper atmosphere and escape rates is also not fully understood. On one hand, a magnetic field may impede the escape of ions at low and middle latitudes and restrict it to polar winds (e.g., Trammell et al. 2011). If the ion–neutral collision frequency is sufficiently high, trapped ions may also impede the escape of the neutral atmosphere that is not otherwise impeded by the magnetic field. On the other hand, a planetary magnetic field, particularly when combined with the incident stellar wind, can lead to strong resistive (Joule) heating in the upper atmosphere (Cohen et al. 2014), not only in extended auroral zones but also in nonauroral regions where electric currents can be generated by a wind dynamo (Koskinen et al. 2014). This additional heating could substantially increase the escape rate, even if escape only takes place at middle and high latitudes.

Koskinen et al. (2013b) used arguments based on plasma β (the ratio of the thermal pressure to magnetic pressure) and the second Cowling number Co (the ratio of magnetic pressure to the planetary wind dynamic pressure) to estimate the maximum magnetic field strength that would allow for the unimpeded escape of ions from a planet like HD209458b. They inferred a limiting magnetic moment of about 0.04 m_J , where m_J = 1.5 × 10²⁷ A m² is the magnetic moment of Jupiter, corresponding to an equatorial magnetic field strength of about 0.1 G at the 1 R_p level. Similar limiting values for hot Jupiters have been inferred in subsequent work (e.g., Owen & Adams 2014). The maximum allowed magnetic field strength, however, increases significantly when Roche lobe overflow sets in. For example, the thermal pressure at the L1 point in our Model B of a hot Uranus at 0.02 au is 40 nbar while the dynamical pressure is 58 nbar. For the magnetic pressure to exceed these values, the required equatorial magnetic field strength at the surface pressure level is about 27 G (corresponding to about 1.2 G at the L1 point). This field is much stronger than the 1 bar equatorial magnetic field of 4.3 G on Jupiter, which is the strongest planetary magnetic field in the solar system.

The effects of different assumptions about composition and radiative transfer on the mass-loss rate could be larger than the effect of multidimensional dynamics or magnetic fields that appear to be limited to a factor of a few or an order of magnitude at most. Generally, the extent of the atmosphere and the mass-loss rate depend on the temperature and mean molecular weight. Our simplified approach reproduces the mean molecular weight of a solar composition atmosphere in chemical equilibrium. It ignores photochemistry in the middle atmosphere than can dissociate molecules more efficiently and produce a more extended atmosphere. The potential of photochemistry to significantly alter our results, however, is debatable, particularly for hotter atmospheres that are expected to be close to chemical equilibrium. The overall effect of the composition, particularly the presence of heavier atoms and molecules, on the heating and cooling rates, however, could alter the temperature profile in the middle and upper atmosphere.

The effect of the composition is particularly important if the envelope metallicity is higher than solar. This is the case for the ice giants Uranus and Neptune in the solar system. Several studies also argue that supersolar metallicity is required to explain observations of transiting super-Earths and warm Neptunes (e.g., Lavvas et al. 2019). Further support for this idea may be inferred from the empirical mass–radius relationship for exoplanets that is largely based on properties of transiting close-in planets (Chen & Kipping 2017). Interior structure models predict that a pure hydrogen and helium envelope of just a few percent of the planet mass around an Earth-like rock-metal core roughly doubles the radius of a super-Earth (e.g., Owen & Wu 2017). The mean slope of the observed mass–radius relationship above the rocky planet limit does not appear to be as steep as this would imply, suggesting that super-Earths typically do not have pure hydrogen and helium or solar-metallicity envelopes.

The potential of an envelope with significantly supersolar metallicity to lose mass to escape could be much lower than that of a pure hydrogen and helium envelope. This would not change the fundamental results or the theory in this work, but it would lead to different quantitative outcomes for mass-loss and stability limits. For example, from solar to 1000 times solar metallicity, the mean molecular weight of the atmosphere increases from 2.3 amu to 15–17 amu (Lavvas et al. 2019). This leads to smaller scale height and extent of the atmosphere. At the same time, temperatures in the middle and upper atmosphere could become cooler. Enhanced metallicity produces larger abundances of radiatively active molecules such as H₂O, CH₄, and CO₂ that can cool the atmosphere and thereby reduce the mass-loss rate. For example, the thermosphere of Venus is much cooler than the thermosphere of the Earth because it is primarily composed of CO₂, which is an efficient coolant unlike the homonuclear diatomic molecules N₂ and O₂ that dominate on the Earth.

In this work, we focus on the consequences of mass loss from the closest-in exoplanets, with an emphasis on the transition from stellar XUV-driven escape of the upper atmosphere to Roche lobe overflow of the middle atmosphere that is primarily heated by the stellar bolometric luminosity. This work was inspired by the NUV observations of WASP-12b and WASP-121b (Fossati et al. 2010; Sing et al. 2019), together with earlier work on related theory (e.g., Lecavelier des Etangs et al. 2004; Li et al. 2010). The visible surface of a planet need not fill the Roche lobe for the planet to lose its envelope over a relatively short time (see also Jackson et al. 2017), and we quantify the limit for significant envelope loss in terms of atmospheric structure and exoplanet system properties. It occurs once the 10⁻⁸–10⁻⁷ bar level in the upper atmosphere reaches the Roche lobe. When this happens, stellar XUV heating and photoevaporation driven by it are no longer important because the optical depth to ionizing radiation in the XUV reaches unity at around the 10⁻⁹ bar level.

The equilibrium shape based on the Roche potential can be used to predict the stability of gaseous envelopes. We represent the equilibrium shape of close-in planets by calculating the ratio of the substellar radius to the polar radius (X/Z) at the level in the atmosphere that would be probed by visual broadband transit observations. At the Roche lobe, this ratio would be about 1.5. For a Jupiter-like planet orbiting a Sun-like star, a ratio of X/Z ∼ 1.2 implies rapid envelope loss. The same limit for a Uranus-like planet orbiting a Sun-like star is only X/Z ∼ 1.02. This is because the Uranus-like planet has a lower surface gravity and relatively more extended envelope than hot Jupiters. The latter conclusion also applies more generally to mini-Neptunes and super-Earths. These results are supported by the properties of known exoplanets. Among the known transiting planets for which we have roughly reliable masses and radii, there are none with implied X/Z ratios that exceed 1.2. In addition, all planets with masses less than 100 M_⊕ have X/Z ratios close to unity, while the few planets with higher ratios are all hot Jupiters.

As expected, planets WASP-12b and WASP-121b have some of the highest X/Z ratios in our sample. This explains the detection of substantial mass loss that includes metals and heavy elements from these planets (Fossati et al. 2010; Sing et al. 2019). In fact, the surface distortion of these planets is so significant that Roche lobe effects should be taken into account in studies of their lower and middle atmospheres and are required for the interpretation of the entire UV–IR transmission spectrum. How fast these planets might lose their envelopes, however, is subject to debate. This is because the mass-loss rate can vary substantially with atmospheric structure and model assumptions for planets that approach the stability limit. What is clear is that the mass-loss rate from these planets exceeds the stellar XUV-driven escape rate and that their upper atmospheres are fundamentally different from more "mainstream" hot Jupiters.

Many studies of exoplanet evolution rely on simple energy-limited escape to represent stellar XUV-driven atmospheric escape (e.g., Owen & Wu 2017). This approach significantly underestimates Roche lobe effects that cannot fully be captured by adding a simple correction to energy-limited escape based on the potential difference between the surface pressure level of the planet and the L1 point (Erkaev et al. 2007). Studies of exoplanet evolution using this approach underestimate the effect of Roche lobe overflow that is likely to play a significant role in clearing out the super-Earth and hot Neptune deserts at the shortest orbits. Our results indicate that Roche lobe overflow can effectively clear out the sub-Jovian desert at orbital periods of less than about two days, even if we do not account for the early evolution of the planets when mass-loss rates are likely to be even higher. We look forward to future studies of exoplanet populations that incorporate these results in more complete, coupled models of mass loss and evolution.

We use a simple approach to represent the extent of the lower and middle atmosphere and calculate the mass-loss rate due to Roche lobe overflow under the assumption that the atmosphere escapes through the L1 point roughly at the speed of sound. Our detailed simulations of atmospheric escape with a 1D multispecies model support this approach. The simulations also indicate that stellar XUV-driven hydrodynamic escape is consistent with energy-limited escape as long as the models properly account for the extent of the upper atmosphere where XUV radiation is absorbed and escape is initiated. By energy-limited escape here we mean mass-loss rates that are directly proportional to the net energy deposited as heat in the upper atmosphere. The heating efficiency of the upper atmosphere persists as a key uncertainty in this regime.

The main limitations of our model are the fact that they are 1D, only apply to hydrogen and helium envelopes, and do not include a possible planetary magnetic field or interaction of the planetary upper atmosphere with the stellar wind. Comparisons of 1D, 2D, and 3D models with similar physics indicate that the differences in the global mass-loss rates predicted by models with different dimensions are at most a factor of a few (Guo 2013; Owen et al. 2020), while the effect of magnetic fields and stellar winds may also be of secondary importance to thermal escape for the closest-in exoplanets. Possible differences in envelope composition may be a more significant source of uncertainty in escape models than multidimensional dynamics at this point. Ice giants in the solar system have supersolar envelope metallicities and observations of warm Neptunes and super-Earths point to the same direction for lower-mass exoplanets. Supersolar metallicity leads to higher mean molecular weight in the atmosphere and likely more efficient radiative cooling, both of which reduce the mass-loss rate. The solar system also provides another warning for modelers of exoplanet escape: the atmosphere of Pluto.

Most models of atmospheric escape that were adapted to Pluto prior to the New Horizons (NH) flyby, including the one used in this study, predicted energy-limited escape with a heating efficiency of about 30% from Pluto's extended N₂ and CH₄ atmosphere (Koskinen et al. 2015). Instead, the flyby observations revealed that the outer layers of the atmosphere are much cooler than expected, and the escape rate is 2 orders of magnitude lower and consistent with the Jeans escape regime (Young et al. 2018). The reason for the lower than expected temperature is still not well understood but it could be due to an unforeseen radiative coolant, such as supersaturated water or aerosols (Strobel & Zhu 2017; Zhang et al. 2017). While close-in exoplanets obviously need not bear any resemblance to Pluto, we ought to allow for the possibility of surprises such as this, especially if the envelope metallicity is supersolar.

In the end, future observations of exoplanet atmospheres will help to elucidate the composition of planetary envelopes and therefore their potential for mass loss over time. The extent of the hydrogen and helium envelope, as observed in the H Lyα and He i lines, appears to be relatively larger for warm Neptunes than for hot Jupiters (e.g., Ehrenreich et al. 2015; Bourrier 2018). Our model of a Uranus-like planet orbiting a Sun-like star agrees with this general idea, showing that the relative extent of the cloud of neutral hydrogen surrounding the planet is much larger than on hot Jupiters. This is because the atmosphere is more extended and outflow has a more pronounced effect on the density structure on planets with a lower gravitational potential. In contrast to hot Jupiters, we also predict that H₂ molecules should be present at large distances from the planet, along with the molecular ions ${{\rm{H}}}_{3}^{+}$ and HeH⁺. While the detection of the molecular ions may prove difficult (e.g., Chadney et al. 2016), the detection of H₂ in the FUV could provide an interesting additional test for future observations (Barthélémy et al. 2007) to determine if a hydrogen and helium envelope accurately represents the envelopes of hot Neptunes, mini-Neptunes, and super-Earths.

T.T.K. and C.H. acknowledge support by the NASA Exoplanet Research Program grant 80NSSC18K0569. R.B.F. and G.B. acknowledge support by NASA under grant 80NSSC21K0593 for the program "Alien Earths". This work benefited from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. M.E.Y. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No 805445. We thank Ilaria Pascucci for helpful discussions on planet populations and evolution.

We calculate equipotential surfaces based on the Roche potential numerically by using the iterative method adapted from Lindal et al. (1985). First, we set the reference potential at the north pole to:

$$\begin{eqnarray}&&{U}_{p}=U(z,0,0).\end{eqnarray} \tag{ A1 }$$

Then we calculate the radii along the equipotential surface at different latitudes and longitudes by using an iterative method. For each point, we choose an initial guess r₁(θ, ϕ). The subsequent error in the potential at a given point is

$$\begin{eqnarray}&&{\rm{\Delta }}{U}_{1}=U({r}_{1},\theta ,\phi )-{U}_{p}.\end{eqnarray} \tag{ A2 }$$

Based on this, the updated value r₂ is

$$\begin{eqnarray}&&{r}_{2}={r}_{1}+{\rm{\Delta }}{r}_{1},\end{eqnarray} \tag{ A3 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}{r}_{1}=\displaystyle \frac{{\rm{\Delta }}{U}_{1}}{{g}_{r}}\end{eqnarray} \tag{ A4 }$$

and

$$\begin{eqnarray}&&{g}_{r}=-\displaystyle \frac{\partial U}{\partial r}.\end{eqnarray} \tag{ A5 }$$

The iteration concludes once ΔU is vanishingly small, usually after just a few steps.

CETIMB is a 1D model of the thermosphere-ionosphere and hydrodynamic escape for close-in exoplanets. The model was described and applied to the hot Jupiter HD209458b by Koskinen et al. (2013a, 2013b). Here, we provide a revised description that includes recent updates and features specific to the current study. The model solves the coupled equations of continuity, momentum, and energy in the radial direction:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{s}}{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left[{r}^{2}{\rho }_{s}(w+{w}_{s})\right]=\sum _{t}{\rho }_{s}{R}_{{st}},\end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho w)}{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\rho {w}^{2}+p\right)=-\rho \displaystyle \frac{\partial U}{\partial r}+{F}_{\mu },\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial (\rho u)}{\partial t}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\rho {uw}\right)=\rho q-p\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}w\right)\\ \quad +\,\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\kappa \displaystyle \frac{\partial T}{\partial r}\right)+{q}_{\mu },\end{array}\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&\rho =\sum _{s}{\rho }_{s},\ w=\displaystyle \frac{1}{\rho }\sum _{s}{\rho }_{s}(w+{w}_{s}),\end{eqnarray} \tag{ B4 }$$

where ρ_s is the mass density of species s, w + w_s is the species velocity, w is the center of mass velocity, R_st is the chemical reaction rate for species s with species t, p is the pressure, U is the gravitational potential, u = c_v T is the specific internal energy, q is the sum of radiative heating and cooling rates, and κ is the coefficient of heat conduction. The viscous momentum flux and dissipation functional, respectively, are:

$$\begin{eqnarray}&&{F}_{\mu }=\displaystyle \frac{4}{3}\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\mu \displaystyle \frac{\partial w}{\partial r}\right)-\displaystyle \frac{\partial \mu }{\partial r}\displaystyle \frac{\partial w}{\partial r}-\displaystyle \frac{16}{3}\mu \displaystyle \frac{w}{{r}^{2}},\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&{q}_{\mu }=\displaystyle \frac{10}{3}\mu {\left(\displaystyle \frac{\partial w}{\partial r}\right)}^{2}-\displaystyle \frac{8}{3}\mu \displaystyle \frac{w}{r}\displaystyle \frac{\partial w}{\partial r}+\displaystyle \frac{16}{3}\mu {\left(\displaystyle \frac{w}{r}\right)}^{2},\end{eqnarray} \tag{ B6 }$$

where μ is the coefficient of viscosity. In this study, the simulations use either the normal, spherical gravitational potential (Model A) or the Roche potential for the substellar streamline (Model B):

$$\begin{eqnarray}&&\phi =-\displaystyle \frac{{{GM}}_{p}}{r}-\displaystyle \frac{{{GM}}_{s}}{a-r}-\displaystyle \frac{1}{2}{{\rm{\Omega }}}^{2}{\left(r-{M}_{{sp}}a\right)}^{2},\end{eqnarray} \tag{ B7 }$$

where M_p is the mass of the planet, M_s is the mass of the host star, a is the orbital distance,

$$\begin{eqnarray}&&{M}_{{sp}}=\displaystyle \frac{{M}_{s}}{{M}_{s}+{M}_{p}},\end{eqnarray} \tag{ B8 }$$

and Ω is the orbital angular frequency.

As indicated above, the model solves separate continuity equations for the different neutral and ion species coupled to common momentum and energy equations for the bulk atmosphere. The use of the diffusion approximation is therefore necessary to solve for the individual species velocity perturbation w_s . The diffusion velocities are obtained by matrix inversion from

$$\begin{eqnarray}&&\begin{array}{l}\left(1+{{\rm{\Lambda }}}_{s}\right)\displaystyle \frac{\partial {x}_{s}}{\partial r}+\left({x}_{s}-\displaystyle \frac{{\rho }_{s}}{\rho }\right)\displaystyle \frac{\partial (\mathrm{ln}p)}{\partial r}-\displaystyle \frac{{n}_{s}{e}_{s}E}{p}\\ \quad =\,-\sum _{t\ne s}\displaystyle \frac{{x}_{s}{x}_{t}}{{D}_{{st}}}\left({w}_{s}-{w}_{t}\right)\end{array}\end{eqnarray} \tag{ B9 }$$

under the condition that

$$\begin{eqnarray}&&\sum _{s}{\rho }_{s}{w}_{s}=0,\end{eqnarray} \tag{ B10 }$$

where x_s is the volume mixing ratios, Λ_s = K_zz /D_s , K_zz is the eddy diffusion coefficient,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{D}_{s}}=\sum _{t\ne s}\displaystyle \frac{{x}_{t}}{{D}_{{st}}}\end{eqnarray} \tag{ B11 }$$

is the net molecular diffusion coefficient, D_st is the mutual diffusion coefficient, e_s is the electric charge, and E is the radial electric field. We assume that the planetary magnetic field is negligible and ignore thermal diffusion since the coefficients for the latter are poorly known at the relevant temperatures. The diffusion coefficients, which are based on momentum transfer collision frequencies, include neutral–neutral, resonant and nonresonant ion–neutral, and electron–ion interactions. The diffusion approximation does not place constraints on the bulk (center of mass) velocity w that can be either subsonic or supersonic. It does place constraints on the magnitude of w_s , but it can be shown that a violation of the conditions in which the approximation is valid is very unlikely to occur and generally requires minor species velocities to be significantly faster than the bulk flow velocity.

In this work, we include the species H₂, H, He, ${{\rm{H}}}_{2}^{+}$, ${{\rm{H}}}_{3}^{+}$, H⁺, He⁺, and HeH⁺. While it is possible and perhaps likely that the envelopes of hot Neptunes, mini-Neptunes, and super-Earths have enhanced metallicities, with molecules and heavy atoms extending to the upper atmosphere, we do not include elements heavier than H and He. Our calculations here are primarily provided to test the common approximations and simplifications in Sections 2.3 and 2.4. The relevant chemical reactions used by the model are given in Table 1 with references. The ionization rates are calculated self-consistently by the model based on the mean solar spectrum and relevant cross sections. We note that our model does not include neutral dissociation of H₂ by photons with wavelengths longer than the ionization threshold of 80.35 nm (15.43 eV). The cross section of H₂ at these longer wavelengths consists of a multitude of electronic bands and has to be computed at high resolution (e.g., Koskinen et al. 2021). Excitation of the electronic bands leads to the emission of fluorescent radiation and dissociation, where the probability for dissociation is 0.1–0.15 (Shematovich 2010).

Table 1. Reaction Rate Coefficients

	Reaction	Rate (cm3 s−1)	References
P1	H + h ν → H+ + e	SC	Hummer (1963)
P2	He + h ν → He+ + e	SC	Yan et al. (1998)
P3	H2 + h ν → ${{\rm{H}}}_{2}^{+}\ +e$	SC	Yan et al. (1998)
P4	H2 + h ν → H+ + H + e	SC	Chung et al. (1993)
P5	H2 + h ν → H+ + H+ + e + e	SC	Dujardin et al. (1987)
R1	H+ + e → H + h ν	4.0 × 10 ${}^{-12}{\left(300/{T}_{e}\right)}^{0.64}$	Storey & Hummer (1995)
R2	He+ + e → He + h ν	4.6 × 10 ${}^{-12}{\left(300/{T}_{e}\right)}^{0.64}$	Storey & Hummer (1995)
R3	H + e → H+ + e + e	2.91 × 10 ${}^{-8}{U}^{0.39}\exp (-U)/(0.232+U),U\,=\,13.6/{E}_{e}({eV})$	Voronov (1997)
R4	He + e → He+ + e + e	1.75 × 10 ${}^{-8}{U}^{0.35}\exp (-U)/(0.180+U),U\,=\,24.6/{E}_{e}({eV})$	Voronov (1997)
R5	${{\rm{H}}}_{2}^{+}+e\to$ H + H	2.3 × 10 ${}^{-8}{\left(300/{T}_{e}\right)}^{0.4}$	Auerbach et al. (1977)
R6	${{\rm{H}}}_{3}^{+}+e\to$ H2 + H	2.16 × 10 ${}^{-8}{\left(300/{T}_{e}\right)}^{0.65}$	Larsson et al. (2008)
R7	${{\rm{H}}}_{3}^{+}+e\to$ H + H + H	5.04 × 10 ${}^{-8}{\left(300/{T}_{e}\right)}^{0.65}$	Larsson et al. (2008)
R8	${{\rm{H}}}_{2}^{+}+$ H2 → ${{\rm{H}}}_{3}^{+}$ + H	2 × 10−9	Theard & Huntress (1974)
R9	${{\rm{H}}}_{2}^{+}+$ H → H+ + H2	6.4 × 10−10	Karpas et al. (1979)
R10	H+ + H2(ν ≥ 4) → H+ + H2	10−9 $\exp (-21900/T)$	Yelle (2004)
R11	${{\rm{H}}}_{3}^{+}+$ H → ${{\rm{H}}}_{2}^{+}$ + H2	2.1 × 10−9 $\exp (-20000/T)$	Harada et al. (2010)
R12	H2 + M → H + H + M	1.5 × 10−9 $\exp (-48350/T)$	Baulch et al. (1992)
R13	H+ + H2 + M → ${{\rm{H}}}_{3}^{+}$ + M	3.2 × 10−29 n	Miller et al. (1968)
R14	H2 + e → H + H + e	1.33 × 10 ${}^{-6}{\left(300/{T}_{e}\right)}^{0.91}$ $\exp (-55800/T)$	Stibbe & Tennyson (1999)
R15	H + H + M → H2 + M	8 × 10 ${}^{-33}{\left(300/{T}_{e}\right)}^{0.6}n$	Ham et al. (1970)
R16	HeH+ + e → He + H	${10}^{-8}{\left(300/{T}_{e}\right)}^{0.6}$	Yousif & Mitchell (1989)
R17	He+ + H2 → H+ + H + He	10−9 $\exp (-5700/T)$	Moses & Bass (2000)
R18	HeH++ H2 → ${{\rm{H}}}_{3}^{+}$ + He	1.5 × 10−9	Bohme et al. (1980)
R19	HeH++ H → ${{\rm{H}}}_{2}^{+}$ + He	9.1 × 10−10	Karpas et al. (1979)
R20	He++ H2 → HeH+ + H	4.2 × 10−13	Schauer et al. (1989)
R21	H + He+ → H++ He	1.2 × 10 ${}^{-15}{\left(300/T\right)}^{-0.25}$	Stancil et al. (1998)
R22	H+ He → H + He+	1.75 × 10 ${}^{-11}{\left(300/T\right)}^{0.75}\exp (-{\rm{128,000}}/T)$	Glover & Jappsen (2007)
R23	H2 + ${\mathrm{He}}^{+}\to {{\rm{H}}}_{2}^{+}$ + He	7.2 × 10−15	Barlow (1984)

Note. SC: Calculated self-consistently based on ionization cross sections, stellar flux, and density profiles in the model.

Download table as:  ASCIITypeset image

In our models of hot Jupiters, we have generally assumed that thermal dissociation of H₂ (reaction R12) and dissociation of H₂ due to photochemistry in the middle atmosphere (Moses et al. 2011; Lavvas et al. 2014) are more important than dissociation of the electronic excited states. The same is not necessarily true for Neptunes with cooler atmospheres. The extension of the model to properly include the high resolution H₂ cross section and detailed calculations of neutral dissociation, and the related heating rate requires more research and is beyond the scope of this work. However, we used a cross section based on the laboratory work of Backx et al. (1976) with a dissociation probability of 0.125 in the model to check the potential of neutral dissociation to change our results. Since this cross section has very low wavelength resolution, it is likely to overestimate the dissociation rate (e.g., see Figure 1 in Kim et al. 2014). Figure 13 shows the resulting composition for Model A at 0.05 au that should be compared with Figure 8 for the model that does not include neutral dissociation of H₂. The abundance of H₂ is somewhat reduced in the revised model and the mass-loss rate is slightly higher but the results are qualitatively similar and the mass-loss rates agree to within a factor of 1.4. We have therefore decided to postpone the inclusion of detailed H₂ calculations to future work.

Zoom In Zoom Out Reset image size

We also do not include direct ionization of H from the n = 2 state. This ionization pathway and related heating by the stellar flux in the Balmer continuum are potentially important for hot Jupiters (Huang et al. 2017; García Muñoz & Schneider 2019). In order to test the possible significance of the n = 2 state to our results, we used the detailed balance model of Huang et al. (2017) to calculate the 2p and 2s excited state populations and the related ionization and heating rates for the simulations presented in Section 3.2. In all cases, the excitation of H is dominated by the diffuse Lyα radiation field that arises from multiple resonant scattering of the stellar Lyα flux in the atmosphere. This poses a complex problem that we solve by using a Monte Carlo approach. We found the heating rate by the Balmer continuum to be negligible. The ionization rate, however, becomes comparable to and exceeds the ionization rate from the ground state by XUV radiation in a narrow region near the bottom of the model at 0.03 au (see Figure 14). The impact of this additional ionization on our models, however, is small, and in reality, heavier molecules that might dominate the chemistry instead are likely to be present in this region. We have therefore also postponed the coupling of the escape models to the detailed balance model for hydrogen level populations to future work.

Zoom In Zoom Out Reset image size

The prevalence of neutral H and outflow of H₂ at high altitudes around hot Neptunes differ from similar models of hot Jupiters. This leads to large transit depths at H Lyα (Ehrenreich et al. 2015; Loyd et al. 2017a) and may enable the future detection of H₂ on exoplanets (e.g., Barthélémy et al. 2007). To investigate the mechanism that enables these features, we analyze the terms of the continuity equation in detail for Model A at 0.05 au. For example, Figure 15 shows the continuity equation terms for neutral H. At the lowest altitudes in the model, the chemical loss and production rates (ionization and recombination reactions) are in rough balance but above the heating peak, the ionization rate exceeds the recombination rate. Instead of recombination, ionization of H is balanced by upward transport of neutral H from below by escaping gas. In other words, outflow replenishes the outer layers of the atmosphere with neutral H from below and keeps the upper atmosphere mostly neutral. The remarkable numerical accuracy of the model is highlighted by the near perfect alignment of the proton production rate (solid black line) with the sum of the proton loss and transport rates (solid blue line).

Zoom In Zoom Out Reset image size

We find a similar result for H₂, as illustrated by Figure 16, which shows the continuity equation terms for H₂. Chemical loss of H₂ (dissociation and dissociative ionization reactions) dominates over production but outflow brings fresh H₂ up from deeper layers of the atmosphere and primarily balances the total loss rate. In the end, the sum of the chemical loss and transport terms equals the production rate, as expected in a steady state model. These remarkable results highlight the fact that escape is not a trivial process that can be ignored or simply parameterized in models that attempt to fit, say, exoplanet transmission spectra. Instead, it can fundamentally alter the ionization balance and relative abundances of different species in the upper atmosphere.

Zoom In Zoom Out Reset image size

The sum of the radiative heating and cooling rates (q in Equation (B3)) includes heating due to photoionization of the atmosphere by stellar XUV radiation and cooling by recombination, infrared emission by ${{\rm{H}}}_{3}^{+}$, and H Lyα emission. We calculate stellar XUV energy deposition self-consistently based on the density profiles in the model every 10 time steps, where one time step is 1 s. Each ionization event by a photon with frequency ν releases an energy E = h ν − I, where I is the ionization potential. In line with many previous studies (e.g., Murray-Clay et al. 2009), we assume that 100% of this energy heats the atmosphere and that the energy I is lost in subsequent radiative recombination events or chemistry. In addition, we account for the loss of electron thermal energy in recombination of H⁺ (the corresponding cooling rate due to the recombination of other ions is small). We use the H Lyα cooling rate from Glover & Jappsen (2007), which is similar to the widely adapted rate coefficient from Murray-Clay et al. (2009), with the exception of a small temperature-dependent correction factor. We assume optically thin infrared cooling to space by ${{\rm{H}}}_{3}^{+}$. The per molecule cooling rate is based on line lists (Neale et al. 1996; Miller et al. 2013) and a correction factor for non-LTE conditions that is derived from detailed balance calculations (Koskinen et al. 2009).

The model extends from the bottom of the thermosphere to several planetary radii. There are 410 radius levels on a stretched grid that mimics the expected spreading of pressure levels with altitude:

$$\begin{eqnarray}&&{r}_{n}={r}_{1}+\displaystyle \frac{\delta {r}_{1}(1-{s}^{n})}{1-s},\end{eqnarray} \tag{ B12 }$$

where n is the grid index of level n, r is the radius, r₁ is the radius at the lower boundary, δ r₁ = 10 km is the grid spacing at the bottom, and s = 1.014 is the stretch factor. The lower boundary conditions are the temperature, pressure (p₁ = 1 μbar ), and species mixing ratios that are based on our lower atmosphere models. The upper boundary conditions are adaptable. If an exobase is found within the grid, the model uses modified Jeans escape conditions for individual species velocities, densities, and energy flux (Koskinen et al. 2015) that incorporate the ambipolar electric field for ions (Koskinen et al. 2013a) and Roche potential. If the exobase is above the sonic point, L1 point, or the grid, the model uses extrapolated boundary conditions appropriate for hydrodynamic escape (e.g., Tian et al. 2005). The model checks for the presence of an exobase every time step based on the collision frequencies between different species. All simulations in this paper are consistent with hydrodynamic escape.