Removal of Hot Saturns in Mass–Radius Plane by Runaway Mass Loss

Mass-loss processes transform large planets into their smaller counterparts and leave imprints in the observed exoplanet population in the shape of deficits. The most famous example is the radius gap that separates the gas-enveloped mini-Neptunes (∼2–4 R_⊕) from the more rocky super-Earths (∼1–1.7 R_⊕), which was predicted to exist by theories of photoevaporative mass loss (e.g., Lopez et al. 2012; Owen & Wu 2013) and later confirmed observationally when the host star samples were confined to those that were well characterized by high-resolution spectra (e.g., Fulton et al. 2017; Fulton & Petigura 2018), asteroseismology (e.g., Van Eylen et al. 2018), or followed up by Gaia and detailed statistical analysis (e.g., Hsu et al. 2019). ⁷

Another deficit we see in exoplanet population is the sub-Jovian desert, a roughly triangular region in mass–period and radius–period space where we see no Saturn size objects inside orbital periods of a few days, first reported in mass space by Szabó & Kiss (2011) then in radius space by Beaugé & Nesvorný (2013) and then later confirmed to exist in both mass–period and radius–period spaces by Mazeh et al. (2016). This sub-Jovian desert subsumes the so-called "photoevaporation desert," which refers to the lack of short-period mini-Neptunes whose incident flux would correspond to ∼650 times that received by the Earth (e.g., Lundkvist et al. 2016), suggesting photoevaporative mass loss to be at play in sculpting at least part of the sub-Jovian desert. Yet another way to view the desert is to consider the orbital period distribution of sub-Saturns (4–8 R_⊕). Using a hydrodynamic mass-loss model of Kubyshkina et al. (2018), Hallatt & Lee (2022) demonstrated that mass loss alone could reproduce the falloff in the sub-Saturn occurrence rate at short periods with the stripped Saturns accounting for a fraction of the observed mini-Neptunes and super-Earths.

Although photoevaporation is expected to be an efficient mechanism to strip small planets (≲8 R_⊕) of their gaseous envelopes, detailed radiative-hydrodynamic calculations have found that gas giants are resilient against mass loss, losing less than 1% of its total mass over the star's main-sequence lifetime (e.g., Murray-Clay et al. 2009; Owen & Jackson 2012). It follows that although the lower boundary of the sub-Jovian desert is likely shaped by photoevaporative mass loss, the upper boundary of the desert may require a different explanation.

Valsecchi et al. (2014) proposed that hot Jupiters that are excited onto high-enough eccentricities by, e.g., Kozai can undergo catastrophic Roche lobe overflow and lose all their envelopes, transforming into rocky super-Earths at ultrashort periods. While some of the short-period rocky planets may be stripped giants, the lack of strong correlation between the occurrence of ultrashort-period planets and the host star metallicity (as opposed to the strong correlation between the occurrence of hot Jupiters and high metallicity of the host star) suggest that this is likely not the main origin channel (Winn et al. 2018).

Instead of mass loss, Matsakos & Königl (2016) considered the boundaries of the desert to be traced by the history of dynamical migration. They proposed that both the upper and the lower boundaries can be reproduced by planets arriving at their closest-in orbits by high-eccentricity migration with the final orbital distance set by either the Roche lobe radius or to twice the pericenter distance (of the initially excited eccentric orbit) as expected from secular chaos (Wu & Lithwick 2011). The different shape of the upper and the lower boundary (in the mass–period space) reflects the different mass–radius relationship followed by smaller versus larger planets. Owen & Lai (2018) combined the expectation from high-eccentricity migration (along with tidal decay) with photoevaporative mass loss and concluded that both mechanisms are required to explain the overall shape with the former being responsible for the upper boundary and the latter being responsible for the lower boundary.

All aforementioned literature relied on their assumed mass–radius relationship for different categories of planets. In this Letter, we investigate directly the sub-Jovian population in the radius–mass space and self-consistently compute the evolution of planetary interiors under mass loss. We are motivated by a clean upper boundary in the radius–mass plane that is well approximated by a constant ∼0.1 g cm⁻³ (see the bottom panel of Figure 1), which also coincides with the region of the parameter space where highly irradiated planets become very large and the mass–radius relationship steepens (Thorngren & Fortney 2018), suggestive of photoevaporative mass loss playing a major role in creating the boundary. By comparison, this density boundary manifests as the very upper edge of the hot Jupiter population in the radius–period space (middle panel) while being coincident with the usual sub-Jovian desert in the mass–period space (top panel).

Zoom In Zoom Out Reset image size

This paper is organized as follows. In Section 2, we describe how we build time-evolving models of interior structure accounting for X-ray and extreme UV irradiation (XUV)-driven mass loss. Our results are presented in Section 3 and we conclude in Section 4.

We start with gathering data on the observed exoplanets, which we collect from the NASA Exoplanet Archive (Akeson et al. 2013; NASA Exoplanet Science Institute 2022) and Exoplanet.eu (Schneider et al. 2011). We combine these data and select confirmed planets with radial-velocity mass measurements more precise than 30% and radius measurements more precise than 20%. Additionally, we exclude Kepler-13b, Kepler-435b, Kepler-87b, and KELT-9b as their quoted stellar properties between listed references are inconsistent or their quoted flux is so high that they cannot be self-consistently described within our model framework. While we are primarily interested in planets between 0.1 and 1M_J and with a < 0.1 au, we include planets outside these ranges for context. It is not our goal to reproduce the absolute abundance of planets and so we do not correct for observational biases. Rather, what is relevant is whether our model is able to account for the paucity of planets in spaces where they are favorable to be observed (i.e., large radii and short period).

The majority of our interior structure evolution model is the same as was used in Thorngren & Fortney (2019). These are one-dimensional models constructed by solving the equations of hydrostatic equilibrium, mass conservation, and an equation of state (EOS):

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial m}=-\displaystyle \frac{{Gm}}{4\pi {r}^{4}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial r}{\partial m}=\displaystyle \frac{1}{4\pi {r}^{2}\rho },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\rho =\rho (P,T),\end{eqnarray} \tag{ 3 }$$

where m and r are the enclosed mass and radii of each layer respectively, P is pressure, ρ is density, T is temperature, and G is the gravitational constant. We use the Chabrier et al. (2019) EOS for a solar-ratio mixture of hydrogen and helium, and a 50–50 mixture of rock and ice for the metals (Thompson 1990).

To set the bulk planet metallicity, we use the mass–metallicity distribution from Thorngren et al. (2016), in which the logarithm of the metal mass M_z is normally distributed with a median M_z = 0.182M^0.61. This median value is our default; in Section 3, we experiment with varying the bulk metallicity. The offset from the mean (in log-metallicity) is parameterized by Z_zpl the z-score (number of sigma away from the mean)—e.g., Z_zpl = 1 selects the 84th percentile. This is the same approach as was used in (Thorngren & Fortney 2019). This is equivalent to multiplying the median by a factor of ${10}^{\sigma {Z}_{\mathrm{zpl}}}$, where σ = 0.26. The advantage of this approach is that drawing Z_zpl from a standard normal distribution correctly reproduces the mass–metallicity dispersion. The core mass is set to half the metal mass, but with a minimum of 5 M_⊕ and a maximum of 15 M_⊕; any remaining metals are mixed into the envelope.

Modeling the thermal evolution of the planet is slightly more complex than in past cases because the mass-loss process can also affect the thermal state of the planet. Our model tracks the thermal state of the planet in terms of the envelope specific entropy s, so we seek to calculate its rate of change with time in order to integrate. We start with conservation of total energy E, split into the change in total energy in the core E_c and the change in total energy in the envelope E_e :

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=\displaystyle \frac{{{dE}}_{c}}{{dt}}+\displaystyle \frac{{{dE}}_{e}}{{dt}},\end{eqnarray} \tag{ 4 }$$

where t is time. We assume an isothermal core, and that the core can release energy into the envelope rapidly enough that the core-envelope boundary has no temperature jump. We adopt c = 7.5 × 10⁶ erg g⁻¹ K⁻¹ (chosen as a typical value for rock, see Waples & Waples 2004) for the specific heat capacity of the core. The core temperature T_c can then change either because the envelope specific entropy changes or because the envelope loses mass and the core-envelope boundary intersects at a lower pressure and therefore temperature (even at constant entropy). We find that this core depressurization heat is a negligible effect, but we include it for completeness. The energy change in the envelope is readily calculated from the change in entropy (see, e.g., Thorngren et al. 2016). Putting this all together and solving for the change in entropy with time, we have

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}={{cM}}_{c}\left(\displaystyle \frac{\partial {T}_{c}}{\partial M}\displaystyle \frac{\partial M}{\partial t}+\displaystyle \frac{\partial {T}_{c}}{\partial s}\displaystyle \frac{\partial s}{\partial t}\right)+\displaystyle \frac{{{dE}}_{e}}{{ds}}\displaystyle \frac{{ds}}{{dt}};\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{ds}}{{dt}}=\left(\displaystyle \frac{{dE}}{{dt}}-{{cM}}_{c}\displaystyle \frac{\partial {T}_{c}}{\partial M}\displaystyle \frac{\partial M}{\partial t}\right)/\left(\displaystyle \frac{\partial {E}_{e}}{\partial s}+{{cM}}_{c}\displaystyle \frac{\partial {T}_{c}}{\partial s}\right),\end{eqnarray} \tag{ 6 }$$

where M is the total mass of the planet, and M_c is the mass of the core. These derivatives are numerically computed by constructing and evaluating the differences between two static models at slightly different envelope masses and at slightly different specific entropy for ∂T_c /∂M and ∂T_c /∂s, respectively. The change in envelope energy with entropy is calculated by integrating ∂s/∂E = 1/T(P, s) with respect to mass across the envelope (as in, e.g., Lopez et al. 2012) with T(P, s) supplied by the adopted EOS.

We now consider the energy sources and sinks and equate them in Equation (4):

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}={L}_{\mathrm{rad}}+\pi {R}^{2}(\epsilon F-4{F}_{\mathrm{int}})\end{eqnarray} \tag{ 7 }$$

where L_rad is radioactive heating, R is the radius of the planet, is the anomalous heating efficiency for which we adopt the median inferred value from Thorngren & Fortney (2018), F is irradiation flux, and F_int is the internal cooling flux of the planet. Radioactive heating is computed from the decay of U235, U238, Th232, and K40 following Nettelmann et al. (2011) and Lopez et al. (2012), using Anders & Grevesse (1989) for the radioactive elemental abundances relative to silicon. We find that L_rad is an extremely minor component of the energy balance. Finally, F_int is calculated for a given envelope specific entropy, surface gravity, and insolation using the nongray atmosphere models of Fortney et al. (2007).

The anomalous heating is a key feature of our models. Their inclusion is motivated by the fact that hot Saturns are above the T_eq = 1000 K inflation threshold (Demory & Seager 2011; Miller & Fortney 2011) and that many of them are too large to be explained by a low metallicity alone (R > 1.2R_J ; see Thorngren et al. 2016). Our adoption of Thorngren & Fortney (2018) inflation power fits—which were based on planets with M > 0.5M_J —is reasonable because most of the planets below that regime fit the models quite well; the cut was only made to avoid the influence of a few dense outliers such as HD 149026 b (see, e.g., Sato et al. 2005; Fortney et al. 2006). ⁸

To model mass loss, we use the following expression:

$$\begin{eqnarray}&&\dot{M}\approx \displaystyle \frac{\eta \pi {F}_{\mathrm{XUV}}{R}_{\mathrm{XUV}}^{3}}{{{GMK}}_{t}}\approx \displaystyle \frac{3}{4}\displaystyle \frac{\eta {F}_{\mathrm{XUV}}}{{{GK}}_{t}{\rho }_{\mathrm{XUV}}};\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{K}_{t}=1-\displaystyle \frac{3}{2\xi }+\displaystyle \frac{1}{2{\xi }^{3}};\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\xi =\displaystyle \frac{{R}_{\mathrm{hill}}}{{R}_{\mathrm{XUV}}},\end{eqnarray} \tag{ 10 }$$

where η is the mass-loss efficiency factor, K_t is a tidal correction factor, R_hill is the Hill radius, R_XUV is the radius below which the planet becomes opaque to XUV irradiation, which we set at 10 nanobars (a conservative limit of Lopez et al. 2012), and F_XUV is the extreme ultraviolet flux from the parent star, which we interpolate from the stellar evolution grid of Johnstone et al. (2021) using the median stellar spin. While Equation (8) mirrors that of energy-limited approximation from Watson et al. (1981), Lopez et al. (2012), and Lopez & Fortney (2013), we do not assume energy-limited mass loss. Instead we use the results of Caldiroli et al. (2022), who apply an empirical fitting to a detailed hydrodynamics code ATES (Caldiroli et al. 2021) that encompasses both energy-limited and non-energy-limited regime of mass loss, reported as an effective mass-loss efficiency η that varies nontrivially with incident flux and planetary potential. We find that hot Saturns fall in their low-gravity regime where η generally ranges in the few tens of percent. The deposition depth of XUV irradiation is another important parameter for XUV-driven mass-loss models; it is often quoted as a few nanobars (Murray-Clay et al. 2009; Caldiroli et al. 2021) for giants, although they may be deeper in depending on the mass of the planet (Ionov et al. 2018; Owen & Lai 2018). We do not need to set this parameter, as ATES distributes the energy according to the ionization cross sections in the model atmosphere. Finally, we have rewritten the mass-loss rate $\dot{M}$ in terms of ρ_XUV the planet bulk density using the XUV radius to highlight the connection between the observed density boundary for hot Saturns (Figure 1) and XUV-driven mass loss.

For planets that grow to overflow their Roche lobes, we model both the mass loss from the overflow as well as the effect on the planets' semimajor axes, which expand in proportion to the mass lost from the conservation of angular momentum (Valsecchi et al. 2014):

$$\begin{eqnarray}&&\displaystyle \frac{\dot{a}}{a}=-2\displaystyle \frac{\dot{M}}{M}\left(1-\displaystyle \frac{M}{{M}_{* }}\right).\end{eqnarray} \tag{ 11 }$$

Here, a is the planet's semimajor axis and $\dot{a}$ is the rate of change thereof, and M_* is the mass of the parent star. During Roche lobe overflow, we assign the mass rate of change $\dot{M}$ to an arbitrary exponential decay with a half-life of 13.8 Myr. We choose this scheme instead of directly modeling the detailed physics of Roche lobe overflow (e.g., Valsecchi et al. 2014) since the true rate of Roche lobe overflow would be fast compared to the lifetime of planets and their thermal evolution. In fact, the net effect of Roche lobe overflow in our models is to cause the planet to rapidly lose as much mass as is necessary to expand the orbit to the point where it is no longer overflowing. This leads to an interesting interplay with XUV-driven mass loss, which we will discuss in the next section.

Figure 2 illustrates our radius and mass evolution tracks. We show three representative cases for a Saturn-mass planet. At Saturn's orbit (9.58 au), the planet experiences no mass loss and cools normally. At 0.07 au, the planet is a fairly typical hot Saturn, which cools to an inflated equilibrium where the luminosity equals the anomalous heating (see, e.g., Thorngren et al. 2021). The planet experiences modest mass loss over the course of its lifetime.

Zoom In Zoom Out Reset image size

Placing a Saturn at 0.05 au has a very different outcome. Although the planet does not initially overflow its Roche lobe nor is the initial mass-loss rate particularly high, the slow erosion of mass accelerates into the catastrophic removal of the envelope at ∼400 Myr. This runaway mass loss is effected by an adiabatic expansion of the planet as it enlarges while losing mass at constant specific entropy. Critically, the rate of mass loss scales as the radius cubed (recall Equation (8)), enhanced further by the tidal factor K_t . Furthermore, at very low densities the mass–radius relationship becomes particularly steep (see Figure 1). The result is a positive feedback loop that reaches a head at around 100 Myr, when the planet grows in size rapidly, overflows its Roche lobe, rapidly losing mass and expanding the semimajor axis. Over the course of 30 Myr., the planet goes from 58 to 17.4 M_⊕ and moves out to 0.95 au. What's left of the envelope continues to be stripped away by XUV-driven mass loss, leaving less than 1 M_⊕ of H/He on top of the 10.67 M_⊕ core by 1 Gyr.

There is a subtle feature of this process which Valsecchi et al. (2014) observed and which we reproduce. After Roche Lobe overflow starts, the orbit quickly expands just far enough to shut it off again (sacrificing some mass in the process). However, unlike in Valsecchi et al. (2014), XUV mass loss is still taking place, and again expands the radius into the overflow regime. Thus the planet is continuously very close to filling its Roche lobe and periodically exceeding it. We do not explore this process in great detail here, as for our purposes it is enough to know that the process is very fast compared to thermal evolution, and that the growth in the semimajor axes is directly proportional to the mass lost.

Our result of runaway mass loss is reminiscent of what was reported by Kurokawa & Nakamoto (2014), which has since then been criticized for their use of fixed η and overestimating the mass-loss rate by fixing the pressure at the bottom of the photoevaporative flow at 1 nanobar (Ionov et al. 2018; Owen & Lai 2018). The latter authors did not find any evidence of runaway mass loss. Our study differs from Kurokawa & Nakamoto (2014) in that we rigorously account for the variable η, but it was mainly our treatment of the thermal evolution that allowed us to recover catastrophic erosion. Compared with, e.g., Owen & Lai (2018), we include the effect of anomalous heating throughout the planet's thermal evolution instead of at the very end, which keeps our planets puffy. We have verified that when this anomalous heating is set to zero, the mass loss never reaches a runaway state.

A related issue we encounter is that during runaway, the gravitational potential at the planet's atmosphere drops below what is modeled by Caldiroli et al. (2022), necessitating extrapolation. The actual efficiency η at low masses may be somewhat smaller than what we have adopted here (see, e.g., Shematovich et al. 2014; Ionov et al. 2018). In the next section we will see that fortunately, our results are robust against even very small values of η.

3.1. Carving Out the Edges

We are particularly interested in determining the overall effect of mass loss, especially the runaway variety, on the hot Saturn population, irrespective of particular initial distribution of host Saturns arising from planet formation, which is outside the scope of this paper. Instead, we seek to answer the more narrow question: if a hot Saturn of some mass is placed at a particular semimajor axis, can it survive?

We have run our mass-loss model on a grid of masses and semimajor axes. The masses were 15 values log-uniformly spaced from 20M_⊕ to 1.5 M _J and the semimajor axes were 0.03, 0.04, 0.05, 0.07, and 0.1 au. The results are shown for several slices in time on Figure 3—an animation is also available. We assume a solar-mass star with a median amount of XUV irradiation (i.e., median spin) for its age. Planets on short orbits experience more anomalous heating and consequently maintain a larger radius than those at longer periods after they cool. Short-period, low-density planets are particularly susceptible to runaway mass loss and we observe a clear boundary at ∼0.1 g cm⁻³ within ∼100 Myr, which persists to ∼1 Gyr, matching the observed density boundary of hot Saturns in radius–mass space.

Download figure:

Video Standard image High-resolution image

We have also considered the effect of stellar mass, the planet metallicity, and the effective XUV flux on the mass-loss rates (Figure 4). More-massive stars have higher bolometric and XUV luminosities so planets around them are more inflated for a given mass and given orbital distance and therefore more susceptible to mass loss. The net result is that the closest-in planets (0.03 au, triangle points) around higher mass stars lose their envelopes and lose their hot Saturn status while those at longer orbital distances survive against mass loss but stay puffy. Similarly, planets of lower bulk metallicity also appear puffy (due to their lower densities) and so the closest-in planets are lost to mass loss (note the lack of triangle points in the bottom center panel). Those at wider separations are stable against mass loss but will be larger than planets of higher bulk metallicity.

Zoom In Zoom Out Reset image size

In Figure 3, we see excess of planets of mass ≳140M_⊕ with radii above 1.5 R_J in the observations compared to our fiducial model. Our parameter study (Figure 4) suggests that such planets are likely around more-massive stars and/or they have lower bulk metallicity. Indeed, planets in our sample above 1.5 R_J all orbit around stars more massive than the Sun, with an average mass of 1.41 M_⊙ (versus ∼1.1 M_⊙ for 1 R_J < R < 1.5 R_J ), corroborating our results.

The right panels of Figure 4 illustrate the effect of varying η. As expected, closest-in planets are lost for higher fixed mass-loss efficiency (note the lack of triangle points in the top right panel). Even at a low η = 0.01, the low-density desert is largely cleared out by 1 Gyr, demonstrating that ∼0.1 g cm⁻³ is the boundary for a runaway mass loss (i.e., when the planet enters the runaway phase, they will be stripped of their sub-Saturn status regardless of the details of how they arrived to such a boundary).

We now explore how the runaway mass loss manifests in the mass–period and radius–period spaces. First, we draw 400 sets of mass, semimajor axis, stellar mass, and metallicity, distributed log-uniformly from 20 to 1.5 M _J in mass, uniformly from 0.8 to 1.2 M_⊙ in stellar mass, log-uniformly from 0.03 to 0.1 au in the semimajor axis, and the metallicity according to the mass–metallicity relation (Thorngren et al. 2021). These planets are evolved out to 5 Gyr and the results are shown in Figure 5. We observe that most of the low-mass, short-period planets have large radii and are quickly reduced to stripped cores. All model planets near or above 2R_J are destroyed by runaway inflation within 100 Myr, shaping the upper edge of the radius distribution. Over 5 Gyr, the sub-Jovian desert in mass–period space is largely carved out by mass loss although it is not fully cleared, mostly due to high-bulk-metallicity planets that have small radii, suggesting the formation of hot Saturns disfavors those with high densities. In some cases, tidal inspiral may lead to runaway mass loss, though careful calculation would be needed to examine this possibility as the XUV luminosity will have dropped significantly on inspiral timescales (∼Gyr).

Download figure:

Video Standard image High-resolution image

Another interesting feature of our evolution models is that the planet population appears to pass through a few phases in time. In the first few tens of Myr, many planets are actively losing significant amounts of mass while undergoing orbital expansion; they have not reached their final mass or size yet. At a few hundred Myr, the majority of planets have lost most of the mass that they are going to lose during their lifetime, but their envelopes are still too hot and remain puffy, giving rise to a large number of puffy planets (R > 0.8R_J ) even at relatively long periods and low mass (see Figure 5, middle column). By around 500 Myr, most planets complete the runaway mass loss, though some continue to runaway late. The population of early-runaway planets has cooled to below 0.5 R_J and will continue to cool toward the evolved radius of around 0.2–0.3 R_J over the subsequent Gyr. The existence of young puffy planets (∼100 Myr, ∼0.5R_J , ≲0.03M_J ) beyond ∼4 days is a unique prediction of our model that could be tested by, e.g., TESS and follow-up mass measurements (such as, e.g., Newton et al. 2019; Zhou et al. 2021).

From our simple initial random distribution of planets, we find that the ones that survive against mass loss to fill in the sub-Jovian desert in mass–period and radius–period spaces are those that have high bulk metallicity (and so higher density). Our result agrees with the discovery of dense Neptune to Saturn size planets found within the desert such as LTT 9799 b (Jenkins et al. 2020) and TOI-849 b (Armstrong et al. 2020), although these two planets could also be stripped massive cores of even larger giants.

Our primary conclusion is that hot Saturns can undergo runaway mass loss caused by positive feedback between the adiabatic expansion of the radius and the mass loss (Figure 2). In many cases, this leads to Roche lobe overflow, though this is not necessary to remove the planets' envelope. Over a range of stellar masses, planet bulk metallicity, and the XUV-mass-loss efficiency, we find that this process can robustly sculpt the low-density boundary at ∼0.1 g cm⁻³ (Figures 3 and 4).

In radius–period space, the radii for our initial distribution already show a relative paucity of planets at ∼0.3–1.0R_J inside ∼3 days, similar to what is observed. This initial radius–period paucity arises from a sudden rise in radius over a small dynamic range in mass at high incident flux (see the model mass–radius relation curves in the bottom panel of Figure 1). Subsequent XUV-driven mass loss quickly removes planets from the sub-Jovian desert in mass–period space, though somewhat too many planets appear to remain in the desert compared to the observations (Figure 5). To completely clear out the desert, we may need to more carefully account for Roche lobe overflow and tidal processes (e.g., Valsecchi et al. 2015) in our models, but perhaps more important is the relatively simple choice of input distribution, especially on planet mass and period. For example, Hallatt & Lee (2022) were able to reproduce the desert (in terms of the decline in sub-Saturn occurrence rate at short periods) through XUV-driven mass loss and a careful consideration of the initial core mass function and gas accretion theory. It also remains possible that some of the short-period giants arrived there via high-eccentricity migration, which contribute to shaping the upper edge of the desert in the mass–period space (Matsakos & Königl 2016; Owen & Lai 2018).

Planets that lose their envelopes collect around ∼0.2 to ∼0.3 R_J based on their initial core mass, which we imposed in Section 2. This assumption does not change the triggering of runaway mass loss, but it does impact the final state. This suggests that if we are able to identify probable stripped cores, they would be a unique window into the deep interiors of giant planets. A few such candidates have already been proposed (Petigura et al. 2017; Armstrong et al. 2020; Jenkins et al. 2020), and a targeted study would likely identify more.

Our results differ somewhat from Owen & Lai (2018) and Ionov et al. (2018), who argue that photoevaporative mass loss is unable to explain the upper edge of the sub-Jovian desert (i.e., hot Saturns are resilient against mass loss). First, we argue that the effect of mass loss is more clearly seen in radius–mass space rather than mass–period or radius–period. Second, and more importantly, we include the anomalous heating of hot giants found in Thorngren & Fortney (2018) throughout the planets' evolution, resulting in larger planets more vulnerable to mass loss than Owen & Lai (2018), who evolve their planets without heating and then inflate them after mass loss has completed. We believe our approach is more consistent with the observations, as in Thorngren et al. (2021), it was found that young planets are as large or larger than their older equivalents once the present-day incident flux is controlled for.

Recently, Vissapragada et al. (2022) examined a set of seven transiting planets on the (mostly upper) edge of the desert in mass–period space. They found that none of the planets were undergoing massive mass loss despite sizable XUV fluxes. We do not find their result to be in tension with ours as all the planets they investigated have densities larger than the 0.1 g cm⁻³ boundary, ranging from 0.172 to 0.717 g cm⁻³. The lowest-density planet, WASP-177 b is still interesting as a planet near the density edge, but they find that its mass-loss rate is quite low, which could be due to its old age of ∼10 Gyr (Turner et al. 2019), as XUV luminosity declines as the star ages.

Future work on this area could investigate the potential mass-loss histories of giant planets near the low-density boundary. One good candidate for this is WASP-17 b (Anderson et al. 2010), the largest known planet with a well-determined mass and radius, and one which is likely to be losing mass at a prodigious rate. Qatar-10 b (Alsubai et al. 2019) and WASP-76 b (West et al. 2016) are also good cases for substantial ongoing mass loss. Also of interest are the smaller and lower-mass planets that are nevertheless near the 0.1 g cm⁻³ boundary, including WASP-20 A b (Anderson et al. 2015) and WASP-153 b (Demangeon et al. 2018). While these planets are not experiencing runaway mass loss at the present day, they likely are experiencing some loss, and they may be nearing runaway. Measuring both the mass-loss rates and XUV fluxes of their parent stars could help to pinpoint some of the fine details of this process.

We thank Andrea Caldiroli, Jonathan Fortney, Ruth Murray-Clay, and Shreyas Vissapragada for helpful discussions, and the anonymous reviewer for thoughtful suggestions and comments. D.P.T thanks the Trottier Institute for Research on Exoplanets (iREx) for support via the Trottier Fellowship and Johns Hopkins University for support via the Davis Fellowship. E.J.L. gratefully acknowledges support from NSERC, FRQNT, the Trottier Space Institute, and the William Dawson Scholarship from McGill. E.D.L. would like to acknowledge support from the GSFC Sellers Exoplanet Environments Collaboration (SEEC), which is funded in part by the NASA Planetary Science Divisions Internal Scientist Funding Model.