Research Article
Open access
Published Online: 15 January 2015

Habitable Evaporated Cores: Transforming Mini-Neptunes into Super-Earths in the Habitable Zones of M Dwarfs

Publication: Astrobiology
Volume 15, Issue Number 1

Abstract

We show that photoevaporation of small gaseous exoplanets (“mini-Neptunes”) in the habitable zones of M dwarfs can remove several Earth masses of hydrogen and helium from these planets and transform them into potentially habitable worlds. We couple X-ray/extreme ultraviolet (XUV)–driven escape, thermal evolution, tidal evolution, and orbital migration to explore the types of systems that may harbor such “habitable evaporated cores” (HECs). We find that HECs are most likely to form from planets with ∼1 M solid cores with up to about 50% H/He by mass, though whether or not a given mini-Neptune forms a HEC is highly dependent on the early XUV evolution of the host star. As terrestrial planet formation around M dwarfs by accumulation of local material is likely to form planets that are small and dry, evaporation of small migrating mini-Neptunes could be one of the dominant formation mechanisms for volatile-rich Earths around these stars. Key Words: Astrobiology—Extrasolar terrestrial planets—Habitability—Planetary atmospheres—Tides. Astrobiology 15, 57–88.

1. Introduction

Due to their small radii and low luminosities, M dwarfs currently offer the best opportunity for the detection of terrestrial planets in the habitable zone (HZ), the region around a star where liquid water can exist on the surface of a planet (Kasting et al., 1993; Kopparapu et al., 2013). Upcoming missions such as the Transiting Exoplanet Survey Satellite (TESS) and the repurposed Kepler spacecraft (K2) will be capable of detecting potentially habitable Earths and super-Earths around M dwarfs (Ricker et al., 2010; Howell et al., 2014). In particular, the detection of potentially habitable planets around M dwarfs is easier because the HZs of these stars can be as close in as ∼0.02 AU (Kopparapu et al., 2013). However, such proximity implies that terrestrial planets forming within the HZs of low-mass stars are likely to be small (≲0.3 M) and form dry (Lissauer, 2007; Raymond et al., 2007). Moreover, M dwarfs are extremely active when young, bombarding their planets with high-energy radiation and bursts of relativistic particles during flaring events, which can erode their atmospheres and potentially sterilize the surface (Scalo et al., 2007; Segura et al., 2010). Strong tidal heating and orbital evolution may further impact the habitability of planets around these stars (Barnes et al., 2013). Many planets formed in situ in the HZs of M dwarfs may therefore be uninhabitable.
However, planets need not form and remain in place. It is now commonly accepted that both disk-driven migration and planet-planet interactions can lead to substantial orbital changes, potentially bringing planets from outside the snow line (the region of the disk beyond which water and other volatiles condense into ices, facilitating the formation of massive planetary cores) to within the HZ (Hayashi et al., 1985; Ida and Lin, 2008a, 2008b; Ogihara and Ida, 2009). Disk-driven migration relies on the exchange of angular momentum between a planet and the surrounding gaseous disk, which induces rapid inward migration for planets in the terrestrial mass range (Ward, 1997). Around solar-type stars, disk dispersal typically occurs on the order of a few million years (Walter et al., 1988; Strom et al., 1989). While there is evidence that disk lifetimes may exceed ∼5 Myr for low-mass stars (Carpenter et al., 2006; Pascucci et al., 2009), disks are no longer present after ∼10–20 Myr around all stellar types (Ribas et al., 2014). Planets migrating via interactions with the disk will thus settle into their new orbits relatively early. Planet-planet scattering, on the other hand, is by nature a stochastic process and may take place at any point during a system's lifetime. However, such interactions are far more frequent during the final stages of planet assembly (up to a few tens of millions of years), after which planets relax into their long-term quasi-stable orbits (Ford and Rasio, 2008).
Given the abundance of ices beyond the snow line, planets that migrate into the HZ from the outer regions of the disk should have abundant water, therefore satisfying that important criterion for habitability within the HZ. However, the situation is complicated by the fact that these planets may have accumulated thick gaseous envelopes, which could render them uninhabitable. Investigating whether these so-called “mini-Neptunes” can lose their envelopes and form planets with solid surfaces is therefore critical to understanding the habitability of planets around low-mass stars.
In this study, we focus on mini-Neptunes with initial masses in the range 1 MMp≤10 M with up to 50% hydrogen/helium by mass that have migrated into the HZs of mid to late M dwarf stars. We investigate whether it is possible for atmospheric escape processes to remove the thick H/He envelopes of mini-Neptunes in the HZ, effectively turning them into volatile-rich Earths and super-Earths (terrestrial planets more massive than Earth), which we refer to as “habitable evaporated cores” (HECs). We consider two atmospheric loss processes: XUV-driven escape, in which stellar X-ray/extreme ultraviolet (XUV) photons heat the atmosphere and drive a hydrodynamic wind away from the planet, and a simple model of Roche lobe overflow (RLO), in which the atmosphere is so extended that part of it lies exterior to the planet's Roche lobe; that gas is therefore no longer gravitationally bound to the planet. We further couple the effects of atmospheric mass loss to the thermal and tidal evolution of these planets. Planets cool as they age, undergoing changes of up to an order of magnitude in radius, which can greatly affect the mass loss rate. Tidal forces arising from the differential strength of gravity distort both the planet and the star away from sphericity, introducing torques that lead to the evolution of the orbital and spin parameters of both bodies.
While many studies have considered the separate effects of atmospheric escape (Erkaev et al., 2007; Murray-Clay et al., 2009; Tian, 2009; Owen and Jackson, 2012; Lammer et al., 2013), RLO (Trilling et al., 1998; Gu et al., 2003), thermal evolution (Lopez et al., 2012), and tidal evolution (Ferraz-Mello et al., 2008; Jackson et al., 2008; Correia and Laskar, 2011) on exoplanets, none have considered the coupling of these effects in the HZ. For some systems, in particular those that may harbor HECs, modeling the coupling of these processes is essential to accurately determine the evolution, since several feedbacks can ensue. Tidal forces in the HZ typically act to decrease a planet's semimajor axis, leading to higher stellar fluxes and faster mass loss. The mass loss, in turn, affects the rate of tidal evolution primarily via the changing planet radius, which is also governed by the cooling rate of the envelope. Changes to the star's radius and luminosity lead to further couplings that need to be treated with care.
Jackson et al. (2010) considered the effect of the coupling between mass loss and tidal evolution on the hot super-Earth CoRoT-7 b, finding that the two effects are strongly linked and must be considered simultaneously to obtain an accurate understanding of the planet's evolution. However, an analogous study has not been performed in the HZ, in great part because both tidal effects and atmospheric mass loss are generally orders of magnitude weaker at such distances from the star. This is not necessarily true around M dwarfs, for two reasons: (a) their low luminosities result in a HZ that is much closer in, exposing planets to strong tidal effects and possible RLO, and (b) M dwarfs are extremely active early on, so that the XUV flux in the HZ can be orders of magnitude higher than that around a solar-type star (see, for instance, Scalo et al., 2007).
In this paper, we present the results of the first model to couple tides, thermal evolution, and atmospheric mass loss in the HZ and show that for certain systems the coupling is key in determining the long-term evolution of the planet. We demonstrate that it is possible to turn mini-Neptunes into HECs within the HZ, providing an important pathway to the formation of potentially habitable, volatile-rich planets around M dwarfs.
In Section 2, we provide a detailed description of the relevant physics. In Section 3, we describe our model, followed by our results in Section 4. We then discuss our main findings and the corresponding caveats in Section 5, followed by a summary in Section 6. We present auxiliary derivations and calculations in the Appendix.

2. Stellar and Planetary Evolution

In the following sections, we review the luminosity evolution of low-mass stars (Section 2.1), the HZ and its evolution (Section 2.2), atmospheric escape processes from planets (Section 2.3), and tidal evolution of star-planet systems (Section 2.4).

2.1. Stellar luminosity

Following their formation from a giant molecular cloud, stars contract under their own gravity until they reach the main sequence, at which point the core temperature is high enough to ignite hydrogen fusion. While the Sun is thought to have spent ≲50 Myr in this pre-main sequence (PMS) phase (Baraffe et al., 1998), M dwarfs take hundreds of millions of years to fully contract; the lowest-mass M dwarfs reach the main sequence only after ∼1 Gyr (e.g., Reid and Hawley, 2005). During their contraction, M dwarfs remain at a roughly constant effective temperature (Hayashi, 1961), so that their luminosity is primarily a function of their surface area, which is significantly larger than when they arrive on the main sequence. M dwarfs therefore remain super-luminous for several hundred millions of years, with total (bolometric) luminosities higher than the main sequence value by up to 2 orders of magnitude. As we discuss below, this will significantly affect the atmospheric evolution of any planets these stars may host.
X-ray/extreme ultraviolet emissions (1 Å≲λ≲1000 Å) from M dwarfs also vary strongly with time. This is because the XUV luminosity of M dwarfs is rooted in the vigorous magnetic fields generated in their large convection zones (Scalo et al., 2007). Stellar magnetic fields are largely controlled by rotation (Parker, 1955), which slows down with time due to angular momentum loss (Skumanich, 1972), leading to an XUV activity that declines with stellar age. However, due to the difficulty of accurately determining both the XUV luminosities (usually inferred from line proxies) and the ages (often determined kinematically) of M dwarfs, the exact functional form of the evolution is still very uncertain (for a review, see Scalo et al., 2007). Further complications arise due to the fact that the processes that generate magnetic fields in late M dwarfs may be quite different from those in earlier-type stars. The rotational dynamo of Parker (1955) involves the amplification of toroidal fields generated at the radiative-convective boundary; late M dwarfs, however, are fully convective and have no such boundary. Instead, their magnetic fields may be formed by turbulent dynamos (Durney et al., 1993), which may evolve differently in time from those around higher-mass stars (Reid and Hawley, 2005).
In a comprehensive study of the XUV emissions of solar-type stars of different ages, Ribas et al. (2005) found that the XUV evolution is well modeled by a simple power law with an initial short-lived “saturation” phase:
\begin{align*} \frac { L_ { { \rm XUV } } } { L_ { { \rm bol } } } = \begin{cases} \left( \frac { L_ { { \rm XUV } } } { L_ { { \rm bol } } } \right) _ { { \rm sat } } \qquad \qquad\quad \ t \le t_ { { \rm sat } } \\ \left( \frac { L_ { { \rm XUV } } } { L_ { { \rm bol } } } \right) _ { { \rm sat } } \left( \frac { t } { t_ { { \rm sat } } } \right) ^ { - \beta } \quad\; \quad t > t_ { { \rm sat } } \end{cases} \tag { 1 } \end{align*}
where LXUV and Lbol are the XUV and bolometric luminosities, respectively, and β=−1.23. Prior to t=tsat, the XUV luminosity is said to be “saturated,” as observations show that the ratio LXUV/Lbol remains relatively constant at early times.
Jackson et al. (2012) found that tsat≈100 Myr and (LXUV/Lbol)sat≈10−3 for K dwarfs. Similar studies for M dwarfs, however, are still being developed (e.g., Engle and Guinan, 2011), but it is likely that the saturation timescale is much longer for late-type M dwarfs. Wright et al. (2011) showed that X-ray emission in low-mass stars is saturated for Prot/τ≲0.1, where Prot is the stellar rotation period and τ is the convective turnover time. The extent of the convective zone increases with decreasing stellar mass; below 0.35 M, M dwarfs are fully convective (Chabrier and Baraffe, 1997), resulting in larger values of τ (see, e.g., Pizzolato et al., 2000). Low-mass stars also have longer spin-down times (Stauffer et al., 1994); together, these effects should lead to significantly longer saturation times compared to solar-type stars. This is consistent with observational studies; West et al. (2008) found that the magnetic activity lifetime increases from ≲1 Gyr for early (i.e., most massive) M dwarfs to ≳7 Gyr for late (least massive) M dwarfs, possibly due to the onset of full convection around spectral type M3. Finally, Stelzer et al. (2013) found that for M dwarfs, β=−1.10±0.02 in the X-ray and β=−0.79±0.05 in the far ultraviolet, suggesting a slightly shallower slope in the XUV compared to the value from Ribas et al. (2005).

2.2. The habitable zone

The HZ is classically defined as the region around a star where an Earth-like planet can sustain liquid water on its surface (Hart, 1979; Kasting et al., 1993). Interior to the HZ, high surface temperatures lead to the runaway evaporation of a planet's oceans, which increases the atmospheric infrared opacity and reduces the ability of the surface to cool in a process known as the runaway greenhouse. Exterior to the HZ, greenhouse gases—gases, like water vapor, that absorb strongly in the infrared—are unable to maintain surface temperatures above the freezing point, and the oceans freeze globally. Recently, Kopparapu et al. (2013) recalculated the location of the HZ boundaries as a function of stellar luminosity and effective temperature using an updated one-dimensional, radiative-convective, cloud-free climate model. Their five boundaries are the (1) Recent Venus, (2) Runaway Greenhouse, (3) Moist Greenhouse, (4) Maximum Greenhouse, and (5) Early Mars HZs.
The first and last limits can be considered “optimistic” empirical limits, since prior to ∼1 and ∼3.8 Gyr ago, respectively, Venus and Mars may have had liquid surface water. The ability of a planet to maintain liquid water and to sustain life at these limits is still unclear and probably depends on a host of properties of its climate. Conversely, the other three limits are the “pessimistic” theoretical limits, corresponding to where a cloud-free, Earth-like planet would lose its entire water inventory due to the greenhouse effect (2 and 3) and to where the addition of CO2 to the atmosphere would be incapable of sustaining surface temperatures above freezing (4).
Because stellar luminosities vary with time, the location of the HZ is not fixed. In Fig. 1, we plot the HZ at 10 Myr (dashed lines) and at 1 Gyr (solid lines), calculated from the HZ model of Kopparapu et al. (2013) and the stellar evolution models of Baraffe et al. (1998). While for K and G dwarfs the change in the HZ is negligible, the HZ of M dwarfs moves in substantially, changing by nearly an order of magnitude for the least-massive stars. Due to this evolution, planets observed in the HZ of M dwarfs today likely spent a significant amount of time interior to the inner edge of the HZ, provided they either formed in situ or migrated to their current positions relatively early. Luger and Barnes (2014) explored in detail the effects of the evolution of the HZ on terrestrial planets.
FIG. 1. Location of the inner habitable zone (red), central habitable zone (green), and outer habitable zone (blue) as a function of stellar mass at 10 Myr (dashed) and 1 Gyr (solid). After 1 Gyr, the evolution of the HZ is negligible for M dwarfs.
Finally, we note that the location of the HZ is also a function of the eccentricity e. This is due to the fact that at a fixed semimajor axis a, the orbit-averaged flux 〈Fbol〉 is higher for more eccentric orbits (Williams and Pollard, 2002):
\begin{align*}\langle F_ { \rm bol } \rangle = \frac { L_ { { \rm bol } } } { 4 \pi a^2 \sqrt { 1 - e^2 } } \tag { 2 } \end{align*}

2.3. Atmospheric mass loss

Planetary atmospheres constantly evolve. Several mechanisms exist through which planets can lose significant quantities of their atmospheres to space, leading to dramatic changes in composition and in some cases complete atmospheric erosion. Early Earth itself could have been rich in hydrogen, with mixing ratios as high as 30% in the prebiotic atmosphere (Tian et al., 2005; Hashimoto et al., 2007). A variety of processes subsequently led to the loss of most of this hydrogen; Watson et al. (1981) argued that on the order of 1023 g of hydrogen could have escaped in the first billion years. Similarly, Kasting and Pollack (1983) calculated that early Venus could have lost an Earth ocean equivalent of water in the same amount of time. Currently, observational evidence for atmospheric escape exists for two “hot Jupiters,” HD 209458b (Vidal-Madjar et al., 2003) and HD 189733b (Lecavelier des Etangs et al., 2010), and one “hot Neptune,” GJ 436b (Kulow et al., 2014), whose proximity to their parent stars leads to the rapid hydrodynamic loss of hydrogen.
Atmospheric escape mechanisms fall into two major categories: thermal escape, in which the heating of the upper atmosphere accelerates the gas to velocities exceeding the escape velocity, and nonthermal escape, which encompasses a variety of mechanisms and may involve energy exchange among ions or erosion due to stellar winds. While nonthermal processes certainly do play a role in the evaporation of super-Earth and mini-Neptune atmospheres, the high exospheric temperatures resulting from strong XUV irradiation probably make thermal escape the dominant mechanism for planets around M dwarfs at early times. However, the escape can be greatly enhanced by flares and coronal mass ejections, which can completely erode the atmosphere of a planet lacking a strong magnetic field (Lammer et al., 2007; Scalo et al., 2007). For a review of the nonthermal mechanisms of escape, the reader is referred to Hunten (1982).

2.3.1. Jeans escape

In the low temperature limit, atmospheric mass loss occurs via the ballistic escape of individual atoms from the collisionless exosphere, where the low densities ensure that atoms with outward velocities exceeding the escape velocity will escape the planet. Originally developed by Jeans (1925), who assumed a hydrostatic, isothermal atmosphere, the mass loss rate of a pure hydrogen atmosphere is given by (Öpik, 1963)
\begin{align*} \frac { dM_p } { dt } = 4 \pi R_ { { \rm exo } } ^2 nm_H \upsilon_t \frac { ( 1 + \lambda_J ) e^ { - \lambda _ { J } } } { 2 \sqrt { \pi } } \tag { 3 } \end{align*}
where Rexo is the radius of the exobase, n is the number density of hydrogen atoms at the exobase, mH is the mass of a hydrogen atom, vt is the thermal velocity of the gas, and λJ is the Jeans escape parameter, defined as the ratio of the potential energy to the kinetic energy of the gas and given by
\begin{align*}\lambda_J \equiv \frac { GM_p m_H } { kT_ { { \rm exo } } R_ { { \rm exo } } } \tag { 4 } \end{align*}
where G is the gravitational constant, Mp is the mass of the planet, and Texo is the temperature in the (isothermal) exosphere. Since in the Jeans regime the thermal velocity of the gas is less than the escape velocity, the bulk of the gas remains bound to the planet, and only atoms in the tail of the Maxwell-Boltzmann distribution are able to escape. Jeans escape is thus slow. As an example, the present Jeans escape flux for hydrogen on Venus is on the order of 10 cm−2 s−1 (Lammer et al., 2006), corresponding to the feeble rate of ∼10−4 g/s, which would remove only one part in 1011 of the atmosphere per billion years.
For higher exospheric temperatures and/or larger values of Rexo, however, corresponding to low values of λJ, the atmosphere enters a regime in which the velocity of the average atom in the exosphere approaches the escape velocity of the planet. In this regime, the bulk of the upper atmosphere ceases to be hydrostatic (and isothermal), (4) is no longer valid, and the escape rate must be calculated from hydrodynamic models.

2.3.2. Hydrodynamic escape

One of the primary mechanisms for heating the exosphere and decreasing λJ is via strong XUV irradiation. XUV photons are absorbed close to the base of the thermosphere, where they deposit energy and heat the gas via the ionization of atomic hydrogen. This heating is balanced by the adiabatic expansion of the upper atmosphere, downward heat conduction, and cooling by recombination radiation. For sufficiently high XUV fluxes, the expansion of the atmosphere accelerates the gas to supersonic speeds, at which point a hydrodynamic wind is established similar to the solar Parker wind (Parker, 1965). Once the gas reaches the exosphere, it will escape the planet provided its kinetic energy exceeds the energy required to lift it out of the planet's gravitational well.
Since the kinetic energy of a hydrogen atom at the exobase is \( \frac { 3 } { 2 } kT_ { { \rm exo } } \), the condition λJ<1.5 implies that the kinetic energy of the gas is greater than the absolute value of its gravitational binding energy, and it should therefore begin to escape in bulk in a process commonly referred to as “blow-off.” Unlike in the Jeans escape regime, where the mass loss occurs on a per-particle basis, blow-off leads to the rapid loss of large portions of the upper atmosphere, irrespective of particle species, as atoms and molecules heavier than hydrogen are carried along by the hydrodynamic wind. However, contrary to what Öpik (1963) suggested, the mass loss in this stage is not arbitrarily high, since once blow-off begins, the upper atmosphere can no longer be treated as isothermal. As the exosphere expands, it also cools, so that in the absence of an energy source the value of λJ will tend to increase, thereby moderating the blow-off. The mass loss is, in this sense, “energy-limited” and may be calculated by equating the energy input to the energy required to drive the escape.
Originally proposed by Watson et al. (1981), the energy-limited model assumes that the XUV flux is absorbed in a thin layer at radius RXUV where the optical depth to stellar XUV photons is unity. Recently updated to include tidal effects by Erkaev et al. (2007), this model approximates the mass loss as
\begin{align*} \frac { dM_p } { dt } \approx \frac { \epsilon_ { { \rm XUV } } \pi F_ { { \rm XUV } } R_p R_ { { \rm XUV } } ^2 } { GM_p K_ { { \rm tide } } ( \xi ) } \tag { 5 } \end{align*}
where εXUV is the heating efficiency parameter (see below), FXUV is the incident XUV flux, Rp is the planet radius, and Ktide is a tidal enhancement factor, accounting for the fact that for sufficiently close-in planets, the stellar gravity reduces the gravitational binding energy of the gas such that it need only reach the Roche radius to escape the planet. Erkaev et al. (2007) showed that
\begin{align*}K_ { { \rm tide } } ( \xi ) = \left( 1 - \frac { 3 } { 2 \xi } + \frac { 1 } { 2 \xi^3 } \right) \tag { 6 } \end{align*}
where the parameter ξ is defined as
\begin{align*}\xi \equiv \frac { R_ { { \rm Roche } } } { R_ { { \rm XUV } } } \tag { 7 } \end{align*}
with
\begin{align*}R_ { { \rm Roche } } \equiv \left( \frac { M_p } { 3 M_\star } \right) ^ { 1 / 3 } a \tag { 8 } \end{align*}
where M is the mass of the star and a is the semimajor axis. For simplicity, as in the work of Lopez et al. (2012), we replace Rp with RXUV in (5), which is approximately valid given that RXUV is typically only 10–20% larger than Rp (Murray-Clay et al., 2009; Lopez et al., 2012).
Since the input XUV energy is balanced in part by cooling radiation (via Lyman α emission in the case of hydrogen) and heat conduction, only a fraction of it goes into the adiabatic expansion that drives escape. Rather than running complex hydrodynamic and radiative transfer models to determine the net heating rate, many authors (Jackson et al., 2010; Leitzinger et al., 2011; Lopez et al., 2012; Koskinen et al., 2013; Lammer et al., 2013) choose to fold the balance between XUV heating and cooling into an efficiency parameter, εXUV, defined as the fraction of the incoming XUV energy that is converted into PdV work. Because of the complex dependence of εXUV on the atmospheric composition and structure, its value is still poorly constrained. Lopez et al. (2012) estimated εXUV=0.2±0.1 for super-Earths and mini-Neptunes based on values found throughout the literature. Earlier work by Chassefière (1996) estimates εXUV≲0.30 for Venus-like planets, but the author argues that the actual value may be closer to 0.15. Recently, Owen and Jackson (2012) found X-ray efficiencies ≲0.1 for planets more massive than Neptune, but higher efficiencies (∼0.15) for terrestrial planets. Moreover, Shematovich et al. (2014) argued that studies that assume efficiencies higher than about 0.2 probably lead to overestimates in the escape rate. On the other hand, some studies suggest higher heating efficiencies: Koskinen et al. (2013) used hydrodynamic and photochemical models of the hot Jupiter HD 209458b to calculate εXUV=0.44.
As we have already implied, unlike Jeans escape, hydrodynamic blow-off is fast. Chassefière (1996) calculated the maximum hydrodynamic escape rate from the early venusian atmosphere to be ∼106 g/s, 10 orders of magnitude higher than the present Jeans escape flux (see Section 2.3.1). Although there has been debate over the validity of the blow-off assumption (see, for instance, the discussion in Tian et al., 2008), recently Lammer et al. (2013) showed that, for super-Earths exposed to high levels of XUV irradiation, the energy-limited approximation yields mass loss rates that are consistent with hydrodynamic models to within a factor of about 2, which is within the uncertainties of the problem.

2.3.3. Controlled hydrodynamic escape

It is also worth noting that there may be an intermediate regime between Jeans escape and blow-off known as “modified Jeans escape” or “controlled hydrodynamic escape” (Erkaev et al., 2013). In this regime, which occurs for intermediate XUV fluxes and/or higher planetary surface gravity, blow-off conditions are not met, but the atmosphere still expands, so that the hydrostatic Jeans formalism is not valid. To calculate the escape rate, one must replace the classical Maxwellian velocity distribution with one that includes the bulk expansion velocity of the atmosphere. This yields escape rates lower than those due to a hydrodynamic flow but significantly higher than those predicted by the hydrostatic Jeans equation (3).

2.3.4. Jeans escape or hydrodynamic escape

Since the location of the HZ is governed primarily by the total (bolometric) flux incident on a planet, the higher ratio of LXUV to Lbol of M dwarfs implies a much larger XUV flux in the HZ compared to solar-type stars. The present-day solar XUV luminosity is LXUV/Lbol≈3.4×10−6 (see Table 4 in Ribas et al., 2005), while for active M dwarfs this ratio is ∼10−3 (e.g., Scalo et al., 2007). Therefore, we should expect planets in the HZ of M dwarfs to experience XUV fluxes several orders of magnitude greater than the present Earth level (FXUV≈4.64 erg/s/cm2). Recent papers (Lammer et al., 2007, 2013; Erkaev et al., 2013) show that terrestrial planets experiencing XUV fluxes corresponding to 10× and 100×FXUV are in the hydrodynamic flow regime, and we may thus expect the same for super-Earths/mini-Neptunes in the HZ of active M dwarfs.
In Fig. 2, we plot the evolution of the XUV flux received by a planet located close to the inner edge of the HZ (defined at 5 Gyr) for three different M dwarf masses and two different XUV saturation times (see Section 2.1). The dashed lines correspond to the critical fluxes in the work of Erkaev et al. (2013) above which hydrodynamic escape occurs, for 1 and 10 M and two values of εXUV. Earth-mass planets remain in the hydrodynamic escape regime for at least 1 Gyr in all cases and in excess of 10 Gyr for active M dwarfs. The duration of this regime is shorter for 10 M planets but still on the order of several 100 Myr.
FIG. 2. Evolution of the XUV flux received by planets close to the inner edge of the HZ (at 1 Gyr) for stars of mass 0.1, 0.2, and 0.3 M. Solid lines correspond to an XUV saturation time of 0.1 Gyr; dashed lines correspond to 1 Gyr. The flux at Earth is indicated by the black line. The dot corresponds to the earliest time for which Ribas et al. (2005) had data for solar-type stars; this is also roughly the time at which Earth formed. An XUV luminosity saturated at 10−3 Lbol is roughly indicated by the dotted line. Finally, the dashed gray lines indicate the minimum XUV fluxes required to sustain blow-off according to the study of Erkaev et al. (2013). Super-Earths in the HZs of M dwarfs remain in the blow-off regime for at least a few 100 Myr; Earths undergo blow-off for much longer. For an XUV saturation time of 1 Gyr, blow-off occurs for several billion years for all planets.
We note also that tidal effects can significantly increase the critical value of λJ below which hydrodynamic escape occurs (see discussion in Erkaev et al., 2007). Hydrodynamic escape ensues when the thermal energy of the gas exceeds its potential energy, which occurs when
\begin{align*}1.5 \ge \frac { GM_p K_ { { \rm tide } } m_H } { kT_ { { \rm exo } } R_ { { \rm exo } } } \end{align*}
or
\begin{align*} \begin{split} \lambda_J \le \frac { 1.5 } { K_ { { \rm tide } } } \\ \equiv \lambda_ { { \rm crit } } \end{split} \tag { 9 } \end{align*}
provided we maintain the original definition of the Jeans parameter (4). Due to the strong tidal forces acting on the planets we consider here, this effect should greatly increase the critical value of the escape parameter, effectively reducing the value of FXUV for which hydrodynamic escape occurs.

2.3.5. Energy-limited or radiation/recombination-limited?

Hydrodynamic escape from planetary atmospheres need not be energy-limited. In the limit of high extreme ultraviolet (EUV) flux (low-energy XUV photons with 100 Å≲λ≲1000 Å), Murray-Clay et al. (2009) showed that the escape is “radiation/recombination-limited,” scaling roughly as \(\dot {M}\)∝(FEUV)1/2. In this regime, the upper atmosphere thermostats to T∼104 K, photoionization is balanced by radiative recombination (as opposed to PdV work), and a large fraction of the gas energy budget is lost to Lyman α cooling. The mass loss rate is found by solving the mass continuity equation, yielding
\begin{align*}\dot{ M} = 4 \pi r_s^2 \rho_s c_s \tag{10}\end{align*}
where rs is the radius of the sonic point (where the wind velocity equals the local sound speed cs) and ρs is the density at rs. Since the bulk of the flow is ionized, the density is fixed by ionization-recombination balance, scaling roughly as (FEUVrs)1/2. For a 0.7 MJ hot Jupiter, the radiation/recombination-limited mass loss rate is (Murray-Clay et al., 2009)
\begin{align*}\dot { M } _ { { \rm RR } } \approx 4 \times 10^ { 12 } { \rm g / s } \left( \frac { F_ { { \rm EUV } } } { 5 \times 10^5 \; { \rm erg / cm^2 / s } } \right) ^ { 1 / 2 } \tag { 11 } \end{align*}
Owen and Jackson (2012) re-derived this expression with explicit scalings on several planetary properties:
\begin{align*} \dot {M} _ {{\rm RR}} \approx \,2.4 \times 10^ {11} {\rm g / s} \ ( 1 + x ) ^ {3 / 2} \ f_ {{\rm parker}} \left( \frac{\Phi_\star} {{10^40} \ {\rm s}^{- 1}} \right)^{1 / 2} \\ \times \left( \frac {a} {0.1 \ {\rm AU}} \right) \left( \frac {R_p} {10 \ {R} _ {\oplus}} \right) ^ {3 / 2} \left( \frac {A} {1 / 3} \right) ^ {1 / 2} \left( \frac {c_ {{\rm EUV}}} {10 \ {\rm km / s}} \right) \tag {12}\end{align*}
where the quantity (1+x) is the radius of the ionization front in units of Rp, fparker is the Mach number of the flow, Φ is the stellar EUV luminosity in photons per second, A is a geometrical factor, and cEUV is the isothermal sound speed of the gas. Taking x=0, fparker=1, A=1/3, cEUV=10 km/s, and an average EUV photon energy =20 eV, this becomes
\begin{align*} \dot { M } _ { { \rm RR } } \approx 7.11 \times 10^7 { \rm g / s } \, \left( \frac {F_{\rm EUV } } {{\rm erg } / {\rm cm}^2 / {\rm s }} \right) ^ { 1 / 2 } \left( \frac { R_p } { R_ { \oplus } } \right) ^ { 3 / 2 } \tag { 13 } \end{align*}
The transition from energy-limited to radiation/recombination-limited escape is found by equating the two expressions, namely, (5) and (13), and solving for the critical EUV flux. For hot Jupiters and mini-Neptunes alike, the transition occurs at roughly FEUV∼104 erg/cm2/s. Below this flux, the escape is energy-limited; above it, the escape is radiation/recombination-limited and thus increases more slowly with the flux.
Mini-Neptunes that migrate early into the HZ of M dwarfs are exposed to EUV fluxes up to an order of magnitude larger than this critical value. During this period, which lasts on the order of a few hundred million years, the mass loss rate may be radiation/recombination-limited.
However, whether the high-flux mass loss rate is more accurately described by (5) or (13) will depend on whether the flux is dominated by X-ray or EUV radiation. Owen and Jackson (2012) showed that the mass loss rate for X-ray-driven hydrodynamic winds scales linearly with the X-ray flux; this is because the sonic point for X-ray flows tends to occur below the ionization front. Even though recombination radiation still removes energy from the flow, it does so once the gas is already supersonic and thus causally decoupled from the planet, such that it cannot bottleneck the escape. Although the authors caution that the dependence of the mass loss rate on planet mass and radius must be determined numerically, this regime is analogous to the energy-limited regime and can be roughly approximated by (5).
For X-ray luminosities LX≳10−6 L, close-in planets undergo X-ray driven hydrodynamic escape (see Fig. 11 in Owen and Jackson, 2012). If X-rays contribute significantly to the XUV emissions of young M dwarfs, their X-ray luminosities may exceed this value for as long as 1 Gyr, and close-in mini-Neptunes will undergo energy-limited escape. Unfortunately, given the lack of observational constraints on the X-ray/EUV luminosities of young M dwarfs, it is unclear at this point whether the hydrodynamic escape will be EUV-driven (radiation/recombination-limited) or X-ray-driven (energy-limited).

2.3.6. The effect of eccentricity

Most of the formalism that has been developed to analytically treat hydrodynamic blow-off only considers circular orbits. For planets on sufficiently eccentric orbits, neither the stellar flux nor the tidal effects may be treated as constant over the course of an orbit. Due to this fact, an expression analogous to (5) for eccentric orbits seems to be lacking in the literature. In this section, we derive such an expression.
There are two separate effects that enhance the mass loss for planets on eccentric orbits. Most papers account for the first effect, which is the increase of the orbit-averaged stellar flux by a factor of 1/\(\sqrt{1 - e^2}\) (see, for instance, Kopparapu et al., 2013). However, for e≲0.3, this effect is quite small. The second, more important effect is that the Roche lobe radius is no longer constant over the course of an orbit, and (8) is not valid. Instead, we must replace a with the instantaneous planet-star separation r(t):
\begin{align*}R_ { { \rm Roche } } ( t ) = \left( \frac { M_p } { 3M_\star } \right) ^ { 1 / 3 } r ( t ) \tag { 14 } \end{align*}
One might wonder whether this replacement is valid. Specifically, if RRoche(t) changes faster than the atmosphere is able to respond to the changes in the gravitational potential, we would expect that the time dependence of the mass loss rate would be a complicated function of the tides generated in the atmosphere. On the other hand, if the orbital period is very large compared to the dynamical timescale of the planet, the atmosphere will have sufficient time to assume the equilibrium shape dictated by the new potential. This limit is known as the quasi-static approximation (Sepinsky et al., 2007). In the Appendix, we show that all the planets in our runs with eccentricities e≲0.4 are in the quasi-static regime and that (14) is therefore valid. In some runs, we allow the eccentricity to increase beyond 0.4. We discuss the implications of this in Section 5.9.
Since RRoche=RRoche(t), ξ, Ktide, and dM/dt (Eqs. 5, 6, and 7) are now also functions of t, varying significantly over a single orbit. To account for this, we may calculate the time-averaged mass loss rate over the course of one orbit, \(\langle \,\dot{ M} \rangle_t\), such that the total amount of mass lost in time Δt is \(\Delta M = \langle\, \dot{ M} \rangle_t \Delta t\). To this end, in the Appendix we derive the eccentric version of Ktide:
\begin{align*} \frac { 1 } { K_ { { \rm ecc } } } \equiv \frac { 1 } { 2 \pi } \int_0^ { 2 \pi } \left[ ( 1 - e \cos E ) - \frac { 3 } { 2 \xi } + \frac { 1 } { 2 \xi^3 ( 1 - e \cos E ) ^2 } \right] ^ { - 1 } dE \tag { 15 } \end{align*}
where ξ is the time-independent parameter given by (7) and (8), and E is the eccentric anomaly. The average mass loss rate is then simply
\begin{align*}\langle\, \dot { M } \rangle_t = \frac { \dot { M } _0 } { K_ { { \rm ecc } } } \tag { 16 } \end{align*}
where
\begin{align*}\dot { M } _0 \equiv \frac { R_ { { \rm XUV } } ^3 \epsilon_ { { \rm XUV } } L_ { { \rm XUV } } } { 4 GM_p a^2 } \tag { 17 } \end{align*}
is the zero-eccentricity mass loss rate in the absence of tidal enhancement. Note that the flux-enhancement factor \(1 / \sqrt{1 - e^2}\) is already folded into Kecc, since it must be incorporated when integrating (33).
For ξ≳10, the integral may be approximated by the analytic expression
\begin{align*}K_ { { \rm ecc } } \approx \sqrt { 1 - \frac { 3 } { \xi_0 } - \frac { 9 } { 4 \xi_0^2 } - e^2 } \tag { 18 } \end{align*}
which greatly reduces computing time.
As we show in the Appendix, the decreased Roche lobe distance for eccentric orbits has a large effect on the amount of mass lost, particularly for low values of ξ and for high e (see Fig. A2). Moreover, higher eccentricities result in RLO at larger values of a compared to the circular case, since the planet may overflow near pericenter, leading to mass loss rates potentially orders of magnitude higher. We discuss RLO in detail in the Model Description section (Section 3.3).

2.4. Tidal theory

The final aspect of planetary evolution we discuss is the effect of tidal interactions with the host star. Below we review two different approaches to analytically calculate the orbital evolution of the system.

2.4.1. Constant phase lag

Classical tidal theory predicts that torques arising from interactions between tidal deformations on a planet and its host star lead to the secular evolution of the orbit and the spin of both bodies. In this paper, we employ the “equilibrium tide” model of Darwin (1880), which approximates the tidal potential as a superposition of Legendre polynomials representing waves propagating along the surfaces of the bodies; these add up to give the familiar tidal “bulge.” Because of viscous forces in the bodies' interiors, the tidal bulges do not instantaneously align with the line connecting the two bodies. Instead, the Nth wave on the ith body will lag or lead by an angle ɛN,i, assumed to be constant in the constant phase lag (CPL) model. In general, different waves may have different lag angles, and it is unclear how the ɛN,i vary as a function of frequency. A common approach (see Ferraz-Mello et al., 2008) is to assume that the magnitudes of the lag angles are equal (see Goldreich and Soter, 1966), while their signs may change depending on the orbital and rotational frequencies involved. This allows us to introduce the tidal quality factor
\begin{align*}Q_i \equiv \frac { 1 } { \varepsilon_ { 0 , i } } \tag { 19 } \end{align*}
which in turn allows us to express the lags (in radians) as
\begin{align*}\varepsilon_ { N , i } = \pm \frac { 1 } { Q_i } \tag { 20 } \end{align*}
The parameter Qi is a measure of the dissipation within the ith body; it is inversely proportional to the amount of orbital and rotational energy lost to heat per cycle, in analogy with a damped-driven harmonic oscillator. The merit of this approach is that the tidal response of a body can be captured in a single parameter. Planets with high values of Qp have smaller phase lags, dissipate less energy, and undergo slower orbital evolution; planets with low values of Qp have larger phase lags, higher dissipation rates, and therefore faster evolution. Measurements in the Solar System constrain the value of Qp for terrestrial bodies in the range 10–500, while gas giants are consistent with Qp∼104 to 105 (Goldreich and Soter, 1966). Values of Q for the Sun and other main sequence stars are poorly constrained but are likely to be ≳105 to 106 (Schlaufman et al., 2010; Penev et al., 2012). Intuitively, this makes sense, given that the dissipation due to internal friction in rocky bodies should be much higher than that in bodies dominated by gaseous atmospheres. One should bear in mind, however, that the exact dependence of Qi on the properties of a body's interior is likely to be extremely complicated. Given the dearth of data on the composition and internal structure of exoplanets, it is at this point impossible to infer precise values of Qp for these planets.
By calculating the forces and torques due to the tides raised on both the planet and the star, one can arrive at the secular expressions for the evolution of the planet's orbital parameters, which are given by a set of coupled nonlinear differential equations; these are reproduced in the Appendix.

2.4.2. Constant time lag

Unlike the CPL model, which assumes the phase lag of the tidal bulge is constant, the constant time lag (CTL) model assumes that it is the time interval between the bulge and the passage of the perturbing body that is constant. Originally proposed by Alexander (1973) and updated by Leconte et al. (2010), this model allows for a continuum of tidal wave frequencies and therefore avoids unphysical discontinuities present in the CPL model. However, implicit in the CTL theory is the assumption that the lag angles are directly proportional to the driving frequency (Greenberg, 2009), which is likely also an oversimplification. We note, however, that in the low eccentricity limit, both the CPL and the CTL models arrive at qualitatively similar results. At higher eccentricities, the CTL model is probably better suited, given that it is derived to eighth order in e (versus second order in the CPL model).
The tidal quality factors Qi do not enter the CTL calculations at any point; instead, the dissipation is characterized by the time lags τi. Although there is no general conversion between Qi and τi, Leconte et al. (2010) showed that provided annual tides dominate the evolution,
\begin{align*}\tau_i \approx \frac { 1 } { n Q_i } \tag { 21 } \end{align*}
where n is the mean motion (or the orbital frequency) of the secondary body (in this case, the planet).
For a planet with Qp=104 in the center of the HZ of a late M dwarf, τp≈10 s; rocky planets with lower Qp may have values on the order of hundreds of seconds. Since τn−1, close-in planets should have much lower time lags. For reference, Leconte et al. (2010) argued that hot Jupiters should have 2×10−3 s≲τp≲2×10−2 s.
The tidal evolution expressions are reproduced in the Appendix. For a more detailed review of tidal theory, the reader is referred to Ferraz-Mello et al. (2008), Heller et al. (2011), and the appendices in Barnes et al. (2013).

3. Model Description

Our model evolves planet-star systems forward in time in order to determine whether HECs can form from mini-Neptunes that have migrated into the HZs of M dwarfs. We perform our calculations on a grid of varying planetary, orbital, and stellar properties to determine the types of systems that may harbor HECs. The complete list is provided in Table 1, where we indicate the ranges of values we consider as well as the default values adopted in the plots in Section 4 (unless otherwise indicated).
Table 1. Free Parameters and Their Ranges
ParameterRangeDefaultNotes
M (M)0.08–0.5Late-mid M dwarfs
Mp (M)1–10
RXUV (Rp)1.0–1.21.2See Section 2.3.2
aIHZ–OHZSee Section 3.1
e0.0–0.95
P0,★ (days)1.0–10030.0Initial rotation period
fH10−6 to 0.5H mass fraction
εXUV0.1–0.40.3
ξmin1+10−5 to 33See Section 3.3
Atmospheric escapeR/R-Lim/E-LimSee Section 3.3
Tidal modelCPL/CTLCTL
Q105 to 106105CPL only
Qp101 to 105104CPL only
τ (s)10−2 to 10−110−1CTL only
τp (s)10−3 to 10310−1CTL only
β0.7–1.231.23See Eq. 1
tsat (Gyr)0.1–1.01.0XUV saturation time
t0 (Myr)10.0–100.010.0Integration start
tstop (Gyr)0.01–5.05.0Integration end
Integrations are performed from t=t0 (the time at which the planet is assumed to have migrated into the HZ) to t=tstop (the current age of the system) with the adaptive timestepping method described in Appendix E of Barnes et al. (2013).

3.1. Stellar model

We use the evolutionary tracks of Baraffe et al. (1998) for solar metallicity to calculate Lbol and Teff as a function of time. We then use (1) to calculate LXUV, given (LXUV/Lbol)sat=10−3 and values of tsat and β given in Table 1.
Using Lbol and Teff, we calculate the location of the HZ from the equations given by Kopparapu et al. (2013), adding the eccentricity correction (2). Given the uncertainty in the actual HZ boundaries and their dependence on a host of properties of a planet's climate, we choose our inner edge of the habitable zone (IHZ) to be the average of the Recent Venus and the Runaway Greenhouse limits and our outer edge of the habitable zone (OHZ) to be the average of the Maximum Greenhouse and the Early Mars limits. Throughout this paper, we will also refer to the center of the habitable zone (CHZ), which we take to be the average of the IHZ and OHZ. Since we are concerned with the formation of ultimately habitable planets, we take the locations of the IHZ, CHZ, and OHZ to be their values at 1 Gyr, at which point the stellar luminosity becomes roughly constant.

3.2. Planet radius model

To determine the planet radius Rp as a function of the core mass Mc, the envelope mass fraction fHMe/Mp, and the planet age, we use the planet structure model described by Lopez et al. (2012) and Lopez and Fortney (2014), which is an extension of the model of Fortney et al. (2007) to low-mass low-density (LMLD) planets. These models perform full thermal evolution calculations of the interior as a function of time. In our runs, the core is taken to be Earth-like, with a mixture of two-thirds silicate rock and one-third iron, and the envelope is modeled as a H/He adiabat. A grid of values of Rp is then computed in the range 1 MMc≤10 M, 10−6fH≤0.5, and 107 years≤t≤1010 years. For values between grid points, we perform a simple trilinear interpolation. For gas-rich planets, Rp is the 20 mbar radius; for gas-free planets, it corresponds to the surface radius. The evolution of Rp with age due solely to thermal contraction is plotted in Fig. 3 for a few different planet masses and values of fH.
FIG. 3. Evolution of the radius as a function of time due to thermal contraction of the envelope, in the absence of tidal effects and atmospheric mass loss. From top to bottom, the plots correspond to planets with initial total masses (core+envelope) of 1, 2, and 5 M. Line styles correspond to different initial hydrogen mass fractions: 1% (red, solid), 10% (green, dashed), 25% (blue, dot-dashed), and 50% (black, dotted). For comparison, the gray shaded regions in the bottom two plots are the spread in radii calculated by Mordasini et al. (2012b) for fH≲0.20. See text for a discussion.
We note that the models of Fortney et al. (2007) and Lopez et al. (2012) are in general agreement with those of Mordasini et al. (2012a, 2012b) and, by extension, Rogers et al. (2011). Mordasini et al. (2012a) presented a validation of their model against that of Fortney et al. (2007), showing that for planets spanning 0.1 to 10 Jupiter masses, the two models predict the same radius to within a few percent. At the lower masses relevant to our study, the two models are also in agreement. To demonstrate this, in Fig. 3 we shade the regions corresponding to the spread in radii at a given mass and age in Fig. 9 of Mordasini et al. (2012b). Since those authors used a coupled formation/evolution code, at low planet mass the maximum envelope mass fraction fH is small; for a total mass of 4 M, Mordasini et al. (2012b) found that all planets have fH<0.2. At 2 M, most planets have fH≲0.1. We can see from Fig. 3 that at these values of fH, the two models predict very similar radius evolution. Note that Mordasini et al. (2012b) did not consider planets less massive than 2 M.
The maximum envelope fraction merits further discussion. Since we do not model the formation of mini-Neptunes, we do not place a priori constraints on the value of fH at a given mass; instead, we allow it to vary in the range 10−6fH≤0.5 for all planet masses. At masses ≲5 M, planets accumulate gas slowly and are typically unable to accrete more than ∼10–20% of their mass in H/He (Rogers et al., 2011; Bodenheimer and Lissauer, 2014); values of fH≈0.5 may thus be unphysical. However, as we argue in Section 5.1, the longer disk lifetimes around M dwarfs (Carpenter et al., 2006; Pascucci et al., 2009) allow more time for gas accretion, potentially increasing the maximum value of fH. Nevertheless, and more importantly, if a planet with fH=0.5 loses its entire envelope via atmospheric escape, any planet with the same core mass and fH<0.5 will also lose its envelope. Below, where we present integrations with fH=0.5, our results are therefore conservative, as planets with fH <<0.5 will in general evaporate more quickly.
While our treatment of the radius evolution is an improvement upon past tidal-atmospheric coupling papers (Jackson et al., 2010, for instance, calculated Rp for super-Earths by assuming a constant density as mass is lost), there are still issues with our approach, as follows: (1) We do not account for inflation of the radius due to high insolation. Instead, we calculate our radii from grids corresponding to a planet receiving the same flux as Earth. While at late times this is justified, since planets in the HZ by definition receive fluxes similar to Earth, at early times this is probably a poor approximation; recall that planets in the HZ around low-mass M dwarfs are exposed to fluxes up to 2 orders of magnitude higher during the host star's PMS phase. The primary effect of a higher insolation is to act as a blanket, delaying the planet's cooling and causing it to maintain an inflated radius for longer. This will result in mass loss rates higher than what we calculate here. (2) Since we are determining the radii from pre-computed grids, we also do not model the effect of tidal dissipation on the thermal evolution of the planet. Planets undergoing fast tidal evolution can dissipate large amounts of energy in their interiors, which should lead to significant heating and inflation of their radii. (3) The radius is also likely to depend on the mass loss rate. Setting Rp equal to the tabulated value for a given mass, age, and composition is valid only as long as the timescale on which the planet is able to cool is significantly shorter than the mass loss timescale. Otherwise, the radius will not have enough time to adjust to the rapid loss of mass, and the planet will remain somewhat inflated, leading to a regime of runaway mass loss (Lopez et al., 2012). While the planets considered here are probably not in the runaway regime [Lopez et al. (2012) found that runaway mass loss occurred only for H/He mass fractions ≳90%], we might still be significantly underestimating the radii during the early active phase of the parent star.
All points outlined above lead to an underestimate of the radius at a given time. Since the mass loss rate is proportional to \(R_{ \rm p}^3\) (5) or \(R_{ \rm p}^{3 / 2}\) (13), calculating the radius in this fashion leads to a lower bound on the amount of mass lost and on the strength of the coupling to tidal effects. Because our present goal is to determine whether it is possible to form HECs via this mechanism, this conservative approach is sufficient. Future work will incorporate a self-consistent thermal structure model to better address the radius evolution.

3.3. Atmospheric escape model

We assume that the escape of H/He from the planetary atmosphere is hydrodynamic (blow-off) at all times, which is valid at the XUV fluxes we consider here (see Erkaev et al., 2013, and Fig. 2). We run two separate sets of integrations: one in which we assume the flow is energy-limited (5) for all values of FXUV, and one in which we switch from energy-limited to radiation/recombination-limited (13) above the critical value of the flux (see Section 2.3.5). For planets whose orbits are eccentric enough that they switch between the two regimes over the course of one orbit, we make use of the expressions derived in Section A.3 in the Appendix. These two sets of integrations should roughly bracket the true mass loss rate.
For eccentric orbits, we calculate the mass loss in the energy-limited regime from (16), with Kecc evaluated from (15). We vary εXUV and RXUV in the ranges given in Table 1. We choose εXUV=0.30 as our default case. While this is consistent with values cited in the literature (see Section 2.3.2), it could be an overestimate. We discuss the implications of this choice in Section 5.3.
Given the large planetary radii at early times, many of the planets we model here are not stable against RLO in the HZ. During RLO, the stellar gravity causes the upper layers of the atmosphere to suddenly become unbound from the planet; this occurs when Rp>RRoche, where RRoche is given by (8). For a planet that forms and evolves in situ, RLO never occurs, since any gas that would be lifted from the planet in this fashion would have never accreted in the first place. However, an inflated gaseous planet that forms at a large distance from the star may initially be stable against overflow and enter RLO as it migrates inward (since RRochea). This is particularly the case for planets in the HZs of M dwarfs, since a and consequently RRoche are small.
Ideally, the tidally enhanced mass loss rate equation (5) should capture this process, but instead it predicts an infinite mass loss rate as RXUVRRoche (or as ξ→1) and unphysically changes sign for ξ<1. This is due to the fact that the energy-limited model implicitly assumes that the bulk of the atmosphere is located at RXUV (the single-layer assumption). Realistically, we would expect the planet to quickly lose any mass above the Roche lobe and then return to the stable hydrodynamic escape regime. However, upon loss of the material above RRoche, the portion of the envelope below the new XUV absorption radius R′XUV will not be in hydrostatic equilibrium; instead, an outward flow will tend to redistribute mass to the evacuated region above, leading to further overflow.
Several models exist that allow one to calculate the mass loss rate due to RLO (e.g., Ritter, 1988; Trilling et al., 1998; Gu et al., 2003; Sepinsky et al., 2007). These often involve calculating the angular momentum exchange between the outflowing gas and the planet, which can lead to its outward migration, given by
\begin{align*} \frac { 1 } { a } \frac { da } { dt } = - \frac { 2 } { M_p } \frac { dM_p } { dt } \tag { 22 } \end{align*}
for a planet on a circular orbit (Gu et al., 2003; Chang et al., 2010). This leads to a corresponding increase in RRoche until it reaches RXUV and the overflow is halted. By differentiating the stability criterion RXUV(Mp)=RRoche(Mp), one may then obtain an approximate expression for dMp/dt in terms of the density profile dM(<R)/dR of the envelope.
However, for mini-Neptunes that migrate into the HZ early on, RLO should occur during the initial migration process, which we do not model in this paper. Instead, we begin our calculations by assuming that our planets are stable to RLO in the HZ. If a planet's radius initially exceeds the Roche lobe radius, we set its envelope mass equal to the maximum envelope mass for which it can be stable at its current orbit; the difference between the two envelope masses is the amount of H/He it must have lost prior to its arrival in the HZ. It is important to note that these planets will initially have RXUV=RRoche, which, as we mentioned above, leads to an infinite mass loss rate in (5). An accurate determination of \(\dot{M}_{\rm p}\) in this case probably requires hydrodynamic simulations. However, the mass loss rate can be approximated by imposing a minimum value ξmin in (5). For ξ<ξmin, we set the mass loss rate equal to \(\dot{M}_p ( \xi = \xi_{ \min} )\). This is equivalent to imposing a maximum mass loss enhancement factor 1/Ktide, preventing the mass loss rate from reaching unphysical values as RXUVRRoche.
In Fig. 4 we show how the time tevap at which complete evaporation takes place scales with ξmin for a typical mini-Neptune in the IHZ (red), CHZ (green), and OHZ (blue). Also plotted is the mass loss enhancement factor 1/Ktide (Eq. 6, black dashed line) as a function of ξ=ξmin. Interestingly, despite the fact that the instantaneous mass loss rate (\(\dot{M}_{\rm p} \propto 1 / K_{{ \rm tide}}\)) approaches infinity as ξ→1, the evaporation time is relatively constant for ξmin≲3. This is because for very large \(\dot{M}_{\rm p}\) the planet loses sufficient mass in an amount of time Δt to decrease Rp substantially and terminate the overflow. In other words, cases in which ξ≈1 are very unstable, and as mass is lost ξ will quickly increase beyond ∼3. Both the net amount of mass lost and the evaporation time are therefore insensitive to the particular value of ξmin, provided it is less than about 3. We therefore choose ξmin=3 as the default value for our runs, noting that this corresponds to a maximum mass loss enhancement of 1/Ktide≈2.
FIG. 4. Complete evaporation time tevap as a function of the cutoff value ξmin for a 2 M mini-Neptune with initial fH=0.5 on a circular orbit around a 0.08 M star. The red, green, and blue lines correspond to planets in the IHZ, CHZ, and OHZ, respectively. Also plotted is the mass loss enhancement factor 1/Ktide (black dashed line), which approaches infinity as ξ→1. Note that for ξmin≲3, the evaporation time is relatively insensitive to the exact cutoff value, despite the fact that 1/Ktide blows up. We therefore choose ξmin=3 as the default cutoff, corresponding to a maximum enhancement factor 1/Ktide≈2.

3.4. Tidal model

We calculate the evolution of the semimajor axis, the eccentricity, and the rotation rates from (60)–(62) and (73)–(75) in the Appendix for the CPL and CTL models, respectively. For simplicity, we set the obliquities of all our planets to zero. Since the tidal locking timescale is very short for close-in planets1 (Gu et al., 2003), we assume that the planet's rotation rate is given by the equilibrium value (65) or (78).
We calculate the stellar spin by assuming different initial periods (see Table 1) and evolving it according to the tidal equations, while enforcing conservation of angular momentum as the star contracts during the PMS phase. We neglect the effects of rotational braking (Skumanich, 1972), whereby stars lose angular momentum to winds and spin down over time. This leads to an overestimate of the spin rate at later times, but tidal effects should only be important early on, when the radii and the eccentricity are higher. Bolmont et al. (2012) recently modeled the coupling between stellar spin and tides, following the evolution of a “slow rotator” star (P0≈10 days) and a “fast rotator” star (P0≈1 day). In both cases, the stars sped up during the first ∼300–500 Myr, after which time angular momentum loss became significant. However, tidal evolution is orders of magnitude weaker at such late times, so we would expect rotational braking to have a minimal effect on the tidal evolution. Moreover, as we show in the Appendix, in general it is the tide raised on the planet that dominates the evolution; as this depends on the planetary rotation rate, and not on the stellar rotation rate, our results are relatively insensitive to the details of the spin evolution of the star.
In the CPL model, we adopt typical gas giant values 104Qp≤105 for gas-rich mini-Neptunes and typical terrestrial values 10≤Qp≤500 for planets that have lost their envelopes; we assume stellar values in the range 105Q≤106. In the CTL model, we consider time lags in the range 10−3 s≤τp≤101 s for gas-rich mini-Neptunes and 10−1 s≤τp≤103 s once they lose their envelopes. Following Leconte et al. (2010), we consider stellar time lags in the range 10−2 s≤τ≤10−1 s.
Once a mini-Neptune loses all its atmosphere, we artificially switch Qp or τp to the terrestrial value adopted in that run. In reality, as the atmosphere is lost, the transition from high to low Qp (or low to high τp) should be continuous. A detailed treatment of the dependence of Qp and τp on the envelope mass fraction is deferred to future work.
Finally, we note that the second-order CPL model described above is valid only at low eccentricity. For \(e \ge \sqrt{1 / 19} \approx 0.23\), the phase lag of the dominant tidal wave discontinuously changes from negative to positive, such that the model then predicts outward migration due to the planetary tide. This effect is unphysical, stemming from the fact that the CPL model considers only terms up to second order in the eccentricity (for a detailed discussion of this, see Leconte et al., 2010). We therefore restrict all our calculations in the CPL framework to values of the eccentricity e≤0.2. For higher values of e, we use the higher-order CTL model.

4. Results

4.1. A typical run

In Fig. 5, we show the time evolution of three mini-Neptunes as a guide to understanding the results presented in the following sections. We plot planet mass and radius (top row), semimajor axis and eccentricity (center row), and stellar XUV luminosity and radius (bottom row), all as a function of time since the planet's initial migration, t0, for 5 Gyr. In the center plots, the post-1 Gyr OHZ and CHZ are shaded blue and green, respectively. At t=t0, the planet is a 1 M core with a 1 M envelope orbiting around a 0.08 M M dwarf with e=0.5 at a semimajor axis of 0.04 AU. We set εXUV=0.30; other parameters are equal to the default values listed in Table 1. Because of the high eccentricity, the tidal evolution is calculated according to the CTL model.
FIG. 5. Three sample integrations of our code. The first-row plots show the envelope mass (left axis) and planet radius (right axis) versus time since formation; the second-row plots show the semimajor axis (left) and eccentricity (right) versus time; and the third row shows the stellar XUV luminosity (left) and stellar radius (right) versus time. The planet is initially a 1 M core with a 1 M envelope orbiting around a 0.08 M M dwarf with e=0.5 at a semimajor axis of 0.04 AU, just outside the OHZ (light blue shading). Unless otherwise noted, all other parameters are set to their default values (Table 1). As it loses mass and tidally evolves, it migrates into the CHZ (light green shading). Note that the evolution of the HZ is not shown; the CHZ and OHZ are taken to be the long-term (>1 Gyr) values. (a) In this run, we force the escape to be energy-limited, as in an X-ray-dominated flow. The planet loses its entire envelope at t≈100 Myr. (b) Same as (a), except the calculation starts at t0=100 Myr, corresponding to a planet that undergoes late migration. While the envelope still completely evaporates, this occurs at a much later time, t≈2 Gyr. (c) Same as (a), except that the escape is radiation/recombination-limited above the critical flux. Here, the envelope does not fully evaporate, and tidal migration is noticeably weaker. See the text for a discussion of the labels AF.
In column (a), we force the escape to be energy-limited (5), corresponding to an X-ray-dominated flow. Prior to the first time step, nearly 90% of the envelope mass is lost to RLO, indicated by the discontinuous drop marked A on the top plot. This is due primarily to the inflated radius shortly after formation, which reaches 30 R for a planet of age t=t0=10 Myr. Once this mass is removed, the planet enters the energy-limited escape regime, which operates on a timescale of ∼10 Myr (B). After ∼100 Myr (C), the planet loses its entire envelope and becomes a HEC.
In the center plot, we see that the planet's orbit steadily decays as it circularizes, with a sharp discontinuity in the slope at ∼100 Myr (D), corresponding to when it transitions from a gaseous (low τp) to a rocky (high τp) body. The tidal evolution from that point forward is dramatically stronger, and e decreases to ∼0.1 at 5 Gyr. The planet's semimajor axis decays by 25%, moving it well into the CHZ. As we noted earlier, the transition from low to high tidal time lags (or, alternatively, from high to low tidal quality factors) is likely to be gradual as the bulk of the energy shifts from being dissipated in the envelope to being dissipated in the core. In this case, the faster inward migration as τp increases is likely to accelerate the rate of mass loss, leading to slightly earlier evaporation times. However, given the large uncertainty in the values of τp and its dependence on planetary and orbital parameters, our current approach should suffice.
In the bottom plot, we see that the bulk of the mass loss occurs when the stellar XUV flux is high. After t≈100 Myr, the XUV luminosity is low enough that a planet with significant hydrogen left (fH≳0.01) is unlikely to completely evaporate. Here, the XUV saturation time is set to 1 Gyr, visible in the kink marked by the label E; prior to that time, the decrease in the XUV luminosity is simply a function of the rate of contraction of the star. After t∼1 Gyr (F), the stellar radius asymptotes to its main sequence value, and the XUV flux decays as a simple power law.
In column (b), we repeat the integration but delay the start time, setting t0=100 Myr. This corresponds to a planet that undergoes a late scattering event, bringing it to a=0.04 AU when both its radius and the XUV flux are significantly smaller. In this case, RLO is somewhat less effective, removing only 50% of the envelope initially (A). However, the planet still loses all its hydrogen at t≈2 Gyr (C). Interestingly, because of its delayed evaporation, the planet's eccentricity at 5 Gyr is significantly higher than in the previous run. This occurs because the transition from low to high τp (D) occurs much later. In this sense, a planet's current orbital properties can yield useful information about its atmospheric history. However, a more rigorous tidal model that accounts for the gradual change in Qp and τp as fH decreases is probably necessary to accurately infer the atmospheric evolution based on a planet's current eccentricity.
Finally, in column (c) we repeat the integration performed in (a), but this time set the flow to be radiation/recombination-limited above the critical flux (Section 2.3.5). Because of the significantly lower escape rate at early times, the envelope does not completely evaporate, and at 5 Gyr this is a super-Earth with slightly less than 1% hydrogen by mass. For a planet to completely lose its envelope in the radiation/recombination-limited regime, it must either migrate into an orbit closer to the star, have a larger eccentricity, have a smaller core, or be stripped by another process.

4.2. Dependence on Mc, MH, M, and t0

In Fig. 6, we plot the initial versus final envelope mass (MH) for planets that end up in the IHZ (red lines), CHZ (green lines), and OHZ (blue lines). Line styles correspond to different values of t0, the time at which the planet migrated into its initial close-in orbit (solid, 10 Myr; dashed, 50 Myr; dash-dotted, 100 Myr). The colored shadings are simply an aid to the eye, highlighting the spread due to different values of t0. For reference, in dotted gray we plot the line corresponding to a planet that undergoes no evaporation. Note that in most plots, the curves approach this line as the initial MH is increased: at constant core mass, it is in general more difficult to evaporate planets with more massive H/He envelopes. Finally, if a curve of a given color/line style is missing, the final hydrogen mass is zero for all values of the initial MH.
FIG. 6. Initial versus final envelope mass (MH) for planets that end up in the IHZ (red), CHZ (green), and OHZ (blue). Line styles correspond to different values of t0 (solid, 10 Myr; dashed, 50 Myr; dash-dotted, 100 Myr). Columns correspond to runs in which the escape mechanism is energy-limited (left) and radiation/recombination-limited (right); rows vary certain parameters as labeled, with all others set to their default values. In the default run, the planet has a 1 M core and orbits an M dwarf with M=0.08 M in a circular orbit. The dotted gray line corresponds to a planet that undergoes no evaporation. An X marks the critical initial envelope mass below which full evaporation occurs within 5 Gyr. In some plots, curves of a given color/line style are missing; for those runs, the entire envelope was lost for all starting values of MH.
The two columns correspond to runs in which the escape was forced to be energy-limited (left) and radiation/recombination-limited at high XUV flux and energy-limited at low XUV flux (right). Rows correspond to different planetary properties, as labeled. In the top (“Default”) row, the planet's core mass is set to 1 M. The planet orbits an M dwarf with M=0.08 M in a circular orbit. All other parameters are set to their default values.
In the second row, we double the core mass. The third row is the same as the first but for a 0.16 M star, and the fourth row is the same as the first but for an initial eccentricity e=0.4. Since planets in this run undergo orbital evolution, the different color curves correspond to planets that end up in the IHZ, CHZ, and OHZ, respectively (their initial semimajor axes are somewhat larger).
We first consider the general trends in the plots. For low-enough MH, all curves become nearly vertical. The critical envelope mass below which all mass is lost is marked with an X. Note that once planets lose sufficient mass such that MH≲10−4, the envelope is very unstable to complete erosion under the XUV fluxes considered here (corresponding to a near-vertical slope in this diagram). Many curves also display a flattening toward high initial envelope masses; some have prominent kinks beyond which the final envelope mass is constant. This is due to RLO, which causes any planets with radii larger than the Roche radius to lose mass prior to entering the HZ. Since increasing the envelope mass increases the planet radius, this results in a maximum envelope mass for some planets.
Planets that migrate in early (small t0) lose significantly more mass than planets that migrate in late. This is a consequence of both the decay in the XUV luminosity of the star with time and the quick decrease of the planet radius as the planet cools. Another interesting trend is that the difference in the evaporated amount is much more pronounced in the energy-limited runs than in the radiation/recombination-limited runs. This is due to the \(R_{ \rm p}^{3 / 2}\) scaling of the mass loss rate in the latter regime (versus the steeper \(R_{ \rm p}^3\) scaling in the former). The difference in the initial radius across runs with different values of t0 is less significant in the radiation/recombination-limited regime, resulting in more comparable mass loss rates. We also note that the \(F_{{ \rm XUV}}^{1 / 2}\) scaling of the mass loss rate in this regime results in a weaker dependence on semimajor axis, as expected: the blue, green, and red curves (for a given t0) are packed more closely together in the right column than in the left column.
Now we focus on individual rows. For the default run (a 1 M core in a circular orbit around a 0.08 M star), all planets in the IHZ lose all their hydrogen and form HECs, regardless of migration time, envelope mass, or escape mechanism. In the CHZ, only planets that migrate within ≲50 Myr and undergo energy-limited escape form HECs for all initial values of MH. However, HECs still form from planets with MH≲0.5–0.9 M in the CHZ. In the OHZ, this is only possible for planets with less than about 1% H/He by mass.
At twice the core mass (second row), all curves shift up and to the left, approaching the zero evaporation line for planets in the OHZ. In the IHZ, HECs still form from planets with any initial hydrogen amount for t0=10 Myr and in the energy-limited regime. In all other cases, HECs only form from planets with MH≲0.1 M. Of all the parameters we varied in our integrations, changing the core mass has the most dramatic effect on whether or not HECs can form. As we discuss below, HECs with masses greater than about 2 M are unlikely.
At higher stellar mass (third row), HECs are again more difficult to form, particularly in the radiation/recombination-limited regime. Due to the more distant HZ, RLO is less effective in removing mass. The shorter superluminous contraction phase of earlier M dwarfs also results in less total XUV energy deposition in the envelope. However, in the energy-limited regime, HECs still form from planets with up to 50% H/He envelopes in the IHZ.
Finally, the effect of a higher eccentricity (bottom row) is much more subtle. In general, these planets lose slightly less mass than in the default run, but the plots are qualitatively similar. At high eccentricity, the orbit-integrated mass loss rate is higher (see Section 2.3.6), but because of the orbital evolution, the planet must start out at larger a in order for it to end up in the HZ at 5 Gyr. These effects roughly cancel out: in general, HECs are just as likely on eccentric as on circular orbits. Note, also, that the green curves in the bottom right plot are the only ones that are non-monotonic. The effect is very small but hints at an interesting coupling between tides and mass loss. At high initial MH, the radius is large enough to drive fast inward migration (see below), exposing the planet to high XUV flux for slightly longer than a planet with smaller MH (and therefore a smaller radius), resulting in a change in the slope of the curve at initial MH≈0.15 M.

4.3. The role of tides

Next, we consider in detail how tides affect the evolution of HECs. We show in the Appendix that the net effect of tides is to induce inward migration and circularization of planet orbits in the HZs of M dwarfs, an effect that couples strongly to the atmospheric mass loss. For e≲0.7, the flux increases with time as planets tidally evolve, accelerating the rate of mass loss; at higher eccentricities, the flux actually decreases due to the circularization of the orbit (see the Appendix for a derivation). The changing mass and radius of the planet can then act back on the tidal evolution, either accelerating it (in cases where |dMp/dt|>>|dRp/dt|) or decelerating it in a negative feedback loop (otherwise).
In Fig. 7, we show the results of an integration of our code on a grid of tidal time lag τp versus the XUV absorption efficiency εXUV. Colors correspond to the final hydrogen mass fraction, f′H; evaporated cores occur in the white regions (f′H=0). Dark gray indicates planets that either migrated beyond the HZ or remained exterior to it and are therefore not habitable. These plots provide an intuitive sense of the relative importance of tidal evolution (y axis) and energy-limited escape (x axis) in determining whether or not a HEC is formed. We note that once a planet loses its gaseous envelope, we switch the time lag to τ′p=max(τp, 64 s), where the latter is a typical (gas-free) tidal time lag of a rocky planet (see, for instance, Barnes et al., 2013).
FIG. 7. Contours of the log of the hydrogen mass fraction f′H at 5 Gyr as a function of the tidal time lag τp and the XUV absorption efficiency εXUV for five integrations of our code. White regions correspond to planets that completely lost their envelopes; dark gray regions indicate planets that are not in the HZ at 5 Gyr. (a) The default run. The planet has a core mass of 1 M with initial fH=0.5, orbiting around a 0.08 M star in an initially highly eccentric orbit (e=0.8) at a=0.07 AU. All other parameters are set to their default values (Table 1), and the escape mechanism is radiation/recombination-limited at high XUV flux and energy-limited at low XUV flux. Note that the final envelope mass highly depends on both εXUV and τp. (b) The same as (a), but for a higher core mass Mc=1.5 M. No evaporated cores form in this scenario. (c) The same as (a), but for a lower eccentricity e=0.7. Again, no HECs form, and the planet remains outside the HZ for a larger range of τp. (d) The same as (a), but for a shorter XUV saturation time of the parent star, tsat=0.1 Gyr. No evaporated cores form. Since the XUV flux drops off much more quickly, energy-limited escape is less effective in removing mass and thus f′H is a weaker function of εXUV. (e) The same as (a), but for energy-limited escape only, which would be the case if the flow is X-ray dominated. Note that in this case whether or not the planet becomes an evaporated core is a much stronger function of both τp and εXUV.
In the default run (a), we consider a planet with a core mass of 1 M and an envelope mass of 1 M (fH=0.5) orbiting at 0.07 AU in a highly eccentric (e=0.8) orbit around a 0.08 M star. The mass loss mechanism is radiation/recombination-limited escape at high XUV flux and energy-limited at low XUV flux. At a given value of log τp, say 0, the final hydrogen fraction is a strong decreasing function of εXUV, as expected: the higher the evaporation efficiency, the smaller the final envelope. Interestingly, the dependence of f′H on the tidal time lag can be nearly just as strong. At a given value of εXUV, say 0.35, f′H depends strongly on τ, varying from 10−1.5≈3% to 0%—that is, the tidal evolution controls whether or not a HEC forms. This is due to the fact that at large τp tidal migration is fast, bringing the planet into the HZ while its radius is still inflated and the stellar XUV luminosity is higher. Planets with lower values of τp migrate in later and undergo slower mass loss. At very low τp, tidal migration is too weak to bring the planets fully into the HZ. The effect is stronger for higher εXUV: these are planets whose radii decrease very quickly (due to the fast mass loss), slowing down the rate of tidal evolution and keeping them outside the HZ for longer.
In plot (b), we increase the core mass slightly to 1.5 M. The effect on the mass loss is significant, and HECs no longer form. At high εXUV and high τp, the lowest value of f′H is about 1%. This reinforces what we argued above: HECs are most likely for the lowest-mass cores.
In plot (c), we instead decrease the eccentricity to 0.7. As in (b), HECs no longer form in these runs due to the decreased strength of the tidal evolution, and again the minimum f′H is about 1%. Note that this does not mean that HECs are more likely at higher eccentricity in general—this is only the case here because we fix the initial semimajor axis at 0.07 AU. At fixed final semimajor axis (i.e., what we can readily observe in actual systems), planets on initially circular orbits will still lose more mass than planets on initially eccentric orbits (which must have formed farther out).
In plot (d), we decrease the saturation time to tsat=0.1 Gyr, which is typical of early M/late K dwarfs (Jackson et al., 2012). No HECs form, and the final hydrogen fraction is less sensitive to both τp and εXUV: in general, mass loss is significantly suppressed.
Finally, in plot (e), we force the escape mechanism to be energy-limited only. Since the escape is now entirely controlled by εXUV, the dependence on this parameter is naturally much stronger, and complete evaporation occurs in this case for εXUV≳0.32 at any τp. Because evaporation occurs more quickly in this case than in the other plots, at any given time the planet radii are smaller, resulting in less efficient migration at a given τp. More planets therefore do not make it into the HZ. Interestingly, for very large εXUV (≳0.4), all planets make it into the HZ. This is due to the fact that these planets transition from gaseous (low τp) to gas-free (high τp, equal to 64 s in these runs) very early on. Despite their lower radii, they benefit from the stronger tidal dissipation of fully rocky bodies and are able to make it into the HZ after 5 Gyr.
In general, the coupling between tides and mass loss is quite complex. For planets with high initial eccentricities, this coupling can ultimately determine whether or not a HEC will form. We discuss this in more detail in Section 4.4.

4.4. Evaporated cores in the habitable zone

Having shown that HECs are possible in certain cases, we now wish to explore where in the HZ we may expect to find them. For a “terrestrial” (we use this term rather loosely2) planet detected in the HZ, it would be very informative to understand whether or not it may be the evaporated core of a gaseous planet, since its past atmospheric evolution may strongly affect its present ability to host life. To this end, we ran eight grids of 2.7×106 integrations each of our evolution code under different choices of parameters in Table 1. For initial semimajor axes in the range 0.01≤a0≤0.5, initial eccentricities in the range 0≤e0≤0.95, and stellar masses between 0.08 MM≤0.5 M, we calculate the final (i.e., at tstop=5 Gyr in the default run) values of a, e, and the envelope mass MH. We then plot contours of MH as a function of the stellar mass and the final semimajor axis and eccentricity (that is, the observable parameters). Since, in principle, the mapping (a0, e0, MH0)→(a, e, MH) is not necessarily a bijection (due to the nonlinear coupling between tides and mass loss and the finite resolution of our grid), the function MH(a, e) may be multivalued at some points. In these cases, we take MH at coordinates (a, e) to be the minimum of the set of all final values of MH that are possible at (a, e). Because of this choice, the MH=0 contour in a versus M plots (Figs. 8–10) separates currently terrestrial planets that could be the evaporated cores of mini-Neptunes (left) from currently terrestrial planets that have always been terrestrial (right). In other words, we are showing where evaporated cores are ruled out.
FIG. 8. Regions of parameter space that may be populated by HECs, for Mc=1 M, initial MH≤1 M, and default values for all other parameters. Terrestrial planets detected today occupying the space to the left of each contour line could be the evaporated cores of gaseous planets with fH≤0.5. Planets detected to the right of the contour lines have always been terrestrial/gaseous. Dark red lines correspond to the conservative mass loss scenario, in which mass loss is radiation/recombination-limited at high XUV flux and energy-limited at low XUV flux. Dark blue lines correspond to mass loss via the energy-limited mechanism only. Planets around stars with significant X-ray emission early on are likely to be in the latter regime. Different line styles correspond to different eccentricities today. Terrestrial planets detected at higher eccentricity (dashed and dash-dotted lines) could be evaporated cores at slightly larger orbital separations than planets detected on circular orbits (solid lines). Note that in the energy-limited regime, all 1 M terrestrial planets in the HZ of low-mass M dwarfs could be HECs. At higher stellar mass, HECs are restricted to planets in the CHZ and IHZ. In the radiation/recombination-limited regime, the accessible region of parameter space is smaller, but around the lowest-mass M dwarfs HECs are still possible in the CHZ.
FIG. 9. The same plot as Fig. 8 but for a conservative choice of the parameters governing atmospheric escape: a short XUV saturation time tsat=0.1 Gyr and a low XUV absorption efficiency εXUV=0.15. In this case, HECs are no longer possible for radiation/recombination-limited escape. For energy-limited escape, HECs are only possible in the IHZ of low-mass M dwarfs and in the CHZ of M dwarfs at the hydrogen-burning limit. At high eccentricity (e=0.5), HECs are only marginally more likely.
FIG. 10. The same as Fig. 8 but for a higher core mass Mc=2 M. HECs are now confined mostly to the IHZ of low-mass M dwarfs in the energy-limited regime. For planets undergoing radiation/recombination-limited escape, super-Earth HECs may not be possible.
In Fig. 8, we perform the calculations described above for planets with Mc=1 M and initial fH=0.5 (i.e., initial MH=1 M). The HZ is plotted in the background for reference, where again the IHZ, CHZ, and OHZ are indicated by the red, green, and blue shading, respectively. Line styles correspond to different final values of the eccentricity (0, solid; 0.25, dashed; and 0.5, dash-dotted). Dark red lines correspond to the MH=0 contours in the radiation/recombination-limited escape model; dark blue lines correspond to the energy-limited model. Note that, though we designate our two models as “energy-limited” and “radiation/recombination-limited,” at low XUV fluxes the escape is energy-limited in both models, since below Fcrit the energy-limited escape rate is always smaller than the radiation/recombination-limited escape rate. In this sense, the “radiation/recombination-limited model” is always conservative, as the mass loss rate is set to the minimum of Eqs. 5 and 13.
As an example of how to interpret this figure, consider a 1 M rocky planet discovered orbiting a 0.2 M star at 0.1 AU, that is, squarely within the CHZ. Since this planet lies to the left of all blue curves, under the assumption of energy-limited escape it could be a HEC. On the other hand, if the atmospheric escape were radiation/recombination-limited, this planet must have always been terrestrial.
Next, consider a putative rocky planet discovered around the same 0.2 M star skirting the IHZ (i.e., at 0.07 AU). Under the energy-limited assumption, this planet could be a HEC. For radiation/recombination-limited escape, however, whether or not it could be a HEC depends on its present eccentricity. If the planet is currently on a circular orbit, we infer that it has always been terrestrial (since it lies to the right of the e=0 contour). However, since the planet lies to the left of the higher eccentricity contours, if e≳0.25, it could be a HEC.
We conclude from Fig. 8 that, if energy-limited escape is the dominant mechanism around M dwarfs, planets with Mp∼1 M in the CHZ and IHZ of these stars could be HECs. For M≲0.15 M, HECs may exist throughout the entire HZ. If, on the other hand, these planets are shaped mostly by radiation/recombination-limited escape, HECs are only possible in the IHZ for M≲0.2 M and in the CHZ of M dwarfs near the hydrogen-burning limit.
The effect of the eccentricity is significant primarily in the radiation/recombination-limited regime, where a lower mass loss rate keeps the planet's radius large for longer than in the energy-limited regime, allowing it to migrate into the HZ faster, thereby enhancing the mass loss rate. In some cases, particularly near the IHZ, the present-day eccentricity yields important information about a planet's evolutionary history: depending on the precise value of e, a given terrestrial planet may or may not be a HEC. However, we urge caution in interpreting the results in Fig. 8 at nonzero eccentricity, given the large uncertainty in the tidal parameters of exoplanets.
On this note, it is important to bear in mind that the curves in Fig. 8 are a function of our choice of parameters in Table 1. To assess the impact of our choice of “default” parameters on these contours, in Fig. 9 we repeat the calculation for more conservative values of the two parameters that govern the mass loss mechanism: the XUV saturation time tsat and the absorption efficiency εXUV. In this figure, we choose tsat=0.1 Gyr, the nominal value for earlier K/G dwarf stars (Section 2.1), and εXUV=0.15.
In this grid, all curves shift quite dramatically to lower a, and HECs are no longer possible under radiation/recombination-limited escape. For energy-limited escape, HECs are confined to the IHZ for M≲0.15 M and the CHZ for M≲0.1 M at all eccentricities considered here.
Given the large difference between the results of Figs. 8 and 9, care must be taken in assessing whether a planet may be a HEC. Since it is likely that the XUV saturation time is much longer for M dwarfs than for K and G dwarfs, and since our goal at this point is to separate regions of parameter space where HECs can and cannot exist, Fig. 8 is probably the more relevant of the two. We discuss this point further in Section 5.
In Fig. 10, we raise the core mass to Mc=2 M. HECs are now possible only in the IHZ and only in the energy-limited regime. At the lowest M, HECs may be possible in the CHZ (energy-limited) and at the very inner edge of the HZ (radiation/recombination-limited). At higher eccentricity, the parameter space accessible to HECs is slightly higher in the energy-limited regime, but it is still a small effect overall. For a conservative choice of escape parameters (tsat=0.1 Myr and εXUV=0.15) with Mc=2 M, no HECs form anywhere in the HZ.

5. Discussion

5.1. Initial conditions

We have shown that HECs can form from small mini-Neptunes (Mp≲2 M) with H/He mass fractions as high as fH=0.5 around mid to late M dwarfs. However, instrumental to our conclusions is the assumption that these planets formed beyond the snow line and migrated quickly into the HZ at t=t0. Planets that form in situ in the HZ are unlikely to accumulate substantial H/He envelopes due to large disk temperatures and relatively long formation timescales. Significant gas accretion can only occur prior to the dissipation of the gas disk, which occurs on a timescale of a few to ∼10 Myr; in particular, Lammer et al. (2014) showed that a 1 M terrestrial planet generally accretes less than ∼3% of its mass in H/He in the HZ of a solar-type star. Moreover, planets that form in situ do not undergo RLO, since accretion of any gas that would lead to Rp>RRoche simply will not take place.
But can mini-Neptunes easily migrate into the HZ? The large number of recently discovered hot Neptunes and hot super-Earths (e.g., Howard et al., 2012) suggests that inward disk-driven migration is a ubiquitous process in planetary systems, as it is highly unlikely that these systems formed in situ (Raymond and Cossou, 2014). Moreover, systems such as GJ 180, GJ 422, and GJ 667C each have at least one super-Earth/mini-Neptune in the HZ (Anglada-Escudé et al., 2013; Tuomi et al., 2014), some or all of which may have migrated to their current orbits. On the theoretical front, N-body simulations by Ogihara and Ida (2009) show that migration of protoplanets into the HZ of M dwarfs is efficient due to the proximity of the ice line and the fact that the inner edge of the disk lies close to the HZ. Mini-Neptunes that assemble early beyond the snow line could in principle also migrate into the HZ, but the migration probability and the distribution of these planets throughout the HZ needs to be investigated further in order to constrain the likelihood of formation of HECs.
A second issue is whether or not Earth-mass cores can accrete large H/He envelopes in the first place. One of our key results is that Earth-mass HECs can form from mini-Neptunes with initial envelope mass fractions of up to 50%; this would require a 1 M planetary embryo to accrete an equivalent amount of gas from the disk. While such a planet would be stable against RLO beyond the snow line, whether or not it could have formed is unclear. Gas accretion takes place in two different regimes: (i) a slow accretion regime, in which the envelope remains in hydrostatic equilibrium and gains mass only as it cools and contracts, evacuating a region that is then filled by nebular gas; and (ii) a runaway accretion regime, in which the rapidly increasing mass of the envelope leads to an increase in the size of the Hill sphere and increasingly faster gas accretion (e.g., Pollack et al., 1996). Since the critical core mass for runaway accretion is thought to be somewhere in the vicinity of 10–20 M (Pollack et al., 1996; Rafikov, 2006), the progenitors of HECs must accrete their gas slowly. As we mentioned above, Lammer et al. (2014) found that Earth-mass planets typically form with fH≲0.03 in the HZ of solar-type stars. This is consistent with recent analyses of Kepler planet data by Rogers (2014) and Wolfgang and Lopez (2014), who found that most planets with radii less than about 1.5 R (corresponding to masses less than ∼5 M) are rocky, with typical H/He mass fractions of about 1%.
Formation beyond the snow line can increase fH: Rogers et al. (2011) showed that a planet with core mass 2.65 M accretes only 0.54 M of gas (corresponding to fH≈17%). Bodenheimer and Lissauer (2014), on the other hand, showed that planets with cores in the range 2.2–2.5 M generally accrete less than 10% of their mass in H/He. However, these authors terminated gas accretion at 2 Myr. For a longer disk lifetime of 4 Myr, Bodenheimer and Lissauer (2014) showed that a planet with core and envelope masses of 2.8 M and ∼1.2 M, respectively, can form (fH≈30%). This is particularly relevant to planet formation around M dwarfs, as these stars may have disk lifetimes in excess of 5 Myr (Carpenter et al., 2006; Pascucci et al., 2009), which could allow for larger initial H/He envelope fractions.
Nevertheless, while a 1 M core with a 1 M envelope is probably an unlikely initial condition, it is important to keep in mind that our results apply to planets that form with smaller H/He fractions as well. Figures 8–10 show where planets with up to 50% H/He can form HECs; planets with lower H/He mass fractions also form HECs at the same values of a and M. See Section 5.6 for a more detailed discussion.

5.2. Are HECs habitable?

Under the core accretion mechanism, mini-Neptunes are likely to form close to, or beyond, the snow line, where disk densities are higher due to the condensation of various types of ices; these planets should therefore have large quantities of volatiles, including water, ammonia, and CO2 ices. If we assume a disk composition similar to that around the young Sun, the bulk composition of their cores would probably be similar to that of comets: roughly equal parts ice and silicate rock, similar to what studies predict for the composition of water worlds (Léger et al., 2004; Selsis et al., 2007). Once these planets have migrated into the HZ and lost their envelopes, it is quite likely that they would be water worlds and therefore not “terrestrial” in the strictest sense of the word. Whether or not such planets are habitable may depend on their ability to sustain active geochemical cycles, which are crucial for life on Earth today.
One concern is a possible interruption of the carbon cycle by a high-pressure ice layer at the bottom of the ocean (Léger et al., 2004). Recently, Alibert (2014) calculated the critical ocean mass for high-pressure ice formation, finding that it lies between 0.02 and 0.03 M for an Earth-mass planet. If HECs are, in fact, cometlike in composition, a deep ice layer would separate their oceans from their mantles, which could inhibit the recycling of carbon and other bioessential elements between the two reservoirs, a process that is critical to life on Earth. However, whether this is the case is far from settled, as processes involving solid state ice convection could mediate the cycling of these elements. In particular, Levi et al. (2013) and Kaltenegger et al. (2013) presented a mechanism that could recycle CH4 and CO2 between the interior and the atmosphere of water worlds, invoking the ability of these molecules to form clathrates that could be convectively transported through the ice to the surface. Other mechanisms could also exist, and without further modeling, it is unclear whether a high-pressure ice layer poses a threat to habitability.
The probable difference in composition between HECs and Earth is likely to have other geophysical implications. Hydrogen-rich compounds such as methane and ammonia could make up a non-negligible fraction of a HEC's mass, leading to extremely reducing conditions at the surface. These could also end up in a secondary atmosphere along with large quantities of CO2, which could lead to strong greenhouse heating, although it is likely that most of the CO2 and NH3 would be sequestered in the ice mantle (Léger et al., 2004). The compositional difference of HECs would also likely lead to mantle convection and tectonic activity different from Earth's, as well as differences in magnetic field generation by a possible core dynamo. Since both an active tectonic cycle and a magnetic field may be necessary for habitability, these issues need to be investigated further.
Also critical to the habitability of a HEC is its ability to outgas a secondary atmosphere once its primordial envelope is lost. In particular, the secondary atmosphere must be stable against erosion. While the power-law decay in XUV emission after tsat could allow for such an atmosphere to form, low-mass M dwarfs may remain active for t >>1 Gyr. Strong planetary magnetic fields could potentially shield these atmospheres; Segura et al. (2010) showed that flares on the extremely active M dwarf AD Leo would have a small effect on the atmosphere of an Earth-like planet in the HZ provided it has an Earth-like magnetic field. However, interactions with the stellar wind could still pose serious problems to these and other planets in the HZs of M dwarfs. In particular, the high XUV/EUV fluxes of active M dwarfs can lead to significant expansion of the upper atmosphere, potentially causing the radius of the exobase to exceed the distance to the stellar magnetopause (Lichtenegger et al., 2010; Lammer et al., 2011). For a pure N2 atmosphere, Lichtenegger et al. (2010) showed that nitrogen ionized by EUV radiation above the exobase is subject to ion pickup by the solar wind, leading to the complete erosion of a 1 bar atmosphere in as short a time as 10 Myr. A stronger planetary magnetic field or large quantities of a heavier background gas such as CO2 may be necessary to suppress ion pickup and preserve secondary atmospheres on HECs.
Given the complex processes governing the habitability of HECs, a detailed investigation of these issues is left to future work. As we discuss below, stronger constraints on the details of the X-ray/EUV evolution of M dwarfs of all masses are essential to understanding the fate and ultimately the habitability of planets around these stars.

5.3. The need for better constraints

Figure 8 shows that the formation of HECs depends strongly on the atmospheric escape mechanism. As we mentioned earlier, whether a flow is closer to radiation/recombination-limited or energy-limited will depend on the ratio of the X-ray to EUV luminosity of the parent star. Low-mass M dwarfs may have XUV luminosities as high as ∼4×1029 erg/s early on (Section 2.1). If X-rays contribute more than a few percent of this luminosity, low-mass low-density planets in the HZs of these stars may undergo energy-limited X-ray-driven escape (Fig. 11 in Owen and Jackson, 2012). Unfortunately, knowledge of the exact age-luminosity relation in the X-ray and EUV bands for M dwarfs is still very poor, in great part because of the large uncertainties on these stars' ages. However, recent studies suggest that X-rays contribute a significant fraction of this luminosity (for a review, see Scalo et al., 2007). In particular, Stelzer et al. (2013) reported high (LX≳1029 erg/s) X-ray luminosities for a sample of early active M dwarfs in the solar neighborhood, with a steeper age dependence than in far ultraviolet and near ultraviolet bands, which dominate the emission for t≳1 Gyr. This is consistent with Owen and Jackson (2012), who argued that close-in mini-Neptunes may undergo a transition from X-ray-driven escape at early times to EUV-driven escape at later times.
FIG. 11. Similar to Fig. 8, but here contours correspond to different choices of final fH (which we denote fH) for a 1 M core on a circular orbit. The solid lines correspond to f′H=0 and are the same as the e=0 contours in Fig. 8. Dashed lines correspond to the transition between planets that have less than (left) and more than (right) 0.1% H/He (f′H=0.001) at 5 Gyr. Dotted lines correspond to the 1% H/He (f′H=0.01) transition. While the f′H=0.001 contours are barely distinguishable from the f′H=0 contours, at final f′H=0.01 the HEC parameter space is significantly larger. However, it is unclear whether planets with f′H=0.01 could be habitable.
Moreover, atmospheric X-ray heating and cooling is primarily done by metals. As Owen and Jackson (2012) pointed out, atmospheric composition is also likely to play a role in determining whether hydrodynamic flows are EUV- or X-ray-driven. Additionally, the presence of dust in the envelopes of these planets could greatly affect their ability to cool via Lyman α radiation. Absorption of recombination radiation by dust particles lifted high into the envelope by vigorous convection could convert a significant fraction of this energy into heating, which could effectively increase the absorption efficiency εXUV and bring the flow closer to the energy-limited regime. Unfortunately, our present parametric escape model is unable to address the effect of dust and metal abundances on the escape rate—this issue needs to be revisited in the future with more sophisticated hydrodynamic models.
The difference between Fig. 8 and Fig. 9 is just as significant. At lower tsat and εXUV, the HEC parameter space is greatly reduced. As we discuss in Section 2.1, the lower value tsat=0.1 Gyr is more representative of K and G dwarfs, while M dwarfs may remain saturated for tsat>1 Gyr. In fact, the energy-limited contours in Fig. 8 closely trace the CHZ/OHZ boundary as M increases, predicting that HECs may even be possible around solar-type stars (not shown). This is of course not the case, since solar-type stars are known to leave the saturation phase around tsat=0.1 Gyr (Ribas et al., 2005). Above a certain value of the stellar mass near the M/K dwarf transition (0.5 MM≲0.7 M), Fig. 9 becomes a better description for HECs; but for low-mass M dwarfs, Fig. 8 is probably more appropriate.
Naturally, the exact value of εXUV will also affect the possible distribution of HECs within the HZ. Recently, Shematovich et al. (2014) performed numerical simulations to solve the kinetic Boltzmann equation for XUV-irradiated hydrogen atmospheres, finding an upper limit to the heating efficiency of εXUV≈0.20. Our nominal value εXUV=0.3 may therefore be an overestimate for many mini-Neptunes. However, given that our goal in this study is to explore where in the HZ HECs are possible and to map regions where the transition from gaseous to terrestrial is not possible, our present approach should suffice. Nevertheless, further studies constraining εXUV are crucial to improving our understanding of this parameter space.

5.4. Eccentricity effects

From Figs. 8–10 we see that as we consider planets with higher current eccentricity, the region where evaporated cores are possible overlaps increasingly more of the HZ, particularly in the radiation/recombination-limited regime. This does not necessarily mean that HECs are more likely at high eccentricity—but rather that planets in circular orbits at certain a cannot be HECs, while planets at higher eccentricity can. However, that being said, there is an interesting negative feedback between tides and mass loss that may enhance the probability of HECs on eccentric orbits. A planet with initially high eccentricity will undergo fast tidal evolution (due to the strong β dependence in Eq. 74), particularly if its radius is large (since de/dt\(R_{ \rm p}^5\)). If a planet loses mass quickly (as is usually the case for gaseous planets with large Rp and at large eccentricity), its radius will shrink, and de/dt will decrease. The earlier a planet sheds its envelope, the more likely it is to maintain a nonzero eccentricity in the long run (i.e., after ∼5 Gyr). Since HECs typically form from such quickly evaporating planets, they are more likely to end up in eccentric orbits than their gaseous counterparts that did not lose their envelopes. There is, of course, a trade-off here in the sense that gas-free planets should have higher τp (and therefore faster de/dt) than gaseous planets, but the dependence on τp is linear and therefore much weaker than the dependence on the radius. What this means is that there is an interesting link between present eccentricity and mass loss history. Translating a planet's present eccentricity into a probability that it is an evaporated core is no easy task, however, and likely requires a detailed understanding of its initial orbital state and the migration mechanism (or mechanisms) it underwent. On the other hand, statistical surveys of planets found to the left of the contours in Fig. 8 could uncover interesting trends.
We note that Figs. 8–10 show eccentricity contours only up to e=0.5. At final eccentricities higher than 0.5 today, the contour lines should move farther to the right, increasing the parameter space populated by HECs. However, since high eccentricities today in general require extremely high eccentricities in the past, we do not calculate contours above this level. We note, finally, that even though we run simulations with e0 as high as 0.95, Figs. 8–10 remain unchanged for a lower cutoff e0≲0.7.

5.5. Orbital effects due to Roche lobe overflow

In the present work, we neglect any orbital effects due to RLO, which could in some cases lead to significant outward migration of the planet; see (22). We argued in Section 3.3 that since RLO should occur during the initial stage of migration (i.e., from beyond the snow line into the HZ), modeling it was outside the scope of our study. However, a careful investigation of this initial migration process could uncover interesting couplings. For instance, mini-Neptunes that initially overshoot the HZ would experience strong RLO, which could in principle cause them to migrate back into the HZ. These planets could have lost significantly more mass than the ones we considered here, since both RLO and XUV-driven escape would be stronger in their closer-in orbits.
A second important point concerning RLO-induced migration is that for nonzero eccentricity (22) does not apply. In this case, it can be shown that angular momentum exchange between the gas and the planet at pericenter (where overflow is strongest) leads to a net increase in the eccentricity. In the limiting case that both bodies may be treated as point masses and the mass transfer rate may be approximated as a delta function at pericenter (which is appropriate for high e), Sepinsky et al. (2009) showed that
\begin{align*} \frac { da } { dt } = - \frac { a } { \pi } \frac { \dot { M } _p } { M_p } ( 1 - e^2 ) ^ { 1 / 2 } \left( 1 - \frac { M_p } { M_\star } \right) \tag { 23 } \end{align*}
\begin{align*} \frac { de } { dt } = - \frac { 1 } { \pi } \frac { \dot { M } _p } { M_p } ( 1 - e^2 ) ^ { 1 / 2 } ( 1 - e ) \left( 1 - \frac { M_p } { M_\star } \right) \tag { 24 } \end{align*}
predicting that both a and e will tend to increase with time (recall that \(\dot{M}_{ \rm p}\) is negative). Their relative rates of change are
\begin{align*} \frac { de } { dt } = ( 1 - e ) \frac { 1 } { a } \frac { da } { dt } \tag { 25 } \end{align*}
In other words, at intermediate values of e, the fractional rates of change in the semimajor axis and the eccentricity are comparable, and the eccentricity will increase proportionally to the semimajor axis.
As e increases, so does the atmospheric mass loss rate (via Kecc), the RLO rate (modulo the increase in a), and the rate of tidal evolution, leading to potentially faster mass loss and rich couplings between these different processes. We will consider all these effects in a future paper.

5.6. Different fH

Figures 8–10 correspond to planets with initial hydrogen fraction fH=0.5. Planets with lower initial fH are naturally more unstable to complete loss of their envelopes—this would move all contours to the right, expanding the region where HECs are possible. However, for a different choice of initial fH in the range 0.1≲fH≲0.5, the figures change very little: at constant core mass, it is only marginally harder to fully evaporate a fH=0.5 envelope than it is to evaporate a fH=0.1 envelope (see discussion below). This is due to both RLO, which dramatically reduces the envelope mass early on, and the fast atmospheric escape for extremely inflated planets. In Fig. 5a, for instance, the envelope fraction decreases by a factor of 10 within the first ∼10 Myr due to energy-limited escape. The escape process is in general very fast for large fH and bottlenecks for fH≲0.1; thus the HEC boundary is relatively insensitive to the exact choice of the initial fH.
Along the same lines, we can ask how the choice of final fH (which here we denote fH) affects our results. In this study, we define an evaporated core as a planet with f′H=0 and no atmosphere. However, mini-Neptunes with substantial gaseous envelopes that evaporate down to f′H∼0.01 become fundamentally solid planets. At such low f′H, additional nonthermal mass loss processes, including flaring events and interactions with the stellar wind, may significantly erode the remaining envelope. Investigating the habitability of planets with nonzero f′H is beyond the scope of this paper, but it is worthwhile to consider how an alternative definition of a “habitable evaporated core” affects our results. To this end, in Fig. 11 we plot contour lines corresponding to three different choices of the critical f′H below which we consider a planet to be a HEC. Solid lines correspond to f′H=0 (the default choice); dashed lines correspond to f′H=0.001; dotted lines correspond to f′H=0.01. Note that in the radiation/recombination-limited regime, this choice does not significantly affect evaporated cores in the HZ. However, for f′H=0.01, all planets undergoing energy-limited escape in the HZ could be HECs.
Earths and super-Earths with up to a few percent of their mass in H/He may still harbor liquid water oceans, so at this point their habitability cannot be ruled out. At higher f′H, however, the surface pressures in excess of ∼1 GPa (e.g., Choukroun and Grasset, 2007) result in the formation of high-pressure ices, at which point the planet will likely no longer be habitable.

5.7. Other atmospheric escape mechanisms

We modeled the XUV emission of M dwarfs as a smooth power law, but active M dwarfs are seldom so well behaved. Frequent flares and coronal mass ejections punctuate the background XUV emission and will lead to erosion of planetary atmospheres beyond what we model here. During flares, the ratio of the X-ray to bolometric luminosity (LX/Lbol) can increase by up to 2 orders of magnitude (Scalo et al., 2007). Moreover, interactions with the stellar wind and nonthermal escape mechanisms such as ion pickup and charge exchange can lead to a few Earth-ocean equivalents of hydrogen from these planets (Kislyakova et al., 2013, 2014). However, as Kislyakova et al. (2013) demonstrated, the escape rate due to these processes generally amounts to only a few percent of the hydrodynamic escape rate. Including these nonthermal mechanisms would thus have a minor effect on our results.
In this study, we also neglect the effects of magnetic fields. The presence of a strong planetary magnetic field may inhibit escape during flares and possibly decelerate the mass loss rate if the planetary wind is ionized. Since a strong magnetic field is likely a requirement of surface habitability around an M dwarf (Scalo et al., 2007), this point must be addressed in future work.
Finally, stars less massive than about 0.1 M may be in a “supersaturation” regime early on, saturating at XUV fluxes 1 or even 2 orders of magnitude below higher-mass M dwarfs (Cook et al., 2014). This could reduce escape rates from planets around M dwarfs close to the hydrogen burning limit.

5.8. Multiplanet systems

In this paper, we restricted our calculations to single-planet systems. However, many of the super-Earths and mini-Neptunes recently detected by Kepler reside in multiplanet systems, such as Kepler-32 (Fabrycky et al., 2012), Kepler-62 (Borucki et al., 2013), and Kepler-186 (Quintana et al., 2014). Moreover, Swift et al. (2013) demonstrated how the five-planet system Kepler-32 could be representative of the entire ensemble of Kepler M dwarf systems, arguing that multiplicity could be the rule rather than the exception for these stars.
While tidal processes still operate in systems with multiple planets, the orbital evolution of individual planets will be much more complex. In general, though tidal dissipation may still lead to a decrease in the semimajor axis, the eccentricity evolution will be governed by secular interactions between the planets. The coupling between mass loss and tidal evolution we investigated in Section 4.3 would likely be weaker, particularly for closely packed coplanar systems where the eccentricities are necessarily low. However, planet-planet interactions could allow for even richer couplings to the atmospheric evolution of planets in the HZ.
Consider, for instance, the case of a terrestrial planet in the HZ and a mini-Neptune interior to the HZ. Both mass loss and tidal dissipation scale inversely with the semimajor axis, so both effects could strongly shape the fate of the inner planet, which could in turn affect both the orbital and atmospheric evolution of the HZ planet in complex ways. The same could happen for a planet overflowing its Roche lobe interior to the HZ, which as we argued in Section 5.5 could experience large changes in a and e.
Similarly, consider a scenario in which a planet in the HZ loses mass at the same time as it is perturbed by a more massive companion exterior to the HZ. As the planet's mass decreases, it will undergo larger swings in eccentricity, which in turn lead to faster mass loss and stronger orbital evolution.
Another interesting scenario involves planets close to mean motion resonances. For a system of two planets in resonance, Lithwick and Wu (2012) and Batygin and Morbidelli (2013) showed that the presence of a dissipative mechanism such as tidal evolution naturally leads to the repulsion of the two planets' orbits: the inner planet's orbit decays, while the outer planet gets pushed outward. In the hypothetical case of two mini-Neptunes in resonant orbits interior to the HZ, this mechanism could result in the migration of the outer planet into the HZ, having lost significantly more mass (due to its initially smaller orbit) than if it had originated exterior to the HZ and migrated inward. A chain of resonant planets, which is a potential outcome of disk-driven migration (Terquem and Papaloizou, 2007; Ogihara and Ida, 2009), could have similarly interesting consequences on planets in the HZ.
All these scenarios, which involve coupling between atmospheric mass loss, tidal evolution, and secular planet-planet interactions, probably occur even for planets outside the HZ. Modeling this coupling could be critical to understanding systems like Kepler-36 (Carter et al., 2012), where two planets with semimajor axes differing by ∼10% have a density ratio close to 8, which Lopez and Fortney (2013) showed could be the result of starkly different mass loss histories. Future work should investigate how mass loss modifies the orbital interactions in multiplanet systems.

5.9. Other caveats

5.9.1. The habitable zone

Yang et al. (2013) argued that cloud feedback on tidally locked planets can greatly increase the planetary albedo and move the IHZ in by a substantial amount. Similarly, Abe et al. (2011) showed that planets with limited surface water are stable against runaway greenhouses at greater insolation than Earth, which also decreases the IHZ distance. Since evaporated cores are more likely at smaller a, a closer-in inner edge to the HZ could greatly increase the parameter space available to these planets. However, as we argued in Section 5.2, HECs are likely to have abundant surface water.
Luger and Barnes (2014) showed that terrestrial planets in the HZ of M dwarfs may experience long runaway greenhouses during their host stars' extended PMS phases, during which time water loss to space can lead to their complete desiccation, rendering them uninhabitable. HECs are naturally more robust against severe water loss, given that their initially dense H/He envelopes can shield the surface and lower atmosphere from XUV radiation. Once the envelope is lost, water loss from the surface could occur, but by that point the stellar XUV flux will be much lower. Future work will address the fate of the water on HECs in detail.

5.9.2. Thermal evolution

As mentioned before, we do not consider how the mass loss rate affects the planet radius, but instead assume the radius instantaneously returns to the “equilibrium” value for the given age, mass, and hydrogen fraction. In reality, the radius is likely to remain inflated for some time, particularly for fast mass loss (Lopez et al., 2012). Because our radii decrease too quickly, the mass loss feedback on the tidal evolution is always negative; mass loss always acts to dampen the tidal evolution. However, accurate modeling of the radius could enable a positive feedback, in which the effect of the decreasing mass (higher da/dt) overpowers the effect of the decreasing radius (lower da/dt) in (60) and (73), leading to faster orbital decay.
Currently, our planet radii are also independent of both the insolation and the degree of tidal heating. Under the extremely high fluxes and strong tidal forces at early times, planets in the HZs of M dwarfs are likely to be significantly more inflated than modeled here, resulting in faster atmospheric mass loss and likely a higher probability of complete envelope loss. We also do not model radiogenic heating, which could further add to the inflation of the planet. In this sense, our results are conservative, and HECs may in fact be possible at larger distances from their parent stars. A self-consistent thermal evolution model must be developed to accurately address this point.
We also emphasize that the quasi-static approximation (see the Appendix) may not be valid for the most inflated planets on eccentric orbits above e≈0.4. Given the large uncertainty concerning the orbital migration due to RLO for planets on eccentric orbits, we urge caution in interpreting our results quantitatively for highly eccentric planets.

5.9.3. Tidal evolution

Given that the CPL tidal model is accurate only to second order in eccentricity, most of our integrations were performed under the CTL framework. At low e, however, these models predict qualitatively similar results for both the orbital migration and the coupling to the mass loss processes. At high e, on the other hand, tidal models lack observational validation, given the relatively low eccentricities of the major bodies in our solar system. We again urge caution in interpreting quantitative results for e≳0.5.
Recall that our values of Qp and τp are fixed throughout the integrations and change discontinuously to the more dissipative rocky values as soon as the envelope is lost; this is certainly an oversimplification. Though we defined Q via its relationship to the phase lag in (19), it may also be defined as the specific dissipation function (Goldreich and Soter, 1966), defined as
\begin{align*}Q^ { - 1 } \equiv \frac { 1 } { 2 \pi E_0 } \int_0^ { 2 \pi } \left( - \frac { dE } { dt } \right) dt \tag { 26 } \end{align*}
where E0 is the maximum tidal energy stored in the system, and the integral is the energy dissipated per cycle. Q is therefore inversely proportional to the dissipation per unit energy stored. Planets with large values of fH, where the dissipation throughout most of the body is small, should therefore have large values of Qp; rocky planets with small fH, conversely, have large specific dissipation rates and correspondingly low values of Qp, in agreement with the remarks in Section 2.4.1. Therefore, we would expect that as planets lose mass, Qp should decrease, corresponding to an increase in τp. Modeling Qp(fH) and τp(fH), however, is outside the scope of this paper and will be considered in future work. We simply note that, as planets lose mass, the decreasing value of Qp/increasing value of τp should offset the decreasing radii, leading to stronger tidal evolution and stronger couplings than reported here.
One must likewise be careful in choosing the planet radius that goes into calculating the tidal evolution. This radius should be the effective radius of the dissipating material and is thus dependent on the chosen value of Qp and τp; for high Qp/low τp, the radius should probably be that of the envelope. We took this to be Rp, the 20 mbar radius shown in the tracks in Fig. 3; other choices of the tidal radius will lead to different evolutions.
We only consider planets with zero initial obliquity. For nonzero obliquity, the tidal evolution equations must be modified and will lead to differences in the evolution; see the works of Heller et al. (2011) and Barnes et al. (2013). We also ignore the spin-up of the planet as it thermally contracts over time, since it should return to the equilibrium spin almost instantaneously. The excess angular momentum would most likely be absorbed into the orbit, but it is a small enough fraction of the orbital angular momentum that it can be safely neglected.
Finally, we note that throughout this paper we set tstop=5 Gyr, based on the age of the Solar System. For systems older than 5 Gyr, our results should change little, since both the tidal and atmospheric evolution will have strongly tapered off by this time. Younger systems, however, may still be in the throes of tidal decay and mass loss processes, and this must be kept in mind when searching for HECs.

6. Conclusions

We have shown that gas-rich mini-Neptunes that migrate early (t ≲ 10 Myr) into the HZs of M dwarfs can naturally shed their gaseous envelopes and form gas-free Earth-mass planets. Together, RLO and hydrodynamic escape can remove up to a few Earth masses of hydrogen and helium from these planets, potentially turning them into “habitable evaporated cores” (HECs). This process is most likely for mini-Neptunes with solid cores on the order of 1 M and up to about 50% H/He by mass, and can occur around all M dwarfs, particularly close to the IHZ. HECs are less likely to form around K and G dwarfs because of these stars' shorter superluminous PMS phases and shorter XUV saturation timescales. Furthermore, we find that HECs cannot form from mini-Neptunes with core masses greater than about 2 M and more than a few percent H/He by mass; thus, massive terrestrial super-Earths currently in the HZs of M dwarfs have probably always been terrestrial. Our results are thus similar to those of Lammer et al. (2014), who showed that planets more massive than ∼1.5 M typically cannot lose their accreted nebular gas in the HZs of solar-type stars.
Whether or not a given mini-Neptune forms a HEC is highly dependent on the early XUV evolution of the host star. In particular, a long XUV saturation timescale (tsat≳1 Gyr) is needed to fully evaporate the envelopes of mini-Neptunes in the HZs of early and mid M dwarfs. While a large tsat is consistent with the long activity timescales (West et al., 2008) and long spin-down times (Pizzolato et al., 2000) of M dwarfs, more observations are needed to pin down tsat for these stars. Moreover, the relative strength of the X-ray and EUV luminosity early on also affects whether or not HECs can form, as this determines whether the atmospheric escape is energy-limited or radiation/recombination-limited. In the energy-limited regime, which occurs for an X-ray-dominated flow, the escape is fast, and HECs can form throughout most of the HZs of all M dwarfs. In the radiation/recombination-limited regime, which applies to an EUV-dominated flow, HECs only form close to the inner edge of low-mass M dwarfs.
We further find that HECs can form from planets on circular as well as eccentric orbits, though they are marginally more likely to have higher e in the long run. While there exist feedbacks between atmospheric mass loss and tidal evolution, we find that these are only significant at e≳0.5; in these cases, whether or not a HEC forms can depend just as strongly on the tidal properties of the planet as on the efficiency of the atmospheric escape.
Many of the Earth-mass terrestrial planets detected in the HZs of M dwarfs in the coming years could be HECs. These planets should have abundant surface water and are likely to be water worlds, whose potential for habitability should be investigated further.

Acknowledgments

We wish to thank Russell Deitrick, Eric Agol, and the rest of the VPL team for their insightful suggestions and many rich discussions. We also wish to thank two anonymous referees for their excellent comments and suggestions. This work was supported by the NASA Astrobiology Institute's Virtual Planet Laboratory under Cooperative Agreement solicitation NNA13AA93A and by a generous fellowship from the ARCS Seattle chapter.

Footnotes

1
Tidal locking refers to the state in which a body's rotation rate is fixed by tidal forces at an equilibrium value. While for circular orbits this implies ωi=n, in the general case of an eccentric orbit in the CPL model, the planet's rotation rate assumes a slightly super-synchronous value. See Barnes et al. (2013) for a discussion.
2
As we explain below, it is quite possible that HECs are all but terrestrial, since it is likely that they have large ice mass fractions and are compositionally distinct from Earth. By “terrestrial,” in this case, we mean planets that are similar in size and mass to Earth and are not gaseous; these may or may not have a rocky surface like Earth. Our definition of “terrestrial” therefore encompasses water worlds.
3
While we refer to (43) as the “exact” expression for the (inverse of the) mass loss enhancement factor, it is important to remember that it is the eccentric version of the third-order expression derived by Erkaev et al. (2007) and is, in this sense, still an approximation to the true enhancement.

Abbreviations

CHZ, center of the habitable zone.
CPL, constant phase lag.
CTL, constant time lag.
EUV, extreme ultraviolet.
HECs, habitable evaporated cores.
HZ, habitable zone.
IHZ, inner edge of the habitable zone.
OHZ, outer edge of the habitable zone.
PMS, pre-main sequence.
RLO, Roche lobe overflow.
XUV, X-ray/extreme ultraviolet.

References

Abe Y., Abe-Ouchi A., Sleep N.H., and Zahnle K.J. (2011) Habitable zone limits for dry planets. Astrobiology 11:443–460.
Alexander M.E. (1973) The weak friction approximation and tidal evolution in close binary systems. Astrophys Space Sci 23:459–510.
Alibert Y. (2014) On the radius of habitable planets. Astron Astrophys 561.
Anglada-Escudé G., Tuomi M., Gerlach E., Barnes R., Heller R., Jenkins J.S., Wende S., Vogt S.S., Butler R.P., Reiners A., and Jones H.R.A. (2013) A dynamically packed planetary system around GJ 667C with three super-Earths in its habitable zone. Astron Astrophys 556:A126
Baraffe I., Chabrier G., Allard F., and Hauschildt P.H. (1998) Evolutionary models for solar metallicity low-mass stars: mass-magnitude relationships and color-magnitude diagrams. Astron Astrophys 337:403–412.
Barnes R., Mullins K., Goldblatt C., Meadows V.S., Kasting J.F., and Heller R. (2013) Tidal Venuses: triggering a climate catastrophe via tidal heating. Astrobiology 13:225–250.
Batygin K. and Morbidelli A. (2013) Dissipative divergence of resonant orbits. Astron J 145.
Bodenheimer P. and Lissauer J.J. (2014) Accretion and evolution of ∼2.5 M planets with voluminous H/He envelopes. Astrophys J 791.
Bolmont E., Raymond S.N., Leconte J., and Matt S.P. (2012) Effect of the stellar spin history on the tidal evolution of close-in planets. Astron Astrophys 544.
Borucki W.J., Agol E., Fressin F., Kaltenegger L., Rowe J., Isaacson H., Fischer D., Batalha N., Lissauer J.J., Marcy G.W., Fabrycky D., Désert J.M., Bryson S.T., Barclay T., Bastien F., Boss A., Brugamyer E., Buchhave L.A., Burke C., Caldwell D.A., Carter J., Charbonneau D., Crepp J.R., Christensen-Dalsgaard J., Christiansen J.L., Ciardi D., Cochran W.D., DeVore E., Doyle L., Dupree A.K., Endl M., Everett M.E., Ford E.B., Fortney J., Gautier T.N. III, Geary J.C., Gould A., Haas M., Henze C., Howard A.W., Howell S.B., Huber D., Jenkins J.M., Kjeldsen H., Kolbl R., Kolodziejczak J., Latham D.W., Lee B.L., Lopez E., Mullally F., Orosz J.A., Prsa A., Quintana E.V., Sanchis-Ojeda R., Sasselov D., Seader S., Shporer A., Steffen J.H., Still M., Tenenbaum P., Thompson S.E., Torres G., Twicken J.D., Welsh W.F., and Winn J.N. (2013) Kepler-62: a five-planet system with planets of 1.4 and 1.6 Earth radii in the habitable zone. Science 340:587–590.
Carpenter J.M., Mamajek E.E., Hillenbrand L.A., and Meyer M.R. (2006) Evidence for mass-dependent circumstellar disk evolution in the 5 Myr old Upper Scorpius OB association. Astrophys J 651:L49–L52.
Carter J.A., Agol E., Chaplin W.J., Basu S., Bedding T.R., Buchhave L.A., Christensen-Dalsgaard J., Deck K.M., Elsworth Y., Fabrycky D.C., Ford E.B., Fortney J.J., Hale S.J., Handberg R., Hekker S., Holman M.J., Huber D., Karoff C., Kawaler S.D., Kjeldsen H., Lissauer J.J., Lopez E.D., Lund M.N., Lundkvist M., Metcalfe T.S., Miglio A., Rogers L.A., Stello D., Borucki W.J., Bryson S., Christiansen J.L., Cochran W.D., Geary J.C., Gilliland R.L., Haas M.R., Hall J., Howard A.W., Jenkins J.M., Klaus T., Koch D.G., Latham D.W., MacQueen P.J., Sasselov D., Steffen J.H., Twicken J.D., and Winn J.N. (2012) Kepler-36: a pair of planets with neighboring orbits and dissimilar densities. Science 337:556–559.
Chabrier G. and Baraffe I. (1997) Structure and evolution of low-mass stars. Astron Astrophys 327:1039–1053.
Chang S.-H., Gu P.-G., and Bodenheimer P. H. (2010) Tidal and magnetic interactions between a hot Jupiter and its host star in the magnetospheric cavity of a protoplanetary disk. Astrophys J 708:1692–1702.
Chassefière E. (1996) Hydrodynamic escape of hydrogen from a hot water-rich atmosphere: the case of Venus. J Geophys Res 101.
Choukroun M. and Grasset O. (2007) Thermodynamic model for water and high-pressure ices up to 2.2 GPa and down to the metastable domain. J Chem Phys 127.
Cook B.A., Williams P.K.G., and Berger E. (2014) Trends in ultracool dwarf magnetism. II. The inverse correlation between X-ray activity and rotation as evidence for a bimodal dynamo. Astrophys J 785.
Correia A.C.M. and Laskar J. (2011) Tidal evolution of exoplanets. In Exoplanets, edited by Seager S., University of Arizona Press, Tucson, AZ, pp 239–266.
Darwin G.H. (1880) On the secular changes in the elements of the orbit of a satellite revolving about a tidally distorted planet. Royal Society of London Philosophical Transactions Series I 171:713–891.
Durney B.R., De Young D.S., and Roxburgh I.W. (1993) On the generation of the large-scale and turbulent magnetic fields in solar-type stars. Solar Physics 145:207–225.
Engle S.G. and Guinan E.F. (2011) Red dwarf stars: ages, rotation, magnetic dynamo activity and the habitability of hosted planets. In 9th Pacific Rim Conference on Stellar Astrophysics, ASP Conference Series Vol. 451, edited by Qain S., Leung K., Zhu L., and Kwok S., Astronomical Society of the Pacific, San Francisco, p 285.
Erkaev N.V., Kulikov Y.N., Lammer H., Selsis F., Langmayr D., Jaritz G.F., and Biernat H.K. (2007) Roche lobe effects on the atmospheric loss from “hot Jupiters.” Astron Astrophys 472:329–334.
Erkaev N.V., Lammer H., Odert P., Kulikov Y.N., Kislyakova K.G., Khodachenko M.L., Güdel M., Hanslmeier A., and Biernat H. (2013) XUV-exposed, non-hydrostatic hydrogen-rich upper atmospheres of terrestrial planets. Part I: atmospheric expansion and thermal escape. Astrobiology 13:1011–1029.
Fabrycky D.C., Ford E.B., Steffen J.H., Rowe J.F., Carter J.A., Moorhead A.V., Batalha N.M., Borucki W.J., Bryson S., Buchhave L.A., Christiansen J.L., Ciardi D.R., Cochran W.D., Endl M., Fanelli M.N., Fischer D., Fressin F., Geary J., Haas M.R., Hall J.R., Holman M.J., Jenkins J.M., Koch D.G., Latham D.W., Li J., Lissauer J.J., Lucas P., Marcy G.W., Mazeh T., McCauliff S., Quinn S., Ragozzine D., Sasselov D., and Shporer A. (2012) Transit timing observations from Kepler. IV. Confirmation of four multiple-planet systems by simple physical models. Astrophys J 750.
Ferraz-Mello S., Rodríguez A., and Hussmann H. (2008) Tidal friction in close-in satellites and exoplanets: the Darwin theory re-visited. Celestial Mechanics and Dynamical Astronomy 101:171–201.
Ford E.B. and Rasio F.A. (2008) Origins of eccentric extrasolar planets: testing the planet-planet scattering model. Astrophys J 686:621–636.
Fortney J.J., Marley M.S., and Barnes J.W. (2007) Planetary radii across five orders of magnitude in mass and stellar insolation: application to transits. Astrophys J 659:1661–1672.
Goldreich P. (1966) Final spin states of planets and satellites. Astron J 71:1–7.
Goldreich P. and Soter S. (1966) Q in the Solar System. Icarus 5:375–389.
Greenberg R. (2009) Frequency dependence of tidal Q. Astrophys J 698:L42–L45.
Gu P.-G., Lin D.N.C., and Bodenheimer P.H. (2003) The effect of tidal inflation instability on the mass and dynamical evolution of extrasolar planets with ultrashort periods. Astrophys J 588:509–534.
Hart M.H. (1979) Habitable zones about main sequence stars. Icarus 37:351–357.
Hashimoto G.L., Abe Y., and Sugita S. (2007) The chemical composition of the early terrestrial atmosphere: formation of a reducing atmosphere from CI-like material. J Geophys Res Planets 112.
Hayashi C. (1961) Stellar evolution in early phases of gravitational contraction. Publ Astron Soc Jpn Nihon Tenmon Gakkai 13:450–452.
Hayashi C., Nakazawa K., and Nakagawa Y. (1985) Formation of the Solar System. In Protostars and Planets II, edited by Black D.C. and Matthews M.S., University of Arizona Press, Tucson, AZ, pp 1100–1153.
Heller R., Leconte J., and Barnes R. (2011) Tidal obliquity evolution of potentially habitable planets. Astron Astrophys 528:A27
Howard A.W., Marcy G.W., Bryson S.T., Jenkins J.M., Rowe J.F., Batalha N.M., Borucki W.J., Koch D.G., Dunham E.W., Gautier T.N. III, Van Cleve J., Cochran W.D., Latham D.W., Lissauer J.J., Torres G., Brown T.M., Gilliland R.L., Buchhave L.A., Caldwell D.A., Christensen-Dalsgaard J., Ciardi D., Fressin F., Haas M.R., Howell S.B., Kjeldsen H., Seager S., Rogers L., Sasselov D.D., Steffen J.H., Basri G.S., Charbonneau D., Christiansen J., Clarke B., Dupree A., Fabrycky D.C., Fischer D.A., Ford E.B., Fortney J.J., Tarter J., Girouard F.R., Holman M.J., Johnson J.A., Klaus T.C., Machalek P., Moorhead A.V., Morehead R.C., Ragozzine D., Tenenbaum P., Twicken J.D., Quinn S.N., Isaacson H., Shporer A., Lucas P.W., Walkowicz L.M., Welsh W.F., Boss A., Devore E., Gould A., Smith J.C., Morris R.L., Prsa A., Morton T.D., Still M., Thompson S.E., Mullally F., Endl M., and MacQueen P.J. (2012) Planet occurrence within 0.25 AU of solar-type stars from Kepler. Astrophys J 201.
Howell S.B., Sobeck C., Haas M., Still M., Barclay T., Mullally F., Troeltzsch J., Aigrain S., Bryson S.T., Caldwell D., Chaplin W.J., Cochran W.D., Huber D., Marcy G.W., Miglio A., Najita J.R., Smith M., Twicken J.D., and Fortney J.J. (2014) The K2 mission: characterization and early results. Publications of the Astronomical Society of the Pacific 126:398–408.
Hunten D.M. (1982) Thermal and nonthermal escape mechanisms for terrestrial bodies. Planet Space Sci 30:773–783.
Ida S. and Lin D.N.C. (2008a) Toward a deterministic model of planetary formation. IV. Effects of Type I migration. Astrophys J 673:487–501.
Ida S. and Lin D.N.C. (2008b) Toward a deterministic model of planetary formation. V. Accumulation near the ice line and super-Earths. Astrophys J 685:584–595.
Jackson A.P., Davis T.A., and Wheatley P.J. (2012) The coronal X-ray-age relation and its implications for the evaporation of exoplanets. Mon Not R Astron Soc 422:2024–2043.
Jackson B., Greenberg R., and Barnes R. (2008) Tidal evolution of close-in extrasolar planets. Astrophys J 678:1396–1406.
Jackson B., Miller N., Barnes R., Raymond S.N., Fortney J.J., and Greenberg R. (2010) The roles of tidal evolution and evaporative mass loss in the origin of CoRoT-7 b. Mon Not R Astron Soc 407:910–922.
Jeans J. (1925) The Dynamical Theory of Gases, 4th ed., Cambridge University Press, Cambridge, UK.
Kaltenegger L., Sasselov D., and Rugheimer S. (2013) Water planets in the habitable zone: atmospheric chemistry, observable features, and the case of Kepler-62e and -62f. Astrophys J 775:L47
Kasting J.F. and Pollack J.B. (1983) Loss of water from Venus. I—Hydrodynamic escape of hydrogen. Icarus 53:479–508.
Kasting J.F., Whitmire D.P., and Reynolds R.T. (1993) Habitable zones around main sequence stars. Icarus 101:108–128.
Kislyakova K.G., Lammer H., Holmström M., Panchenko M., Odert P., Erkaev N.V., Leitzinger M., Khodachenko M.L., Kulikov Y.N., Güdel M., and Hanslmeier A. (2013) XUV-exposed, non-hydrostatic hydrogen-rich upper atmospheres of terrestrial planets. Part II: hydrogen coronae and ion escape. Astrobiology 13:1030–1048.
Kislyakova K.G., Johnstone C.P., Odert P., Erkaev N.V., Lammer H., Lüftinger T., Holmström M., Khodachenko M.L., and Güdel M. (2014) Stellar wind interaction and pick-up ion escape of the Kepler-11 “super-Earths.” Astron Astrophys 562.
Kopparapu R.K., Ramirez R., Kasting J.F., Eymet V., Robinson T.D., Mahadevan S., Terrien R.C., Domagal-Goldman S., Meadows V., and Deshpande R. (2013) Habitable zones around main-sequence stars: new estimates. Astrophys J 765.
Koskinen T.T., Harris M.J., Yelle R.V., and Lavvas P. (2013) The escape of heavy atoms from the ionosphere of HD209458b. I. A photochemical-dynamical model of the thermosphere. Icarus 226:1678–1694.
Kulow J.R., France K., Linsky J., and Loyd R.O.P. (2014) Lyα transit spectroscopy and the neutral hydrogen tail of the hot Neptune GJ 436b. Astrophys J 786:132
Lammer H., Lichtenegger H.I.M., Biernat H.K., Erkaev N.V., Arshukova I.L., Kolb C., Gunell H., Lukyanov A., Holmstrom M., Barabash S., Zhang T.L., and Baumjohann W. (2006) Loss of hydrogen and oxygen from the upper atmosphere of Venus. Planet Space Sci 54:1445–1456.
Lammer H., Lichtenegger H.I.M., Kulikov Y.N., Grießmeier J.-M., Terada N., Erkaev N.V., Biernat H.K., Khodachenko M.L., Ribas I., Penz T., and Selsis F. (2007) Coronal mass ejection (CME) activity of low-mass M stars as an important factor for the habitability of terrestrial exoplanets. II. CME-induced ion pick up of Earth-like exoplanets in close-in habitable zones. Astrobiology 7:185–207.
Lammer H., Kislyakova K.G., Odert P., Leitzinger M., Schwarz R., Pilat-Lohinger E., Kulikov Y.N., Khodachenko M.L., Güdel M., and Hanslmeier A. (2011) Pathways to Earth-like atmospheres. Extreme ultraviolet (EUV)-powered escape of hydrogen-rich protoatmospheres. Orig Life Evol Biosph 41:503–522.
Lammer H., Erkaev N.V., Odert P., Kislyakova K.G., Leitzinger M., and Khodachenko M.L. (2013) Probing the blow-off criteria of hydrogen-rich ‘super-Earths.’ Mon Not R Astron Soc 430:1247–1256.
Lammer H., Stökl A., Erkaev N.V., Dorfi E.A., Odert P., Güdel M., Kulikov Y.N., Kislyakova K.G., and Leitzinger M. (2014) Origin and loss of nebula-captured hydrogen envelopes from ‘sub’- to ‘super-Earths’ in the habitable zone of Sun-like stars. Mon Not R Astron Soc 439:3225–3238.
Lecavelier des Etangs A., Ehrenreich D., Vidal-Madjar A., Ballester G.E., Désert J.-M., Ferlet R., Hébrard G., Sing D.K., Tchakoumegni K.-O., and Udry S. (2010) Evaporation of the planet HD 189733b observed in H I Lyman-α. Astron Astrophys 514.
Leconte J., Chabrier G., Baraffe I., and Levrard B. (2010) Is tidal heating sufficient to explain bloated exoplanets? Consistent calculations accounting for finite initial eccentricity. Astron Astrophys 516:A64
Léger A., Selsis F., Sotin C., Guillot T., Despois D., Mawet D., Ollivier M., Labèque A., Valette C., Brachet F., Chazelas B., and Lammer H. (2004) A new family of planets? “Ocean-Planets.” Icarus 169:499–504.
Leitzinger M., Odert P., Kulikov Y.N., Lammer H., Wuchterl G., Penz T., Guarcello M.G., Micela G., Khodachenko M.L., Weingrill J., Hanslmeier A., Biernat H.K., and Schneider J. (2011) Could CoRoT-7b and Kepler-10b be remnants of evaporated gas or ice giants? Planet Space Sci 59:1472–1481.
Levi A., Sasselov D., and Podolak M. (2013) Volatile transport inside super-Earths by entrapment in the water-ice matrix. Astrophys J 769.
Lichtenegger H.I.M., Lammer H., Grießmeier J.-M., Kulikov Y.N., von Paris P., Hausleitner W., Krauss S., and Rauer H. (2010) Aeronomical evidence for higher CO2 levels during Earth's Hadean epoch. Icarus 210:1–7.
Lissauer J.J. (2007) Planets formed in habitable zones of M dwarf stars probably are deficient in volatiles. Astrophys J 660:L149–L152.
Lithwick Y. and Wu Y. (2012) Resonant repulsion of Kepler planet pairs. Astrophys J 756.
Lopez E.D. and Fortney J.J. (2013) The role of core mass in controlling evaporation: the Kepler radius distribution and the Kepler-36 density dichotomy. Astrophys J 776.
Lopez E.D. and Fortney J.J. (2014) Understanding the mass-radius relation for sub-Neptunes: radius as a proxy for composition. Astrophys J 792.
Lopez E.D., Fortney J.J., and Miller N. (2012) How thermal evolution and mass-loss sculpt populations of super-Earths and sub-Neptunes: application to the Kepler-11 system and beyond. Astrophys J 761.
Luger R. and Barnes R. (2014) Extreme water loss and abiotic O2 buildup on planets throughout the habitable zones of M dwarfs. Astrobiology 15.
Mordasini C., Alibert Y., Klahr H., and Henning T. (2012a) Characterization of exoplanets from their formation. I. Models of combined planet formation and evolution. Astron Astrophys 547.
Mordasini C., Alibert Y., Georgy C., Dittkrist K.-M., Klahr H., and Henning T. (2012b) Characterization of exoplanets from their formation. II. The planetary mass-radius relationship. Astron Astrophys 547.
Murray-Clay R.A., Chiang E.I., and Murray N. (2009) Atmospheric escape from hot Jupiters. Astrophys J 693:23–42.
Ogihara M. and Ida S. (2009) N-Body simulations of planetary accretion around M dwarf stars. Astrophys J 699:824–838.
Öpik E.J. (1963) Selective escape of gases. Geophysical Journal International 7:490–506.
Owen J.E. and Jackson A.P. (2012) Planetary evaporation by UV and X-ray radiation: basic hydrodynamics. Mon Not R Astron Soc 425:2931–2947.
Parker E.N. (1955) Hydromagnetic dynamo models. Astrophys J 122:293–314.
Parker E.N. (1965) Dynamical theory of the solar wind. Space Sci Rev 4:666–708.
Pascucci I., Apai D., Luhman K., Henning T., Bouwman J., Meyer M.R., Lahuis F., and Natta A. (2009) The different evolution of gas and dust in disks around Sun-like and cool stars. Astrophys J 696:143–159.
Penev K., Jackson B., Spada F., and Thom N. (2012) Constraining tidal dissipation in stars from the destruction rates of exoplanets. Astrophys J 751.
Pizzolato N., Maggio A., Micela G., Sciortino S., Ventura P., and D'Antona F. (2000) Determination of convective turnover times in young stars. In Stellar Clusters and Associations: Convection, Rotation, and Dynamos, Astronomical Society of the Pacific Conference Series, Vol. 198, edited by Pallavicini R., Micela G., and Sciortino S., Astronomical Society of the Pacific, San Francisco, p 71.
Pollack J.B., Hubickyj O., Bodenheimer P., Lissauer J.J., Podolak M., and Greenzweig Y. (1996) Formation of the giant planets by concurrent accretion of solids and gas. Icarus 124:62–85.
Quintana E.V., Barclay T., Raymond S.N., Rowe J.F., Bolmont E., Caldwell D.A., Howell S.B., Kane S.R., Huber D., Crepp J.R., Lissauer J.J., Ciardi D.R., Coughlin J.L., Everett M.E., Henze C.E., Horch E., Isaacson H., Ford E.B., Adams F.C., Still M., Hunter R.C., Quarles B., and Selsis F. (2014) An Earth-sized planet in the habitable zone of a cool star. Science 344:277–280.
Rafikov R.R. (2006) Atmospheres of protoplanetary cores: critical mass for nucleated instability. Astrophys J 648:666–682.
Raymond S.N. and Cossou C. (2014) No universal minimum-mass extrasolar nebula: evidence against in situ accretion of systems of hot super-Earths. Mon Not R Astron Soc 440:L11–L15.
Raymond S.N., Scalo J., and Meadows V.S. (2007) A decreased probability of habitable planet formation around low-mass stars. Astrophys J 669:606–614.
Reid I.N. and Hawley S.L. (2005) New Light on Dark Stars: Red Dwarfs, Low-Mass Stars, Brown Dwarfs, 2nd ed., Springer, Berlin.
Ribas Á., Merín B., Bouy H., and Maud L.T. (2014) Disk evolution in the solar neighborhood. I. Disk frequencies from 1 to 100 Myr. Astron Astrophys 561.
Ribas I., Guinan E.F., Güdel M., and Audard M. (2005) Evolution of the solar activity over time and effects on planetary atmospheres. I. High-energy irradiances (1–1700 Å). Astrophys J 622:680–694.
Ricker G.R., Latham D.W., Vanderspek R.K., Ennico K.A., Bakos G., Brown T.M., Burgasser A.J., Charbonneau D., Clampin M., Deming L.D., Doty J.P., Dunham E.W., Elliot J.L., Holman M.J., Ida S., Jenkins J.M., Jernigan J.G., Kawai N., Laughlin G.P., Lissauer J.J., Martel F., Sasselov D.D., Schingler R.H., Seager S., Torres G., Udry S., Villasenor J.N., Winn J.N., and Worden S.P. (2010) Transiting Exoplanet Survey Satellite (TESS) [abstract 450.06]. In 215th Meeting of the American Astronomical Society, Vol. 42, Bulletin of the American Astronomical Society, Washington, DC.
Ritter H. (1988) Turning on and off mass transfer in cataclysmic binaries. Astron Astrophys 202:93–100.
Rogers L.A. (2014) Most 1.6 Earth-radius planets are not rocky. arXiv:1407.4457v1.
Rogers L.A., Bodenheimer P., Lissauer J.J., and Seager S. (2011) Formation and structure of low-density exo-Neptunes. Astrophys J 738.
Scalo J., Kaltenegger L., Segura A.G., Fridlund M., Ribas I., Kulikov Y.N., Grenfell J.L., Rauer H., Odert P., Leitzinger M., Selsis F., Khodachenko M.L., Eiroa C., Kasting J., and Lammer H. (2007) M stars as targets for terrestrial exoplanet searches and biosignature detection. Astrobiology 7:85–166.
Schlaufman K.C., Lin D.N.C., and Ida S. (2010) A population of very hot super-Earths in multiple-planet systems should be uncovered by Kepler. Astrophys J 724:L53–L58.
Segura A., Walkowicz L.M., Meadows V., Kasting J., and Hawley S. (2010) The effect of a strong stellar flare on the atmospheric chemistry of an Earth-like planet orbiting an M dwarf. Astrobiology 10:751–771.
Selsis F., Chazelas B., Bordé P., Ollivier M., Brachet F., Decaudin M., Bouchy F., Ehrenreich D., Grießmeier J.-M., Lammer H., Sotin C., Grasset O., Moutou C., Barge P., Deleuil M., Mawet D., Despois D., Kasting J.F., and Léger A. (2007) Could we identify hot ocean-planets with CoRoT, Kepler and Doppler velocimetry? Icarus 191:453–468.
Sepinsky J.F., Willems B., and Kalogera V. (2007) Equipotential surfaces and Lagrangian points in nonsynchronous, eccentric binary and planetary systems. Astrophys J 660:1624–1635.
Sepinsky J.F., Willems B., Kalogera V., and Rasio F.A. (2009) Interacting binaries with eccentric orbits. II. Secular orbital evolution due to non-conservative mass transfer. Astrophys J 702:1387–1392.
Shematovich V.I., Ionov D.E., and Lammer H. (2014) Heating efficiency in hydrogen-dominated upper atmospheres. arXiv:1409.0730v1.
Skumanich A. (1972) Time scales for Ca II emission decay, rotational braking, and lithium depletion. Astrophys J 171:565–567.
Stauffer J.R., Caillault J.-P., Gagne M., Prosser C.F., and Hartmann L.W. (1994) A deep imaging survey of the Pleiades with ROSAT. Astrophys J Suppl Ser 91:625–657.
Stelzer B., Marino A., Micela G., López-Santiago J., and Liefke C. (2013) The UV and X-ray activity of the M dwarfs within 10 pc of the Sun. Mon Not R Astron Soc 431:2063–2079.
Strom K.M., Strom S.E., Edwards S., Cabrit S., and Skrutskie M.F. (1989) Circumstellar material associated with solar-type pre-main-sequence stars—a possible constraint on the timescale for planet building. Astron J 97:1451–1470.
Swift J.J., Johnson J.A., Morton T.D., Crepp J.R., Montet B.T., Fabrycky D.C., and Muirhead P.S. (2013) Characterizing the Cool KOIs. IV. Kepler-32 as a prototype for the formation of compact planetary systems throughout the Galaxy. Astrophys J 764.
Terquem C. and Papaloizou J.C.B. (2007) Migration and the formation of systems of hot Super-Earths and Neptunes. Astrophys J 654:1110–1120.
Tian F. (2009) Thermal escape from super Earth atmospheres in the habitable zones of M stars. Astrophys J 703:905–909.
Tian F., Toon O.B., Pavlov A.A., and De Sterck H. (2005) A hydrogen-rich early Earth atmosphere. Science 308:1014–1017.
Tian F., Kasting J.F., Liu H.-L., and Roble R.G. (2008) Hydrodynamic planetary thermosphere model: 1. Response of the Earth's thermosphere to extreme solar EUV conditions and the significance of adiabatic cooling. J Geophys Res Planets 113.
Trilling D.E., Benz W., Guillot T., Lunine J.I., Hubbard W.B., and Burrows A. (1998) Orbital evolution and migration of giant planets: modeling extrasolar planets. Astrophys J 500:428–439.
Tuomi M., Jones H.R.A., Barnes J.R., Anglada-Escudé G., and Jenkins J.S. (2014) Bayesian search for low-mass planets around nearby M dwarfs—estimates for occurrence rate based on global detectability statistics. Mon Not R Astron Soc 441:1545–1569.
Vidal-Madjar A., Lecavelier des Etangs A., Désert J.-M., Ballester G.E., Ferlet R., Hébrard G., and Mayor M. (2003) An extended upper atmosphere around the extrasolar planet HD209458b. Nature 422:143–146.
Walter F.M., Brown A., Mathieu R.D., Myers P.C., and Vrba F.J. (1988) X-ray sources in regions of star formation. III—Naked T Tauri stars associated with the Taurus-Auriga complex. Astron J 96:297–325.
Ward W.R. (1997) Protoplanet migration by nebula tides. Icarus 126:261–281.
Watson A.J., Donahue T.M., and Walker J.C.G. (1981) The dynamics of a rapidly escaping atmosphere: applications to the evolution of Earth and Venus. Icarus 48:150–166.
West A.A., Hawley S.L., Bochanski J.J., Covey K.R., Reid I.N., Dhital S., Hilton E.J., and Masuda M. (2008) Constraining the age-activity relation for cool stars: the Sloan Digital Sky Survey Data Release 5 low-mass star spectroscopic sample. Astron J 135:785–795.
Williams D.M. and Pollard D. (2002) Earth-like worlds on eccentric orbits: excursions beyond the habitable zone. International Journal of Astrobiology 1:61–69.
Wolfgang A. and Lopez E. (2014) How rocky are they? The composition distribution of Kepler's sub-Neptune planet candidates within 0.15 AU. arXiv:1409.2982v1.
Wright N.J., Drake J.J., Mamajek E.E., and Henry G.W. (2011) The stellar-activity-rotation relationship and the evolution of stellar dynamos. Astrophys J 743:48
Yang J., Cowan N.B., and Abbot D.S. (2013) Stabilizing cloud feedback dramatically expands the habitable zone of tidally locked planets. Astrophys J 771.

Appendix

Appendix A. Mass Loss for Eccentric Orbits

A.1. The quasi-static approximation

In Section 2.3.6, we argued that for planets on eccentric orbits, the Roche lobe radius may be calculated by replacing a with the instantaneous planet-star distance r in (8), leading to tidally enhanced mass loss close to pericenter. However, this is true only if the orbital period is large compared to the dynamical timescale of the planet, thus allowing sufficient time for the atmosphere to assume the equilibrium shape due to the changing potential. Sepinsky et al. (2007) called this limit the quasi-static approximation. Below we show that this approximation is valid for most of our planets.
As in Sepinsky et al. (2007), we begin by introducing the parameter
\begin{align*}\alpha ( e , \,f ) = \frac { ( 1 + e ) ^ { 1 / 2 } } { ( 1 - e ) ^ { 3 / 2 } } \times \mid 1 - f \mid \tag { 27 } \end{align*}
where
\begin{align*}f = \frac { P_ { { \rm orb } } \, ( 1 - e ) ^ { 3 / 2 } } { P_ { { \rm rot } } \, ( 1 + e ) ^ { 1 / 2 } } \tag { 28 } \end{align*}
is the ratio of the rotational angular velocity to the orbital angular velocity at periastron. Sepinsky et al. (2007) show that provided the condition
\begin{align*} \frac { P_ { { \rm orb } } } { \tau_ { \rm dyn } } > > \alpha ( e , \,f ) \tag { 29 } \end{align*}
holds, the system may be treated as quasi-static. The orbital period is given by Kepler's law,
\begin{align*}P_ { { \rm orb } } = 2 \pi \sqrt { \frac { a^3 } { GM_\star } } \tag { 30 } \end{align*}
while τdyn, the dynamical timescale of the planet, is
\begin{align*}\tau_ { { \rm dyn } } \sim \frac { 1 } { \sqrt { GM_p / R_p^3 } } \tag { 31 } \end{align*}
The ratio of the two will be smallest for close-in, low-mass planets with large radii. The minimum value in our runs occurs for a super-inflated 2 M, 30 R planet in the IHZ of a 0.08 M M dwarf, for which Porb/τdyn≈2.4.
We must now compare this to α(e, f). The equilibrium rotational period for a synchronously rotating planet is obtained from (65) and (78) in the CPL and CTL models, respectively, so that we may write α(e) as
\begin{align*} \alpha ( e , \,f ) = \frac { ( 1 + e ) ^ { 1 / 2 } } { ( 1 - e ) ^ { 3 / 2 } } \times \mid 1 - \frac { ( 1 - e ) ^ { 3 / 2 } \,\omega_ { 2 , { \rm eq } } } { ( 1 + e ) ^ { 1 / 2 } \ \,n } \mid \tag { 32 } \end{align*}
This equation is plotted in Fig. A1. Note that α only begins to approach Porb/τdyn for e≳0.6 in the CPL model and e≳0.4 in the CTL model, and only for the most inflated planets in the IHZ. We therefore urge caution in interpreting results above e≳0.5, where the quasi-static approximation may not hold for some planets.

A.2. The mass loss enhancement factor

Since the flux FXUV is not constant over the course of one orbit, (5) must be modified slightly:
\begin{align*} \dot { M } ( t ) = \frac { R_ { { \rm XUV } } ^3 \epsilon_ { { \rm XUV } } L_ { { \rm XUV } } } { 4 GM_p } \left[ r ( t ) ^2 \left( 1 - \frac { 3 } { 2 \xi ( t ) } + \frac { 1 } { 2 \xi ( t ) ^3 } \right) \right] ^ { - 1 } \tag { 33 } \end{align*}
where r(t) is the instantaneous separation between the centers of mass of the star and the planet and where we have plugged in for FXUV in terms of LXUV and made use of (6). The parameter ξ must also be modified, as the value of the Roche radius will also change as the planet's distance from the star changes during one orbit:
\begin{align*} \begin{split} \xi ( t ) & \equiv \frac { R_ { { \rm Roche } } } { R_ { { \rm XUV } } } \\ & = \left[ \left( \frac { M_p } { 3M_\star } \right) ^ { 1 / 3 } \frac { 1 } { R_ { { \rm XUV } } } \right] r ( t ) = \frac { r ( t ) } { a } \xi \end{split} \tag { 34 } \end{align*}
FIG. A1. The Sepinsky et al. (2007) parameter α as a function of e for a synchronously rotating planet in the CPL model (black) and in the CTL model (gray). The dashed line corresponds to the minimum value of the ratio Porb/τdyn across all our runs. Note that for e≲0.5, the quasi-static approximation is probably valid.
To avoid confusion, henceforth ξ will denote the original expression for circular orbits (Eqs. 7 and 8), while ξ(t) is the time-dependent parameter given by the expression above. The time average of \(\dot{M}\) over one orbit is
\begin{align*}\langle \dot { M } \rangle_t = \frac { 1 } { 2 \pi / n } \int_0^ { 2 \pi / n } \dot { M } ( t ) dt \tag { 35 } \end{align*}
where n is the mean motion. Now, we know that the relationship between n and the eccentric anomaly E is
\begin{align*}nt = E - e \sin E \tag{36}\end{align*}
and that
\begin{align*}r ( E ) = a ( 1 - e \cos E ) \tag{37}\end{align*}
so
\begin{align*} \frac { dt } { dE } = \frac { 1 } { n } ( 1 - e \cos E ) = \frac { r ( E ) } { an } \tag { 38 } \end{align*}
Substituting into (35), we have
\begin{align*} \begin{split} \langle \dot { M } \rangle_t = \frac { 1 } { 2 \pi / n } \int_0^ { 2 \pi } \dot { M } ( E ) \frac { dt } { dE } dE \\ = \frac { 1 } { 2 \pi a } \int_0^ { 2 \pi } \dot { M } ( E ) r ( E ) dE \end{split} \tag { 39 } \end{align*}
Now, introducing the mass loss rate for e=0 and Ktide=1 from (5),
\begin{align*}\dot { M } _0 \equiv \frac { R^3_ { { \rm XUV } } \epsilon_ { { \rm XUV } } L_ { { \rm XUV } } } { 4 GM_p a^2 } \tag { 40 } \end{align*}
it follows from (33) and (34) that
\begin{align*}\dot { M } ( E ) = \,\dot { M } _0 a^2 \left[ r ( E ) ^2 \left( 1 - \frac { 3 } { 2 \xi } \left( \frac { a } { r ( E ) } \right) + \frac { 1 } { 2 \xi^3 } \left( \frac { a } { r ( E ) } \right) ^3 \right) \right] ^ { - 1 } \tag { 41 } \end{align*}
Plugging this into (39) and using (37),
\begin{align*} \langle \dot { M } \rangle_t & = \frac { \dot { M } _0 } { 2 \pi } a \int_0^ { 2 \pi } \left( r ( E ) - \frac { 3 } { 2 \xi } a + \frac { 1 } { 2 } \left( \frac { a } { \xi } \right) ^3 r ( E ) ^ { - 2 } \right) ^ { - 1 } dE \\ & = \frac { \dot { M } _0 } { 2 \pi } \int_0^ { 2 \pi } \left[ ( 1 - e \cos E ) - \frac { 3 } { 2 \xi } + \frac { 1 } { 2 \xi^3 ( 1 - e \cos E ) ^2 } \right] ^ { - 1 } dE \\ & = \frac { \dot { M } _0 } { K_ { { \rm ecc } } } \tag { 42 } \end{align*}
where we define the eccentric mass loss enhancement factor
\begin{align*} \begin{split} { 1 / K_{{ \rm ecc}} \equiv \frac {1} {2 \pi} \int_0^{2 \pi} \Bigg [ ( 1 - e \cos E ) - \frac {3} (2 \xi} + \frac {1} {2 \xi^3 ( 1 - e \cos E ) ^2} \Bigg ] ^{ - 1} dE \end{split} \tag{43} \end{align*}
Note that for e=0, Kecc reduces to the circular version of Ktide, (6). Unfortunately, the integral (42) is not analytic. However, if the last term in the integral is small compared to the first two, an analytic solution is possible. In particular, if we write
\begin{align*} \frac { 1 } { 2 \xi^3 ( 1 - e \cos E ) ^2 } = \eta \ \left[ ( 1 - e \cos E ) - \frac { 3 } { 2 \xi } \right] \tag { 44 } \end{align*}
then we may neglect the last term provided η <<1. Since the term is largest at pericenter (E=0), we may instead require that
\begin{align*}0 < \frac { 1 } { 2 \xi^3 ( 1 - e ) ^2 ( 1 - \frac { 3 } { 2 \xi } - e ) } < < 1 \tag { 45 } \end{align*}
If this holds, (42) simplifies to
\begin{align*} \langle \dot { M } \rangle_t & \approx \frac { \dot { M } _0 } { 2 \pi } \int_0^ { 2 \pi } \left( 1 - \frac { 3 } { 2 \xi } - e \cos E \right) ^ { - 1 } dE \\ & \approx \frac { \dot { M } _0 } { 2 \pi ( 1 - \frac { 3 } { 2 \xi } ) } \int_0^ { 2 \pi } \left[1 - \left( \frac { e } { 1 - \frac { 3 } { 2 \xi } } \right) \cos E \right] ^ { - 1 } dE \\ & \approx \frac { \dot { M } _0 } { 1 - \frac { 3 } { 2 \xi } } \frac { 1 } { \sqrt { 1 - \left( \frac { e } { 1 - \frac { 3 } { 2 \xi } } \right) ^2 } } \\ & \approx \frac { \dot { M } _0 } { \left(1 - \frac { 3 } { \xi } - \frac { 9 } { 4 \xi^2 } - e^2 \right) ^ { 1 / 2 } } \tag { 46 } \end{align*}
Thus, provided condition (45) holds (typically for ξ≳10), we may write
\begin{align*}K_ { { \rm ecc } } \approx \sqrt { 1 - \frac { 3 } { \xi } - \frac { 9 } { 4 \xi^2 } - e^2 } \tag { 47 } \end{align*}
We note that for e=0 and keeping only first-order terms in ξ, the above expression agrees with that of Erkaev et al. (2007) for ξ >>1.
In Fig. A2 we plot the loss rate enhancement factor of Erkaev et al. (2013), 1/Ktide (blue lines), along with the eccentric version, 1/Kecc (exact3 expression in black, approximate expression in red). Recall that there are two distinct effects contributing to the extra mass loss for eccentric orbits: the \(1 / \sqrt{1 - e^2}\) flux enhancement and the smaller Roche lobe distance during part of the orbit. In order to better distinguish between these two, we display both 1/Ktide (dotted blue lines) and \(1 / ( K_{{ \rm tide}} \sqrt{1 - e^2}\) (dashed blue lines). Thus, any difference between 1/Ktide and 1/Kecc in the figure is due solely to the Roche lobe effect.
FIG. A2. Loss rate enhancement factor 1/K as a function of the normalized Roche lobe distance ξ=RRoche/Rp for e=0, 0.25, 0.5, and 0.75. The Erkaev et al. (2007) model (Ktide) is plotted in blue, with the flux enhancement factor \(1 / \sqrt{1 - e^2}\) (dashed) and without it (dotted). The exact expression Kecc derived in this work (42) is plotted in black, and the analytic approximation (46) is plotted in red. For low values of ξ, condition (45) is not satisfied, and the approximation diverges significantly. However, at high eccentricity, as ξ increases, the analytic approximation converges to the exact value of Kecc sooner than Ktide. Finally, we note that for sufficiently large ξ, all curves accounting for the flux enhancement converge to the same value.
We first compare the curves corresponding to the exact expressions (blue and black). For e=0, Kecc=Ktide, as expected. As the eccentricity increases to 0.25, the flux enhancement effect remains small (i.e., the dashed and dotted curves are similar), while the Roche lobe effect begins to become significant (the dashed black curve exceeds the dashed blue curve), particularly for low values of ξ. For e≳0.5, the effect becomes even more important, leading to an enhancement of a factor of several for ξ≲10.
A particularly important effect is that high eccentricities will cause the planet to undergo RLO at larger values of ξ. In the circular Erkaev et al. (2007) model, RLO occurs when ξ=1 by definition. However, it is straightforward to show that for
\begin{align*}\xi = \xi_ { { \rm crit } } = \frac { 1 } { 1 - e } \tag { 48 } \end{align*}
the expression in the integral (43) diverges, resulting in an infinite mass loss rate. This occurs because at pericenter, where r=a(1−e), ξ(r)=1. In other words, for ξξcrit, the planet will overflow its Roche lobe during at least part of its orbit, leading to rapid mass loss.
It is important to note that the approximate expression (dashed red curve) converges only for ξ≳10. For smaller values of ξ, it overestimates the mass loss enhancement significantly and should therefore not be used. However, for high eccentricities (see the last panel, for instance), it converges to 1/Kecc quicker than 1/Ktide. Provided condition (45) is satisfied, the approximate analytic expression (46) does a good job at modeling the actual mass loss.

A.3. Orbits crossing the critical radius

As we described in Section 2.3.5, the escape regime for an EUV-dominated flow may switch from energy-limited to radiation/recombination-limited above a certain critical value of the XUV flux, Fcrit, given by
\begin{align*}F_ { { \rm crit } } = \left( \frac { B } { A } \right)^ { 2 } \tag { 49 } \end{align*}
where
\begin{align*}A \equiv \frac { \pi \epsilon_ { { \rm XUV } } R_ { { \rm XUV } } ^3 } { GM_p K_ { { \rm tide } } } \tag { 50 } \end{align*}
and
\begin{align*}B \equiv 7.11 \times 10^7 { \rm g } ^ { \frac { 1 } { 2 } } { \rm s } ^ { \frac { 1 } { 2 } } \left( \frac { R_ { \rm p } } { R_ \oplus } \right) ^ { \frac { 3 } { 2 } } \tag { 51 } \end{align*}
are the coefficients multiplying the flux in the energy-limited and radiation/recombination-limited equations, respectively (Eqs. 5 and 13). At fixed XUV luminosity, this corresponds to a certain critical orbital radius,
\begin{align*}r_ { { \rm crit } } = \sqrt { \frac { L_ { { \rm XUV } } } { 4 \pi F_ { { \rm crit } } } } \tag { 52 } \end{align*}
Planets on low-eccentricity orbits that do not cross rcrit are therefore safely within either the energy-limited (a >>rcrit) or radiation/recombination-limited (a <<rcrit) regime. In this case, the orbit-averaged mass loss rate is determined simply by replacing F in (5) and (13) by its orbit-averaged value, 〈F〉, given by replacing Lbol with LXUV in (2).
However, when an orbit is sufficiently eccentric that the atmospheric escape regime switches from energy-limited to radiation/recombination-limited over the course of a single orbit, the method outlined above is no longer rigorously correct. We must instead integrate the two mass loss rate expressions over the portions of the orbit where they apply. Fortunately, as we demonstrate below, these integrals are analytic. This allows us to calculate the time-averaged value of the mass loss rate, \(\langle \dot{M} \rangle_t\), much as we did in Section 2.3.6:
\begin{align*}\langle \dot { M } \rangle_t = \frac { A } { P } \int_ { { \rm EL } } F ( t ) dt + \frac { B } { P } \int_ { { \rm RRL } } F^ \frac { 1 } { 2 } ( t ) dt \tag { 53 } \end{align*}
where P is the orbital period, F(t) is the instantaneous XUV flux, and EL and RRL correspond to the energy-limited and radiation/recombination-limited portions of the orbit, respectively. By noting as in (38) that
\begin{align*}dt = \frac { r ( E ) } { an } dE = \frac { r ( E ) P } { 2 \pi a } dE \tag { 54 } \end{align*}
where E is the eccentric anomaly, we may write
\begin{align*} \langle \dot {M} \rangle_t =& \frac{B} {2 \pi a} \int_0^{E_ {\rm crit}} F^{\frac {1}{2}} (E) r (E) dE \\& + \frac {A} {2 \pi a} \int_ {E_ {\rm crit}}^{2 \pi - E_ {\rm crit}} F ( E ) r ( E ) dE \\& + \frac {B} {2 \pi a} \int_ {2 \pi - E_ {{\rm crit}}} ^ {2 \pi} F^{\frac {1} {2}} (E) r (E) dE \tag {55} \end{align*}
where
\begin{align*}E_ { \rm crit } = \cos^ { - 1 } \left( \frac { 1 } { e } - \frac { r_ { { \rm crit } } } { ae } \right) \tag { 56 } \end{align*}
is the value of the eccentric anomaly when r=rcrit. The three integrals above follow from the fact that starting from pericenter (E=0), the planet is (by assumption) in the RRL regime up until E=Ecrit, switches to EL until E=2πEcrit, and completes the orbit in the RRL regime. By symmetry of the orbit, the first and last integrals are identical, so we may simplify:
\begin{align*} \langle \dot { M } \rangle_t & = \frac { A } { 2 \pi a } \int_ { E_ { \rm crit } } ^ { 2 \pi - E_ { \rm crit } } F ( E ) r ( E ) dE \\& \quad + \frac { B } { \pi a } \int_0^ { E_ { \rm crit } } F^ \frac { 1 } { 2 } ( E ) r ( E ) dE \tag {57} \end{align*}
Now, noting that F(E)=LXUV/4πr2(E) and r(E)=a(1−e cos E), we have
\begin{align*} \langle \dot {M} \rangle_t & = \frac{AL_ {{\rm XUV}}} {8 \pi^2 a^2} \int_ {E_ {\rm crit}} ^ {2 \pi - E_ {{\rm crit}}} \frac {dE} {1 - e \cos E} + \frac {B} {\pi a} \sqrt {\frac {L_ {{\rm XUV}}} {4 \pi}} \int_0^ {E_ {\rm crit}} dE \\& = \frac {AL_ {{\rm XUV}}} {8 \pi^2 a^2} \int_ {E_ {\rm crit}} ^ {2 \pi - E_ {\rm crit}} \frac {dE} {1 - e \cos E} + \frac {BE_ {\rm crit}} {\pi a} \sqrt {\frac {L_ {{\rm XUV}}} {4 \pi}} \tag{58}\end{align*}
By evaluating the integral, we may finally write
\[\openup4pt\begin{align*} \begin{split} \langle \dot{M} \rangle_t = & A\,\langle {F}\rangle \left[ 1 - {2 \over \pi} \tan^{ - 1} \left( { ( 1 + e ) \tan ( E_{ \rm crit} / 2) \over \sqrt{1 - e^2}} \right) \right] \\ & + B\,\langle {F}\rangle^{{1 \over 2}} \left[{ ( 1 - e^2 ) ^{{1 \over 4}}\ E_{ \rm crit} \over \pi} \right] \end{split} \tag{59}\end{align*}\]
where 〈F〉 is the orbit-averaged XUV flux. Note, importantly, that the expressions above are valid only for planets that cross the critical radius during their orbit; that is, planets for which the expression a(1 - e)<rcrit<a(1+e) holds. The mass loss rate for planets outside this region must be calculated via the method described at the beginning of the section.
Finally, note that the formalism derived here makes use of Ktide rather than its eccentric version Kecc, which we derived in Section A.2. It is in principle possible to account for the tidal enhancement during the energy-limited portion of the orbit, but the resulting integral would not be analytic. Moreover, the tidal enhancement due to the eccentricity is important primarily near pericenter, where the mass loss is radiation/recombination-limited and therefore independent of Ktide. Thus, while the method outlined above may underpredict the mass loss rate in some cases, the effect will be small.

Appendix B. Tidal Evolution Expressions: CPL

For zero inclination and zero obliquity, the tidal evolution equations for a, e, and ωi in a two-body system are (Barnes et al., 2013)
\begin{align*} \frac {da} {dt} &= \frac {a^2} {4GM_\star M_p} \mathop\sum_{i \neq j} Z_i^{\prime} \bigg(4 \epsilon_{0,i} + e^2 \bigg[ - 20 \epsilon_{0 , i} \\& + \frac {147} {2} \epsilon_{1 , i} + \frac {1} {2} \epsilon_{2,i} - 3 \epsilon_{5,i} \bigg]\bigg) \tag{60}\end{align*}
\begin{align*} \frac { de } { dt } = - \frac { ae } { 8GM_\star M_p } \mathop\sum_ { i \neq j } Z_i { \prime } \left( 2 \epsilon_ { 0 , i } - \frac { 49 } { 2 } \epsilon_ { 1 , i } + \frac { 1 } { 2 } \epsilon_ { 2 , i } + 3 \epsilon_ { 5 , i } \right) \tag { 61 } \end{align*}
\begin{align*} \frac { d \omega_i } { dt } = - \frac { Z_i { \prime } } { 8M_i r_ { g , i } ^2 \,R_i^2 n } \left( 4 \epsilon_ { 0 , i } + e^2 [ - 20 \epsilon_ { 0 , i } + 49 \epsilon_ { 1 , i } + \epsilon_ { 2 , i } ] \right) \tag { 62 } \end{align*}
where the sums are taken over the two bodies. Here, G is the gravitational constant, Ri is the radius of the ith body (planet or star), and rg,i are the radii of gyration. The parameter Z′i is defined as
\begin{align*}Z_i { \prime } \equiv 3 G^2 k_ { 2 , i } M_j^2 ( M_i + M_j ) \frac { R_i^5 } { a^9 } \frac { 1 } { nQ_i } \tag { 63 } \end{align*}
where k2,i is the Love number of degree 2 of the ith body, which is of order unity and quantifies the contribution of the tidal deformation to the total potential, n is the mean motion of the secondary body (i.e., the planet), and Qi are the tidal quality factors. The parameters εN,i are the signs of the phase lags (assumed equal in magnitude) of the Nth wave on the ith body, calculated from
\begin{align*} \epsilon_{0 , i} = { \rm sgn} ( 2 \omega_i - 2n ) \\ \epsilon_{1 , i} = { \rm sgn} ( 2 \omega_i - 3n ) \\ \epsilon_{2 , i} = { \rm sgn} ( 2 \omega_i - n ) \\ \epsilon_{5 , i} = { \rm sgn} ( n ) \\ \epsilon_{8 , i} = { \rm sgn} ( \omega_i - 2n ) \\ \epsilon_{9 , i} = { \rm sgn} ( \omega_i ) ( 64 ) \end{align*}
Since short-period planets around M dwarfs are likely to be tidally locked, one need not calculate the planetary spin from (62). Instead, the planet's rotation rate may be calculated from (Goldreich, 1966)
\begin{align*}\omega_{p , { \rm eq}}^{{ \rm CPL}} = n ( 1 + 9.5e^2 ) \tag{65}\end{align*}
where n is the mean motion.

B.1. The typical case

Because of the fast rotation of M dwarfs at early times, planets in the HZ are likely to be far outside the corotation radius of their parent stars. It is useful at this point to consider as an example the specific case of a tidally locked planet for which n <<ω. In this case, the stellar phase lags (64) are all positive, and εN, = +1. For a tidally locked planet, ε2,p=ε5,p=ε9,p=1, and ε1,p=ε8,p=−1. The parameter ε0,p, however, is less straightforward to calculate. If tidal locking is to be maintained, the average angular acceleration over one period must be zero. Ferraz-Mello et al. (2008) showed that the only self-consistent way to ensure this is if ε0,p has a different magnitude than the other lags, equal to
\begin{align*}\varepsilon_{0 , p} = 12 e^2 \varepsilon_{2 , p} \tag{66}\end{align*}
which in our notation corresponds to
\begin{align*}\epsilon_{0 , p} = + 12 e^2 \tag{67}\end{align*}
Note that if e=0, ε0,p=0, consistent with ωp=n in (64).
In the limit M >>Mp and keeping only terms up to order e2, Eqs. 60–62 now reduce to
\begin{align*} \frac{1}{a} \frac{da}{dt} = & \bigg[ 3 \sqrt{\frac{G}{M_\star}} \frac{k_{2 , \star} R_{\star}^5 M_p}{Q_\star} \bigg( 1 + \frac{51}{4} e^2 \bigg) \\& - 21 \sqrt{GM_\star^3} \frac{k_{2 , p} R_p^5}{Q_p M_p} e^2 \bigg] a^{- 13 / 2} \tag{68} \end{align*}
\begin{align*} \frac { 1 } { e } \frac { de } { dt } = \left[ \frac { 57 } { 8 } \sqrt { \frac { G } { M_\star } } \frac { k_ { 2 , \star } R_\star^5 M_p } { Q_\star } - \frac { 21 } { 2 } \sqrt { GM_\star^3 } \frac { k_ { 2 , p } R_p^5 } { Q_p M_p } \right] a^ { - 13 / 2 } \tag { 69 } \end{align*}
\begin{align*} \frac { d \omega_\star } { dt } = - \left[ \frac { 3 } { 32 } \frac { G^2 M_p^2 k_ { 2 , \star } R_\star^3 } { r_ { g , \star } ^2 n^2 Q_\star } \left( 1 + \frac { 15 } { 2 } e^2 \right) \right] a^ { - 9 } \tag { 70 } \end{align*}
We need not calculate p/dt, since we assume the planet's spin is instantaneously set to the equilibrium value. The first term in (68) and (69) is the orbital effect of the tide raised by the planet on the star. In both equations this term is positive, implying that the stellar tide acts to increase both a and e. This makes sense under the assumption that the planet is outside of corotation; the stellar bulge therefore leads the planet in the orbit, torquing the planet such that it speeds up and migrates outward. For an eccentric orbit, the strongest impulse occurs at pericenter; since the planet must return to that point in the orbit, the pericenter distance is preserved, but the faster orbital speed results in a more distant apocenter and thus higher eccentricity. A similar analysis for the tide raised by the star on the planet (the second term in each of the above equations) yields the result that such a tide acts to decrease both a and e.
The evolution of the orbit will therefore depend on the relative magnitudes of the stellar and planetary tides. A simple order-of-magnitude calculation will show that the ratio of the tide generated by the star on the planet to the tide generated by the planet on the star is
\begin{align*}\mid \frac { \dot { a } _ { \star \rightarrow p } } { \dot { a } _ { p \rightarrow \star } } \mid \approx 7 \left( \frac { M_\star } { M_p } \right) ^2 \left( \frac { R_p } { R_\star } \right) ^5 \left( \frac { Q_\star } { Q_p } \right) \left( \frac { e^2 } { 1 + \frac { 51 } { 4 } e^2 } \right) \tag { 71 } \end{align*}
For a 10 M, 2 R mini-Neptune with Qp=104 orbiting a 0.1 M, 0.15 R star with Q=106, this becomes
\begin{align*}\mid \frac { \dot { a } _ { \star \rightarrow p } } { \dot { a } _ { p \rightarrow \star } } \mid \approx 3 \times 10^4 \left( \frac { e^2 } { 1 + \frac { 51 } { 4 } e^2 } \right) \tag { 72 } \end{align*}
which is >>1 for all e≳0.01. For an Earth-like planet with Qp=102, the ratio is >>1 for all e≳0.001. Therefore, for planets starting with eccentricities ≳0.1, we would expect the planetary tide to dominate the evolution, such that the planet's orbit will shrink.

Appendix C. Tidal Evolution Expressions: CTL

The expressions for the evolutions of the orbital parameters are (Barnes et al., 2013)
\begin{align*} \frac { da } { dt } = \frac { 2a^2 } { GM_\star M_p } \mathop\sum_ { i \neq j } Z_i \left( \frac { f_2 ( e ) } { \beta^ { 12 } ( e ) } \frac { \omega_i } { n } - \frac { f_1 ( e ) } { \beta^ { 15 } ( e ) } \right) \tag { 73 } \end{align*}
\begin{align*} \frac { de } { dt } = \frac { 11ae } { 2GM_\star M_p } \mathop\sum_ { i \neq j } Z_i \left( \frac { f_4 ( e ) } { \beta^ { 10 } ( e ) } \frac { \omega_i } { n } - \frac { 18 } { 11 } \frac { f_3 ( e ) } { \beta^ { 13 } ( e ) } \right) \tag { 74 } \end{align*}
\begin{align*} \frac { d \omega_i } { dt } = \frac { Z_i } { M_i r_ { g , i } ^2 R_i^2 n } \left( \frac { f_2 ( e ) } { \beta^ { 12 } ( e ) } - \frac { f_5 ( e ) } { \beta^ { 9 } ( e ) } \frac { \omega_i } { n } \right) \tag { 75 } \end{align*}
where
\begin{align*}Z_i \equiv 3 G^2 k_ { 2 , i } M_j^2 ( M_i + M_j ) \frac { R^5_i } { a^9 } \tau_i \tag { 76 } \end{align*}
and
\begin{align*}&\beta (e) \equiv \sqrt {1 - e^2} \\& f_1 ( e ) \equiv 1 + \frac{31} {2} e^2 + \frac {255} {8 } e^4 + \frac {185} {16} e^6 + \frac {25} {64} e^8 \\& f_2 ( e ) \equiv 1 + \frac {15} {2} e^2 + \frac {45} {8} e^4 + \frac {5} {16} e^6 \\& f_3 ( e ) \equiv 1 + \frac {15} {4} e^2 + \frac {15} {8} e^4 + \frac {5} {64} e^6 \\& f_4 ( e ) \equiv 1 + \frac {3} {2} e^2 + \frac {1} {8 } e^4 \\& f_5 ( e ) \equiv 1 + 3 e^2 + \frac {3} {8} e^4 \tag{77}\end{align*}
As before, if we assume the planet is tidally locked, its rotation rate is given by
\begin{align*}\omega_ { p , { \rm eq } } ^ { \rm CTL } = n \left( \frac { f_2 ( e ) } { \beta^3 ( e ) f_5 ( e ) } \right)\tag { 78 } \end{align*}

C.1. The typical case

As in the CPL model, we expect that, for the typical case, the tides raised on the planet will dominate the evolution. Plugging the result of (78) into (73) and (74), we find that the second (negative) terms dominate, and the orbit should therefore shrink and circularize. Proceeding as before, we have
\begin{align*} \mid \frac { \dot { a } _ { \star \rightarrow p } } { \dot { a } _ { p \rightarrow \star } } \mid \approx \left( \frac { M_\star } { M_p } \right) ^2 \left( \frac { R_p } { R_\star } \right) ^5 \left( \frac { \tau_p } { \tau_\star } \right) \mid F \left(e , \frac { \omega_\star } { n } \right) \mid \tag { 79 } \end{align*}
where
\begin{align*} F \left(e , \frac{\omega_\star} {n} \right)& \equiv \frac {\frac {f_2^2 (e)} {f_5 (e) f_1 (e)} - 1} {\frac {f_2 ( e ) \beta^3 ( e )} {f_1 ( e )} \frac {\omega_\star} {n} - 1} \\& \approx - \frac {241} {16 \ ( \frac {\omega_\star} {n} - 1 )} e^2 \tag{80} \end{align*}
for small e. For the same mini-Neptune considered in the CPL case (72), this becomes
\begin{align*}\left| \frac { \dot { a } _ { \star \rightarrow p } } { \dot { a } _ { p \rightarrow \star } } \right| \approx 5 \times 10^5 \left( \frac { e^2 } { \frac { \omega_\star } { n } - 1 } \right)\tag { 81 } \end{align*}
which is >>1 for all e≳0.001\(\,\sqrt{ \omega_\star / n - 1}\), implying that the CTL model also predicts a net inward migration due to tides.

Appendix D. Rate of Change of the Flux

Conservation of angular momentum requires that, for tidal evolution, the rates of change of the eccentricity and the semimajor axis must be related through
\begin{align*} \frac { 1 } { e } \frac { de } { dt } = \frac { 1 } { 2a } \frac { da } { dt } \tag { 82 } \end{align*}
The rate of change of the flux due to orbital changes is
\begin{align*} \frac { dF } { dt } = & \frac { L_ { \rm bol } } { 4 \pi } \frac { d } { dt } \left( a^2 \sqrt { 1 - e^2 } \right) ^ { - 1 } \\ & = - \frac { L_ { \rm bol } } { 4 \pi } \left( a^2 \sqrt { 1 - e^2 } \right) ^ { - 2 } \Bigg [ 2a \frac { da } { dt } \sqrt { 1 - e^2 } \\& \quad + \left( \frac { a^2 } { 2 \sqrt { 1 - e^2 } } \right) \left( - 2e \frac { de } { dt } \right) \Bigg ] \\ & = \frac { L_ { \rm bol } } { 4 \pi } \frac { 1 } { a^2 } \left[ - \frac { 2 } { a } \frac { da } { dt } \frac { 1 } { ( 1 - e^2 ) ^ { 1 / 2 } } + \frac { e } { ( 1 - e^2 ) ^ { 3 / 2 } } \frac { de } { dt } \right] \\& = - F \frac { de } { dt } \left[ \frac { 4 } { e ( 1 - e^2 ) ^ { 1 / 2 } } - \frac { e } { ( 1 - e^2 ) ^ { 3 / 2 } } \right] \\& = - F \frac { de } { dt } \left[ \frac { 4 - 5e^2 } { e ( 1 - e^2 ) ^ { 3 / 2 } } \right] \tag { 83 } \end{align*}
For e≲0.7, this simplifies to
\begin{align*} \frac { dF } { dt } = - 4F \frac { 1 } { e } \frac { de } { dt } \tag { 84 } \end{align*}
which is equivalent to the trivial result
\begin{align*} \frac { 1 } { F } \frac { dF } { dt } = - \frac { 2 } { a } \frac { da } { dt } \tag { 85 } \end{align*}
Thus, at low eccentricities, the flux will always increase as the orbit shrinks. However, at very high eccentricities (e≳0.8), (83) predicts that dF/dt is negative when da/dt<0: the decrease in the flux due to the circularization of the orbit overpowers the increase due to the shrinking orbit.

Information & Authors

Information

Published In

cover image Astrobiology
Astrobiology
Volume 15Issue Number 1January 2015
Pages: 57 - 88
PubMed: 25590532

History

Published online: 15 January 2015
Published in print: January 2015
Accepted: 27 October 2014
Received: 22 August 2014

Permissions

Request permissions for this article.

Topics

Authors

Affiliations

R. Luger
Astronomy Department, University of Washington, Seattle, Washington.
NASA Astrobiology Institute – Virtual Planetary Laboratory Lead Team, USA.
R. Barnes
Astronomy Department, University of Washington, Seattle, Washington.
NASA Astrobiology Institute – Virtual Planetary Laboratory Lead Team, USA.
E. Lopez
Department of Astronomy and Astrophysics, University of California, Santa Cruz, California.
J. Fortney
Department of Astronomy and Astrophysics, University of California, Santa Cruz, California.
B. Jackson
Carnegie Institute of Washington, Washington, DC.
V. Meadows
Astronomy Department, University of Washington, Seattle, Washington.
NASA Astrobiology Institute – Virtual Planetary Laboratory Lead Team, USA.

Notes

Address correspondence to:Rodrigo LugerAstronomy DepartmentUniversity of WashingtonBox 351580Seattle, WA 98195-1580E-mail: [email protected]

Metrics & Citations

Metrics

Citations

Export citation

Select the format you want to export the citations of this publication.

View Options

View options

PDF/EPUB

View PDF/ePub

Get Access

Access content

To read the fulltext, please use one of the options below to sign in or purchase access.

Society Access

If you are a member of a society that has access to this content please log in via your society website and then return to this publication.

Restore your content access

Enter your email address to restore your content access:

Note: This functionality works only for purchases done as a guest. If you already have an account, log in to access the content to which you are entitled.

Media

Figures

Other

Tables

Share

Share

Copy the content Link

Share on social media

Back to Top