# Habitability of Terrestrial-Mass Planets in the HZ of M Dwarfs.

I. H/He-Dominated Atmospheres.

###### Abstract

The ubiquity of M dwarfs, combined with the relative ease of detecting terrestrial-mass planets around them, has made them prime targets for finding and characterising planets in the “Habitable Zone” (HZ). However, Kepler finds that terrestrial-mass exoplanets are often born with voluminous H/He envelopes, comprising mass-fractions () %. If these planets retain such envelopes over Gyr timescales, they will not be “habitable” even within the HZ. Given the strong X-ray/UV fluxes of M dwarfs, we study whether sufficient envelope mass can be photoevaporated away for these planets to become habitable. We improve upon previous work by using hydrodynamic models that account for radiative cooling as well as the transition from hydrodynamic to ballistic escape. Adopting a template active M dwarf XUV spectrum, including stellar evolution, and considering both evaporation and thermal evolution, we show that: (1) the mass-loss is (considerably) lower than previous estimates that use an “energy-limited” formalism and ignore the transition to Jeans escape; (2) at the inner edge of the HZ, planets with core mass M can lose enough H/He to become habitable if their initial envelope mass-fraction is 1%; (3) at the outer edge of the HZ, evaporation cannot remove a 1% H/He envelope even from cores down to 0.8 M. Thus, if planets form with bulky H/He envelopes, only those with low-mass cores may eventually be habitable. Cores 1 M, with 1% natal H/He envelopes, will not be habitable in the HZ of M dwarfs.

## 1 Introduction

M dwarfs comprise the bulk (75%) of the stellar population of our galaxy. Moreover, their low masses and small radii, compared to Sun-like stars, make planets of a given mass, size and orbital separation much easier to detect around them via both the transit and Doppler methods. In addition, the Habitable Zone (HZ) of such low luminosity red dwarfs lies considerably closer to the central star than in the case of solar-types; this further enhances the transit and Doppler signatures of HZ planets, and also allows a larger number of planetary orbits to be observed in a given time, making the detection and characterisation of such planets even easier. Consequently, the next generation of missions investigating exoplanets are aimed at later-type stars rather than the mainly solar-type ones targeted by Kepler, and discussions of exoplanet habitability increasingly focus on M dwarf systems.

There are many complications pertaining to the habitability of planets around M-dwarfs, even if they possess surface temperature and pressure conditions favourable to liquid water, such as tidal locking (e.g. Heller, Leconte, & Barnes, 2011; Leconte et al., 2015), run-away greenhouse effects (e.g. Kopparapu, 2013) and water loss (Luger et al., 2015b). The first major uncertainty, however, is whether they can possess such conditions at all, given what we currently know about exoplanets (e.g. Scalo et al., 2007).

The last decade of exoplanet discoveries has dramatically altered our understanding of a “standard” low-mass planet. Kepler has detected thousands of planetary candidates (e.g. Mullally et al., 2015), the majority of which are small ( R: Howard et al., 2012; Petigura et al., 2013; Morton & Swift, 2014; Silburt et al., 2015) and close to their host stars. Correcting for observational biases, one finds that most stars probably harbour at least one such “Kepler” planet (Fressin et al., 2013). Masses have been now been obtained for a significant number (albeit a small fraction) of these exoplanets, by combining Kepler transit data with either transit timing variations (TTVs; e.g., Wu & Lithwick, 2013; Hadden & Lithwick, 2014; Jontof-Hutter et al., 2014; Jontof-Hutter et al., 2015) or radial velocity measurements (e.g., Weiss & Marcy, 2014; Marcy et al., 2014; Dressing et al., 2015). A very surprising result to emerge from this work is that a large fraction of low-mass, close-in Kepler planets are enshrouded by voluminous H/He envelopes, which contain a non-negligible fraction of the planet’s mass (as demonstrated, for instance, by the extremely low densities of some: g cm; e.g., Wu & Lithwick, 2013; Jontof-Hutter et al., 2014; Jontof-Hutter et al., 2015): very unlike the low-mass planets in our own solar system.

In particular, statistical analyses of the properties of the low-mass exoplanet population reveal that the latter are inconsistent with a completely rocky composition (Rogers, 2015; Wolfgang et al., 2015). Instead, comparisons of the observed exoplanets with measured masses and radii to structural models (Wolfgang & Lopez, 2015) indicate that the dominant structure is a solid core overlaid with a 1% (by mass) H/He envelope. Since close-in planets are subject to intense irradiation by high energy stellar photons, which can drive an outflow by heating the upper layers of the atmosphere to close to the escape temperature (e.g. Lammer et al., 2003; Tian et al., 2005a; Murray-Clay et al., 2009; Owen & Jackson, 2012; Erkaev et al., 2016), the evaporation of such H/He envelopes becomes important for the planets’ evolution (e.g. Baraffe et al., 2005; Lopez et al., 2012). Forward modelling studies suggest that the majority of close-in exoplanets were born with H/He envelopes, but about half have subsequently lost these through evaporation (Owen & Wu, 2013; Lopez & Fortney, 2013). In summary, current exoplanet demographics lead us to infer that a dominant mode of low-mass planet formation produces a rocky Earth-like core surrounded by a H/He envelope with a mass-fraction of 1% (Wolfgang & Lopez, 2015). This observational result is supported by recent theoretical calculations within the framework of core accretion, which suggest that low-mass cores M will acquire envelope mass-fractions of order a few percent (Rogers et al., 2011; Bodenheimer & Lissauer, 2014; Lee et al., 2014; Lee & Chiang, 2015). For example, the scalings of Lee & Chiang (2015) imply that a 1 M core will accrete a H/He envelope mass-fraction of 0.1-1% during the gas disk’s lifetime.

We note that various studies of gas accretion onto Earth-like cores find that the results are sensitive to a number of uncertain parameters, such as the opacity and disc properties (e.g. Ikoma
& Hori, 2012; Lammer et al., 2014; Stökl, Dorfi,
& Lammer, 2015), and, especially, the assumed core accretion rate. Thus, these studies suggest a range of initial H/He envelope masses for a given core mass. However, this theoretical prejudice should not blind us to the empirical fact that current exoplanet data indicate that envelope mass fractions % are preferred^{1}^{1}1see Owen
& Wu (2016) for a discussion of why much larger envelopes may be scarce, as the data suggest..

For a terrestrial-mass planet, orbiting either a solar-type star or an M dwarf and located within the classical HZ (defined here as the range of orbital separations where an Earth-mass planet, with roughly Earth-like composition and atmosphere, can harbour liquid surface water; eg., see Kopparapu et al., 2013), a H/He envelope with mass-fraction of order a percent would preclude habitability, by yielding very high surface temperatures and pressures incompatible with liquid water. However, habitable conditions may be achieved if evaporation can reduce the H/He envelope mass-fraction to (e.g. Pierrehumbert & Gaidos, 2011); strip away this primordial envelope entirely so that it is replaced by a tenuous secondary atmosphere, such as those of the solar system terrestrial planets (e.g. Kasting et al., 1993; Kopparapu et al., 2013); or separate H and He in planets with a low initial H/He fraction to leave a habitable He atmosphere (Hu, Seager, & Yung, 2015). In other words, evaporation may turn a large population of uninhabitable low-mass planets with voluminous H/He envelopes into habitable ones. Such an effect is unlikely in the HZ of Sun-like stars, since their X-ray/UV fluxes over Gyr timescales at 1 AU are too low to remove a massive H/He envelope (though several works have shown it is possible to remove a small H/He atmosphere; e.g., Tian et al., 2005b; Erkaev et al., 2013; Lammer et al., 2014). M dwarfs, on the other hand, are far more active, and remain so for much longer, than solar-types, leading to much higher XUV fluxes within their HZs over Gyr timescales (e.g. Güdel, 2004; Jackson et al., 2012). Moreover, low-mass planets appear abundant around these red dwarfs, which are themselves the most common stars in the galaxy; evaporation, if efficient enough, could thus lead to a plethora of habitable planets in M dwarf systems.

In this work, therefore, we consider terrestrial-mass planets in the HZ of M dwarfs, with an initial H/He envelope mass-fraction of 1%, and investigate whether evaporation over a Gyr can remove a sufficient portion of this atmosphere to render the planet habitable at the end. We account for both stellar evolution and the thermal evolution of the planet. In order to make advances in this important area, we also abjure various simplifying – and, as it turns out, erroneous – assumptions about the evaporative flow made in previous studies of this subject, as described below, in order to obtain a more rigorous and realistic estimate of the mass loss.

## 2 Overview

### 2.1 Comparison to Previous Work

Given the Kepler results above, and the interest in habitable planets, several recent studies have investigated the evaporation of voluminous H/He envelopes around terrestrial-mass solid cores. However, previous studies have typically made two key simplifying assumptions: (i) radiative cooling in the flow is either neglected, or accounted for by assuming a fixed energy efficiency for driving the flow - either in a global or local sense; and (ii) evaporation is assumed to occur in the hydrodynamic limit at all times.

Both assumptions are likely to lead to an overestimation of the amount of H/He a planet can lose. First, radiative cooling becomes important when the timescale for a fluid element to advect its heat outwards (the ‘flow timescale’) becomes comparable to the timescale for radiative losses to cool that fluid element (e.g. Owen
& Alvarez, 2016). This lowers the temperature in the flow, which in turn pushes the sonic point to greater heights (lower densities) and thus reduces the mass-loss rate. Second, in order to be in the hydrodynamic limit, one requires that the gas remain collisional up to the sonic point (e.g., Owen
& Jackson, 2012). If this is not satisfied, the flow will collapse^{2}^{2}2There may be a small transition region where mass loss occurs sub-sonically (Tian et al., 2005a); however, these “breeze” solutions are exponentially sensitive to density and possibly unstable (Velli, 1994). and switch to Jeans escape, which has a much reduced mass-loss rate (as we demonstrate later, Jeans escape cannot remove a significant H/He atmosphere on Gyr time-scales, and other non-thermal processes are also unlikely to play a significant role). Finally, these two points are inter-related: as radiative cooling pushes the sonic point to higher heights and lower densities, it can also trigger the transition to Jeans escape earlier than in calculations that do not include cooling.

Lammer et al. (2014) and Johnstone et al. (2015) studied the evolution of terrestrial mass cores in the HZ of solar-type stars, using evaporation rates calculated in the “locally energy-limited” approach, where a fixed fraction of the absorbed photons’ energy is locally deposited into heat. This fraction is calibrated to detailed calculations by Shematovich et al. (2014) that explicitly solve the micro-physics of photon absorption; as such, it is more accurate than the standard “energy-limited” formalism, wherein an ad hoc fraction of the total incoming radiative flux goes into heating the gas. However, these calculations still neglect radiative cooling, and as such will overestimate the mass-loss rate when radiative cooling is important^{3}^{3}3Shematovich et al. (2014) caculate the fraction of the incoming XUV flux that goes into heating the gas (via Coulomb collisions of gas particles with photoelectrons liberated by the XUV), instead of being diverted into exciting, ionising or dissociating the atomic and molecular species. However, their calculation does not address radiative cooling; as such, this fraction does not equal the net heating rate, which is the difference of the heating and cooling rates, and the quantity of importance here.. Lammer et al. (2014) further argue that if the core accretion rate is very large (mass doubling times years), resulting in a high-entropy bloated atmosphere, evaporation may be able to completely remove initially massive H/He envelopes from low-mass cores and leave behind a potentially habitable planet. However, such large core accretion rates require that the time of planet formation to be fine-tuned to occur just before disc dispersal, otherwise the core would end up with a mass M. With more reasonable accretion rates, they find cores are likely to retain a large fraction of their original envelope. Johnstone et al. (2015), expanded on this work to demonstrate that Earth-mass planets at 1 AU around solar-type stars can lose a massive 1% H/He envelope over a Gyr only if the star is an unusually fast rotator (i.e., unusually active, with X-ray luminosity in the upper 90th percentile of the observed spread in X-ray luminosities in these stars). In the majority of cases, i.e., for solar-type stars with more standard rotation rates and actitivity levels, they find that such planets can only lose H/He envelopes with initial mass-fractions 0.1%.

Luger et al. (2015a) have recently investigated photoevaporation around M dwarfs. Using the standard “energy-limited” formalism and further assuming a hydrodynamic flow at all times, they argue that evaporation can completely strip H/He envelopes with a mass-fraction % from Earth-mass cores in the HZs of M dwarfs. If true, this would imply a potentially vast number of worlds with habitable conditions around M dwarfs. However, the problems noted above with the adopted assumptions call this result into question. Specifically, below an envelope mass-fraction of 1%, the radius of the planet is an extremely strong function of the remaining envelope mass-fraction (Lopez & Fortney, 2014); at the same time, the flow timescale increases strongly with decreasing planetary radius (as the atmosphere falls deeper into the planet’s gravity well). Thus, one expects radiative cooling to increase significantly as the H/He mass-fraction descends below 1%, resulting first in suppressed hydrodynamic mass-loss rates, and then a transition to Jeans escape with greatly diminished loss rates. This process can quench mass loss from planets with a hefty H/He envelope still remaining; indeed, the effect has already been identified at closer separations and higher planetary masses (Owen & Wu, 2013; Lopez & Fortney, 2013), and is the origin of the so-called “evaporation valley” (Jin et al., 2014).

In this paper, therefore, we eschew the assumptions made in earlier studies. Instead, we explicitly account for the effect described above, by: (i) smoothly transitioning from the regime where heating is balanced by outflow (as in the “locally energy-limited” calculations of Lammer et al., 2013) to the regime where heating is balanced by radiative cooling, using a parametrisation calibrated to detailed Monte-Carlo radiative transfer simulations; and (ii) using our hydrodynamic models to determine when the transition to Jeans escape occurs, triggering a strong suppression in mass-loss rates. In general, with our more appropriate treatment, we find that whether or not a planet can lose enough H/He to be considered habitable depends strongly on the core mass and orbital separation. In particular, the mass loss we derive for 1 M cores with initially 1% He/He envelopes in the HZ of M dwarfs is much lower than previous estimates, implying that such planets will not be habitable.

### 2.2 Model Outline

Stellar XUV-driven evaporation is particularly important around M dwarfs because these stars are extremely active: for instance, an average early- to mid-M dwarf remains at saturated levels of coronal activity, with / 10–10, for a few 100 Myr, compared to only a few10 Myr of saturated activity in an average solar-type star (e.g., Güdel, 2004)^{4}^{4}4These timescales can be extended, for both M dwarfs and solar-type stars, if the star is initially a very rapid rotator; e.g. Johnstone et al. (2015); Tu et
al. (2015). The point is that in general M dwarfs evince saturated activity for far longer than solar-types.. Consequently, a planet with a given equilibrium temperature receives orders of magnitude more high-energy radiation, integrated over its lifetime, around an M dwarf than around a star like the sun.

In this study, we are concerned with the evaporation of low-mass planets. Owen & Jackson (2012) demonstrated that at high X-ray fluxes, mass-loss from such planets is predominantly X-ray-driven, and EUV heating can be neglected. We thus confine ourselves to X-ray-driven flows here. To study the latter, we need a fiducial X-ray spectrum; we use that of AD Leo, an M3.5 dwarf of mass 0.4 M, and one of the most active nearby stars.

Th very strong coronal activity on AD Leo arises not because it is anomalous, but simply because it is relatively young compared to most nearby field stars. At a median age of 100 Myr (estimated age range 25–300 Myr; Shkolnik et al., 2009), it is still in its saturated phase of activity, with / = 10 ( = 710 erg s and = 910 erg s; Delfosse et al., 1998). Its high levels of activity, combined with its proximity ( pc), have made AD Leo a touchstone for understanding the coronae of M dwarfs.

Nevertheless, the intrinsic faintness of M dwarfs (AD Leo included: 10 L) means that, even when / is very high, the star is still very faint in X-rays, and it is not yet possible to acquire a full high-quality X-ray spectrum. Instead, one usually reconstructs the spectrum using the technique of Emission Measure Distributions (EMDs), based upon observations of bright emission lines. We use an EMD-derived synthetic X-ray spectrum of AD Leo supplied by J. Sanz-Forcada (pvt. comm. 2014). The methodology for constructing it is described in detail by Sanz-Forcada et al. (2011) and Chadney et al. (2015); this spectrum is also used by Chadney et al. (2015), who show that it agrees very well with current X-ray (and UV) data for AD Leo. We plot the spectrum in Fig. 1, and compare it to a scaled-solar spectrum (representing a young solar-type star); the comparison solar-like spectrum was calculated in Ercolano et al. (2008b) using the method of Kashyap & Drake (2000), and was designed to be representative of a young, saturated solar-type star. It is immediately clear that the AD Leo spectrum is considerably softer; thus, since X-ray heating is mainly due to photo-electrons liberated by soft 0.1-1 keV photons, we expect significantly more heating using our realistic M dwarf spectrum than with a scaled-solar proxy for it (e.g., as done in Owen & Jackson, 2012).

Since AD Leo is currently still in the saturated regime, we assume it has maintained the same / = 10 since birth to the present (the standard behaviour of saturated low-mass stars, before stellar rotation slows sufficiently for the activity to become unsaturated). However, its bolometric luminosity certainly has evolved with time, as the star has contracted towards the Main Sequence; we assume that its follows the theoretical evolutionary track of Baraffe et al. (1998) for a 0.4 M star, and scale accordingly in time. Finally, for ease of computation, we assume that AD Leo will remain saturated, i.e., its / will stay unchanged, till 1 Gyr (when we terminate our evaporation calculations). This is not physically strictly accurate, since a 0.4 M star should start spinning down after a few 100 Myr (e.g., Reiners & Mohanty, 2012), thus entering the unsaturated regime and evincing decreasing /; assuming saturation at these late times then formally implies that our mass-loss rates will be overestimated here. However, as we shall see, for the core-mass and H/He envelope mass-fraction regime we study in this paper, our calculations show that the planets are either completely stripped of their envelopes well before 100 Myr, or have entered the Jeans escape regime, with tiny mass-loss rates, by this age. As such, our simplifying assumption of saturated activity beyond a few 100 Myr has no discernible impact on our results.

In order to derive evaporative mass-loss rates, we must specify the planetary structure. As our starting point, we assume planets with Earth-like solid cores, composed of 2/3 rock + 1/3 iron, swathed in a H/He envelope with a mass-fraction (defined as /) up to 1%. As discussed above, this is consistent with both observations and theoretical calculations of the conditions at birth of “Kepler” planets (Wolfgang et al., 2015). Current theories of planet formation, though, do not strongly constrain the initial entropy (or equivalently, radius) of such planets, with both “hot start” planets (those with short initial cooling time scales) and “cold start” ones (with long initial cooling time scales) remaining viable. We thus choose initial planetary radii (entropies) corresponding to cooling timescales in the range 10–10 yr, to comfortably span the plausible range from “hot” to “cold start” scenarios (further discussed in §5).

Finally, we are interested in the potential habitability, under the effects of evaporation, of low-mass planets in the HZ of M dwarfs. The “classical” HZ is defined as the range of orbital separations around a star where an Earth-mass planet with a CO–HO–N atmosphere can sustain liquid water on the surface; the inner edge of this zone is (conservatively) set by the moist greenhouse effect, and the outer edge (again, conservatively) by the maximum greenhouse effect (Kasting et al., 1993; Kopparapu et al., 2013). The position of the HZ of course changes as the star evolves in temperature and luminosity; what we are really interested in is the fate of planets located within the stable HZ of an M dwarf on the Main Sequence (MS). Using the fitting equations supplied by Kopparapu et al. (2013) for the HZ boundaries, and the MS values of the stellar temperature and luminosity for a 0.4 M dwarf from the Baraffe et al. (1998) tracks ( 3500 K, 1.8610 L; the star arrives on the MS after 500 Myr), yields an inner edge of the classical HZ for our M dwarf at 0.15 AU, and an outer edge at 0.28 AU. For ease of calculation, and discussions for M dwarfs with masses around 0.4 M, we choose here to parametrise the HZ in terms of the blackbody temperature () of the planet instead, which we define as the equilibrium temperature of a planet with zero albedo, orbiting a star of luminosity at a radial separation of : [/. We set the inner edge of the classical HZ around our M dwarf to be at a MS blackbody temperature of K, and the outer edge at K. These imply radii of = 0.12 AU and 0.26 AU respectively on the MS: very close to the actual separations from the Kopparapu et al. (2013) equations. Since we are interested in framing the general theory of evaporation driven habitability here, the 10–20% difference has no appreciable effect on our conclusions. Detailed calculations of the habitability of systems with specific parameters (i.e. future detected exoplanet systems), can wait until a later date.

We discuss our theory of X-ray driven evaporation, in both the hydrodynamic and ballistic limits, in §3; our numerical implementation of this in §4; and our treatment of thermal evolution in §5. We present our results in §6, and explore the implications for habitability around M dwarfs in §7; our main conclusions are summarised in §8.

## 3 Evaporation: Theory

### 3.1 X-ray Heating

Our interest lies in low mass planets; in these, atmospheric heating and evaporation are dominated by X-rays (Owen & Jackson, 2012). It is well known that, in radiative equilibrium, heating due to high energy photons can be described using a relationship between temperature () and ionization parameter (, where is the received X-ray flux and the number density of particles) (Owen et al., 2010). Owen & Jackson (2012) used this principle to calculate solutions to the evaporation problem. However, the precise form of the relation is sensitive to the input X-ray spectrum and the metallicity of the gas. The X-ray spectrum of M dwarfs is quite different from that of solar type stars, being in general much softer in the keV range important for photo-electric heating. We must therefore calculate a new relation, specific to our template M dwarf, to understand the behaviour of the X-ray irradiated planetary atmosphere.

We use the radiative transfer code mocassin (Ercolano et al., 2003, 2005, 2008a) to solve for the relation appropriate to the AD Leo X-ray spectrum discussed in §2.2. We irradiate a plane parallel atmosphere, with densities spanning g cm, the expected range in the evaporative flow. For specificity, we scale the X-ray flux to that expected at a planet orbiting at K (0.12 AU); however, since we are concerned with mapping out the profile, rather than with the X-ray flux itself, the exact flux chosen to perform the calculation does not matter very much. The resulting profile is shown in Fig. 2, along with a functional fit which we adopt in our further calculations.

The soft M dwarf X-ray spectrum has a higher fraction of 0.1-1 keV photons compared to solar-type spectra (see Fig. 1); consequently, the M dwarf produces considerably higher gas temperatures at low ionization parameters (). We first attempt to solve for the evaporative flow by inserting our relation into the semi-analytical methodology of Owen & Jackson (2012), which assumes radiative equilibrium in the flow. This, however, leads to an inconsistency: over much of the parameter range of interest, the mechanical luminosity of the resulting evaporative wind becomes comparable to the radiative energy input rate, implying that radiative equilibrium is not a good approximation; instead, energy losses due to work may dominate over radiative losses much of the time. In hindsight, this is not too surprising, given the higher heating efficiency of our M dwarf X-ray spectrum relative to that in the original Owen & Jackson (2012) calculations. As such, we must derive a full numerical solution to the radiation-hydrodynamic problem. Nevertheless, the fact that, when radiative equilibrium holds, the gas temperature should equal that given by the relation enables us to simplify the numerics, fully span the pertinent parameter space and, crucially, smoothly transition to the situation where radiative losses do become important (which, as discussed in §2.1, ultimately determines when our hydrodynamic wind is quenched). Our specific numerical technique for accomplishing this is detailed in §4.1.

### 3.2 Ballistic Mass-loss Estimate

To assess the importance of non-hydrodynamic mass loss, we calculate the ballistic (Jeans escape) mass-loss rate, given by

(1) |

where is the planetary radius, the density at the exobase and the distribution of particle velocities. We assume that the exobase is in local thermodynamic equilibrium with the X-ray irradiation, so that its temperature is specified by our relation; is then the Maxwell-Boltzmann distribution at this temperature. The exobase density is

(2) |

where is the scale height of the atmosphere; the mean molecular weight of the gas particles, set to 1.35 for an atomic solar-abundance H/He mixture; and the collision cross-section, for which we adopt cm (Hunten, 1973; Tian et al., 2005a). Equation (1) can be solved numerically for the mass-loss rate as a function of core mass. Since we expect the Jeans escape rates to be low, we assume the extent of the H/He atmosphere is small, and approximate the planetary radius by that of the solid core. The core mass-radius relation is from Fortney et al. (2007a), for a core comprising 2/3 rock and 1/3 iron by mass.

The resulting are plotted in Fig 3 as a function of core mass and are compatible with those found previously for highly irradiated terrestrial mass planets (e.g. Tian et al., 2008; Erkaev et al., 2013). These mass-loss rates are indeed very low; in particular, they are significantly below the rates one might expect from hydrodynamic evaporation (as we will see). We find that Jeans escape can only remove a H/He envelope mass-fraction () of order 10 over a Gyr timescale: completely negligible compared to the 1% initial H/He mass-fractions we are concerned with here. Thus, in order to strip any significant H/He envelope from solid terrestrial-mass cores in the HZ of an M dwarf, efficient hydrodynamic evaporation appears essential.

## 4 Evaporation: Numerical Calculations

As discussed in §3.1, we must resort to full numerical hydrodynamic calculations to obtain the evaporation rates for the M dwarf case. Furthermore, to avoid the problems faced by Murray-Clay et al. (2009) in finding steady-state solutions to the wind problem, we explicitly integrate the time-dependent hydrodynamic equations until a steady-state is reached after several sound-crossing times (e.g Owen
& Alvarez, 2016): in all our calculations, we run the simulations for 30 crossing times to ensure that steady-state is achieved^{5}^{5}5This is confirmed by checking that a pseudo-Bernoulli function is conserved..

We employ the zeus hydrodynamical code (Stone & Norman, 1992; Hayes et al., 2006), which has been successfully tested and used for planetary evaporation problems in the past (e.g. Stone & Proga, 2009; Owen & Adams, 2014; Trammell, Li, & Arras, 2014; Owen & Alvarez, 2016). Our modified zeus algorithm includes heating by X-rays and their radiative transfer (c.f. Owen et al., 2010; Haworth et al., 2016).

### 4.1 X-ray Heating and Radiative Transfer

As noted in §3.1, our evaporative flow may be in a state far from radiative equilibrium. Specifically, we expect the gas to behave as follows. In steady-state, since the flow is expanding and we do not anticipate any heating by shocks, the gas will always be cooler than the temperature at radiative equilibrium (, given by the relation). When the gas is considerably cooler than - due to a large amount of work - the photo-heating rate will greatly exceed the radiative cooling rate, and we may ignore radiative cooling. As the gas temperature approaches , however, radiative cooling will become comparable to the photo-heating rate, and the temperature will asymptote to . With this qualitative picture in mind, we build a simplified thermal model for X-ray heating / radiative cooling, using a “Newtonian-heating” approach. We begin by defining a local heating time as , where is the gas internal energy per unit mass at , and is the local heating rate per unit mass due to absorption of high energy photons, given by

(3) |

where is the optical depth (from above) at a distance from the centre of the planet:

(4) |

Here, represents the fraction of the photon energy that is converted into heating rather than used up in ionization. We caution that this is not to be confused with the efficiency parameter commonly deployed in “energy-limited” estimates of the mass-loss rate (e.g. Lammer et al., 2003; Lammer et al., 2009; Erkaev et al., 2007; Lopez et al., 2012), which (arbitrarily) sets how much of the photon energy goes into work instead of being radiatively lost. In our analysis, merely accounts for the energy diverted to ionization processes; radiative energy losses, on the other hand, are explicitly accounted for by incorporating our derived relationship into the calculations (via in the energy equation, as described further below), and not via an ad hoc efficiency parameter .

Integrating a frequency-dependent radiation-hydrodynamics problem is beyond the scope of this work. Instead, we use the fact that the heating in our case is dominated by soft X-ray photons, mainly at energies of 1 keV (Ercolano et al., 2008b). We can thus approximate our frequency-dependent integral in Equation (3) by a monochromatic calculation at 1 keV, with an absorption cross-section of cm. Heating by 1 keV photons is primarily due to secondary photo-electrons; then depends on the ionization state of the gas, with a value of 0.1 in neutral gas and 0.9 in highly ionized plasma (Xu & McCray, 1991). While we do not attempt to solve for the full ionization structure of the gas, our mocassin results indicate the gas is partially ionized (e.g. Murray-Clay et al., 2009; Owen & Adams, 2014; Owen & Alvarez, 2016), so we pick an intermediate value of consistent with the ionization fraction. Our calculations should not be exceptionally sensitive to our choice, as the thermostatting of the flow towards the relationship will weaken the sensitivity of the flow to the choice of .

These considerations enable us to include radiative heating and cooling explicitly in our calculations, by writing the internal energy equation as

(5) |

where is the specific internal energy of a fluid element and is the ratio of specific heats, set to 5/3 for atomic gas. The first term on the R.H.S. of Equation (5) is the sink in internal energy due to cooling, while the second is the sum of the source and sink terms due to radiative heating and cooling respectively. When , this term is dominated by heating, with negligible cooling by radiation; when , radiative heating and cooling rates begin to balance each other (they must be equal when ); and finally, when , this term dominates over cooling and thermostats the gas temperature to .

Equation (5) allows us to satisfy energy conservation, so that we do not overestimate the fraction of incoming X-ray flux that is converted to mechanical luminosity, while also accounting for the burgeoning amount of energy lost as radiation as the gas temperature draws closer to the radiative equilibrium value. This approach is much more satisfactory than the typical assumption of a constant “heating efficiency” everywhere (e.g. Watson et al., 1981; Tian et al., 2005a; Erkaev et al., 2013; Lammer et al., 2014), since it allows regions to locally account for radiative losses and thus varying the net heating efficiency in both a local and global manner. Hence, in addition to a monochromatic ray-tracing calculation for 1 keV photons, we solve Equation (5) along with the time-dependent equations of hydrodynamics. We note that our model explicitly assumes that , i.e., that we approach thermal equilibrium from below. This is satisfied here since the flow is expanding (and hence cooling via work) everywhere. If this were not true – e.g. if there was heating due to compression or shocks – then Equation (5) would need to be generalised to include a cooling timescale, in addition to the heating timescale .

We solve our evaporation problem in 1D along a streamline connecting the star and the planet, within a spherical co-ordinate system (we thus implicitly assume that the velocity divergence is that of spherically symmetric outflow). We account for stellar gravity by using an effective gravitational potential of the form

(6) |

where and are the planetary and stellar masses respectively, is the radial distance from the centre of the planet, and the orbital semi-major axis. We do not however include a contribution from the Coriolis force, since it is negligible in setting the mass-loss rates (Stone & Proga, 2009; Murray-Clay et al., 2009; Owen & Jackson, 2012; Tripathi et al., 2015).

### 4.2 Hydrodynamic vs Ballistic Mass-loss

While our code solves the hydrodynamic problem over the entire input parameter space, the applicability of the equations of hydrodynamics must be verified a posteriori, once the steady-state flow solution is obtained. Whether a flow is hydrodynamic or ballistic depends on whether the gas particles are collisional or not on a length-scale shorter than the flow-scale (e.g. Johnson, Volkov, & Erwin, 2013). This property is characterised by the Knudsen number:

(7) |

where is the mean free path of the gas particles. The relevant flow length-scale for the problem is the pressure scale-height, given by

(8) |

Clearly, increases with distance from the planet (since the particle number density decreases as an exponential function of the pressure scale-height; note that the exobase is the radius at which ). At the same time, our hydrodynamic flow solutions are valid only if the gas remains in the hydrodynamic limit () all the way to the sonic surface; otherwise, the flow loses causal contact with the planet (and thus the ability to be influenced by upstream conditions) before the sonic point is reached (rendering the concept of a sonic point meaningless). Putting these two facts together, we follow Owen & Jackson (2012) in defining the evaporation to occur in the hydrodynamic regime if at the sonic point.

### 4.3 Mass-loss Rates

Using the framework described above, we construct a grid of evaporative mass-loss rates as a function of planetary mass, planetary radius and stellar X-ray luminosity, for a specified orbital distance. In order to comfortably cover the desired parameter space, we span planet masses ranging from 0.5 to 20 M; radii ranging from a value corresponding to a density of 10 g cm (for a given planet mass) to half a Hill radius at the given orbital separation; and X-ray luminosities (evenly spaced logarithmically) ranging from 10 to 10 erg s. This grid is then coupled to the planetary thermal evolution code, to finally yield the detailed evolution of a given planet and its atmosphere. Before proceeding to a discussion of thermal evolution (§5) and the final results for individual planets (§6), however, it is instructive to examine a slice through our mass-loss grid for a fiducial orbit and stellar X-ray luminosity.

Fig. 4 shows mass-loss rates as a function of planetary mass and radius, for an X-ray flux of 1.23 erg s: the irradiation expected from AD Leo at 100 Myr (roughly its current age) at an orbital separation of AU (corresponding to K, i.e., the inner edge of the HZ at 1 Gyr). Over-plotted is the mass-radius relationship for a solid core comprising 2/3 rock and 1/3 iron (dashed magenta line). A planet with a H/He envelope mass-fraction 1%, for which the radius is dominated by the solid core, will lie only slightly to the right of this line ( factor 2, Lopez & Fortney 2014), while a planet with a larger envelope mass-fraction, with radius dominated by the envelope, will lie much further to the right.

Also over-plotted are lines of constant Knudsen number at the sonic point (which we will denote by Kn): Kn = 10 (red dot-dash line), 1 (solid black line) and 0.1 (blue dotted line). For a fixed planetary mass, we see that Kn decreases with increasing radius. Thus (recalling that hydrodynamic escape requires Kn 1), it is easier to hydrodynamically evaporate a more distended atmosphere; an intuitive result, since more puffed-up atmospheres have a lower gravitational binding energy for a given planetary mass. Indeed, planets with radii to the right of the thin white line in the plot have atmospheres that are so loosely bound that they are unlikely to remain tethered to the planet for more than a Myr: our hydrodynamic calculations indicate that the evaporative flows on these planets are driven entirely by stellar bolometric heating rather than by X-ray heating, leading to a rapid “boil-off” phase during which most of the atmosphere is lost shortly after disc dispersal (see Owen & Wu, 2016, for a discussion of this process), see we are considering planets with relatively low envelope mass fractions this effect is only prominent in one case.

Finally, we also see that the lines of constant Kn are shallower than the mass-radius relationship for our solid core, with the Kn=1 line intersecting the latter locus at M. Since hydrodynamic flows can only occur to the right of the Kn=1 line (because they require Kn 1), only planets less massive than 2 M can be plausibly stripped entirely of their H/He envelopes via such flows. Higher mass planets can still undergo hydrodynamic evaporation as long as their envelopes remain very extended, but such flows will stall once the planetary radius shrinks sufficiently to hit the Kn=1 limit, which occurs while these planets still retain H/He mass-fractions 1%. Subsequent mass loss only happens ballistically, which can hardly put a dent in such atmospheres over Gyr timescales (as demonstrated in §3.2); thus, planets more massive than 2 M will retain substantial H/He envelopes over their lifetimes. In short, it is easier to completely strip a lower mass planet of its H/He envelope than a higher mass one: not surprising, since escape velocity increases with planetary mass.

The threshold mass of 2 M identified here, separating the envelope retention / stripping outcomes, is of course quantitatively valid only for the specific X-ray flux adopted in this snapshot plot, and only in the absence of thermal evolution (which is not accounted for in this plot). The true threshold will depend on the integrated history of the X-ray irradiation (since a higher X-ray flux, expected at earlier times, will push the sonic point deeper into the atmosphere, raising the threshold mass), as well as on the thermal evolution of the planet (since the atmosphere, contracting as it cools, sinks ever deeper into the planet’s gravity well, making hydrodynamic escape harder and decreasing the threshold mass). By coupling our full mass-loss grid (of which Fig. 4 is only an instantaneous slice at one X-ray flux) to our thermal evolution code, we will account for both of these effects in §5 and §6 below. Nevertheless, the physical effects summarised here will remain qualitatively valid, and only particularly low-mass cores can be stripped entirely by evaporation.

## 5 Thermal Evolution

In order to fully model the evolution in evaporative mass loss, we must account for how the H/He envelope cools and contracts. We solve the coupled thermal evolution and evaporation problem for a planet with a solid core and H/He envelope using the mesa stellar evolution code adapted to planets (Paxton et al., 2011, 2013). The method is described at length in Owen & Wu (2013); here we only summarise the salient inputs to the code.

As discussed in §2.2, the stellar bolometric luminosity () is assumed to follow the Lyon evolutionary track (Baraffe et al., 1998) for a 0.4 M star (appropriate for AD Leo), while the planet’s orbital distance is fixed at either AU or 0.26 AU, so that the blackbody temperature at the planet’s location at 1 Gyr is K or 200 K respectively (corresponding to the inner and outer edges of the classical HZ around an M dwarf; see §2.2). For reference, the evolution in blackbody temperature at 0.12 AU is shown in Fig. 5.

The radius of the solid core is assumed to be constant throughout the planet’s evolution, and is specified, for a given core mass and our adopted core composition of 2/3 rock + 1/3 iron, by the mass-radius relationship from Fortney et al. (2007a, b). We do account for the thermal content of the core, due to both radioactive decay and thermal heat capacity, by adopting an Earth-like value (see Lopez et al., 2012). The effects of both external bolometric irradiation by the star and internal release of thermal energy from the core are accounted for in mesa, using the relationship formulated by Guillot (2010) for a gray atmospheric model with both internal and external heat fluxes, where one uses an incoming stream with a frequency set to the peak of the black-body spectrum for the star and an outgoing IR stream in local thermodynamic equilibrium with the upper atmosphere.

Current planet formation models are unable to stringently constrain the initial thermal properties of newborn planets. Consequently, we consider planets with a large range of initial radii (or, more accurately, initial entropies). The latter are parametrised in terms of an initial cooling (i.e., Kelvin-Helmholz) timescale , defined as the ratio of a planet’s initial internal energy to its initial luminosity. Naively, one expects this Kelvin-Helmholz timescale to be of order the “age” of the planet at birth, i.e., the time it takes to form, and thus comparable to the mean lifetime of primordial disks, 10 yr (e.g. Lee et al., 2014). A ten-fold uncertainty in this estimate is not infeasible, due to variations in disc lifetimes (e.g. Hernández et al., 2007; Ercolano et al., 2011; Owen, Ercolano, & Clarke, 2011) or post-formation processes (e.g. Owen & Wu, 2016; Liu et al., 2015).; as such, we use starting models with in the range 10–10 yr, to encompass a plausible maximal difference between “hot start” planets (with short initial cooling times) and “cold start” ones (with long initial cooling times) (c.f. Marley et al., 2007; Spiegel & Burrows, 2012).

Finally, evaporation is included in our mesa calculations by coupling the code to our grid of newly calculated mass-loss rates as a function of planetary mass, radius and stellar X-ray flux (§4.3).

## 6 Planetary Evolution: Results

We present here our results for mass loss (accounting for both evaporation and thermal evolution) in the HZ of our fiducial M dwarf, for three core masses – = 0.8, 0.9 and 1 M – at two different orbital separations: AU (roughly the inner edge of the HZ) and 0.26 AU (roughly the outer edge). We investigate initial H/He mass-fractions (denoted henceforth by ) of 1% and 0.04% (i.e., and 410). We stop our calculations either after 1 Gyr of evolution, or once the envelope mass-fraction falls to , the amount that can be efficiently removed by Jeans escape on less than Gyr timescales (see §3.2).

### 6.1 Evolution Without Evaporation

We start by considering the evolution of the atmospheric temperature, density and pressure at the planetary surface (, and respectively), over Gyr timescales, in the absence of evaporation (i.e., with evaporation artificially turned off). These results constitute a baseline from which to judge whether evaporation will aid habitability or not.

We first carry out this analysis for = 300 K (i.e., at AU). Figs. 6, 7 and 8 show our results for = 0.8, 0.9 and 1.0 M respectively, for a H/He envelope with . With evaporation turned off, of course, this mass fraction remains constant over time. In each case, we present results for both a relatively long initial cooling timescale () for the planetary atmosphere, corresponding to a “cold start” planet formation scenario, and a relatively short , corresponding to a “hot start” (see §5).

We see that, for all three core masses, the final surface temperature after a Gyr is K, and the final pressure is nearly 10 bar (with variations in making negligible difference to these results). These are well above the critical point of water ( K, bar); as such, any surface water would exist as a supercritical fluid, not as a liquid (see phase diagram of water, Fig. 20; constructed from phase transition equations provided by Martin Chaplin (pvt. comm. 2015)^{6}^{6}6More detailed phase diagram for water by M. Chaplin available online at:

http://www1.lsbu.ac.uk/water/water_phase_diagram.html#bd. The equations used for the phase boundaries can be seen by hovering over the boundaries in the latter plot.). Hence, as stated earlier, removal of a large fraction of the H/He envelope via hydrodynamic escape, which will reduce and , is required for habitability when is of order a percent. The degree to which this can happen is the focus of our analysis in §§6.1 and 6.2 below, since, as discussed, such relatively high initial mass fractions are commensurate with current data.

Next, in Figs. 9, 10 and 11, we plot the evolution of surface conditions for a much lower initial atmospheric mass: 4. Again, in the absence of evaporation, this fraction remains constant. Now the surface conditions after a Gyr are much more benign, with K and bar. As the phase diagram in Fig. 20 shows, surface water should exist as a liquid under these conditions. However, realistic levels of evaporation will in fact remove most or all of this tenuous atmosphere: as we will show below, strong hydrodynamic escape very quickly reduces the atmospheric mass-fraction to 10; even Jeans escape can then remove the remainder on Gyr timescales.

Any potential long-term habitability will then depend on the quantity, properties and evolution of (outgassed) secondary atmospheres, similar to the case on Earth and Mars (e.g. Elkins-Tanton, 2008). As such, for these low initial atmospheric mass fractions, we will plot the evolution of the mass fraction down to 10, to make the point that the primordial H/He atmosphere is likely to be entirely lost, but will abstain from any further discussion of the final surface conditions, which will depend on secondary atmospheres whose modeling is beyond the scope of this paper (although several authors have begun to investigate this stage theoretically; e.g., Elkins-Tanton & Seager, 2008; Elkins-Tanton, 2011). Nonetheless, since secondary atmospheres are at least present on all the terrestrial planets and massive moons in the solar system, we will consider planets whose H/He envelopes have been completely stripped to still be “potentially habitable”.

The results for K ( AU) without evaporation (not plotted) are similar: while the surface temperatures after a Gyr are now somewhat lower, and are still far too high for liquid surface water when ; for , the surface conditions are again conducive to liquid water, but with the same caveats as above.

### 6.2 Evaporation at Inner Edge of HZ ( K)

We now turn the evaporation on, and once again examine the results for core masses of 0.8, 0.9 and 1.0 M at K.

0.8 M: Fig. 12 shows the evolution of the atmospheric mass-fraction and planetary radius for = 0.8 M, for (top panels) and 4 (bottom panels). In both cases, evaporation strips the atmosphere down to a mass fraction of 10 extremely rapidly, in 10 Myr. This steep decline is due to runaway evaporation, wherein the cooling timescale exceeds the evaporation timescale: the atmosphere cannot contract quickly enough to avoid evaporation by falling deeper into the planet’s gravitational potential well (c.f. Baraffe et al., 2005). Note that the runaway occurs much earlier than the median age of 10 Myr estimated for AD Leo (and for other early to mid-M field dwarfs evincing similar saturated X-ray activity; see §2.2); consequently, this evaporation result is unaffected by our simplistic assumption that activity remains saturated beyond a few 100 Myr up to a Gyr (instead of declining at later times).

For a 0.8 M core, therefore, we expect any primordial H/He atmosphere with initial mass fraction 1% to be completely lost: eroded down to 10 in 10 Myr by runaway hydrodynamic evaporation, as shown, and the minuscule remainder efficiently peeled off by continuing hydrodynamic and/or Jeans escape over a less than Gyr timescale.

0.9 M: Fig. 13 shows the evolution in atmospheric mass fraction and planetary radius for = 0.9 M. For , we see that the “hot” and “cold” start cases behave very differently. The atmosphere in the former is initially relatively bloated, resulting in a large fraction of it being blown off almost instantaneously when evaporation is switched on (a phenomenon termed “boil-off” by Owen & Wu 2015; see §4.3). For reasons of clarity, this initial boil-off, which occurs over the first few time-steps of our calculations (in 1.5 Myr), is masked out in Fig. 13 and subsequent figures. Essentially one cannot build stable models with envelope mass fractions and cooling times as short as desired – they are not globally thermodynamically stable in the presence of mass-loss. Specifically, in the pressure confining environment of the parent disc the planet can acquire a % envelope, but once this pressure confinement is removed such a massive envelope cannot be gravitationally bound to the planet (see Owen & Wu, 2015). Therefore, the bolometrically driven mass-loss allows the planet to readjust to a stable state that is in quasi thermodynamic equilibrium. The plots trace the evolution of the planet after boil-off ends. The remainder then undergoes runaway hydrodynamic escape, reducing the mass fraction to 10 in 10 Myr. For a “cold start”, on the other hand, the initial atmosphere is much more compact, and the pace of hydrodynamic evaporation is thus far more leisurely; by 30 Myr, evaporation is too slight to make any significant difference, and the atmospheric mass fraction settles down at 10. Jeans escape, which only becomes efficient at mass fractions 10, has no discernible effect on the surviving atmosphere. Note as well that the 30 Myr it takes for the mass fraction to reach approximate steady-state is at the lower limit of the 10 Myr age estimate for AD Leo; our result is thus minimally affected by our naive assumption that activity continues to be saturated beyond a few 100 Myr instead of declining at such late times. Finally, cores with 10 are evaporated down to 10 in only a few Myr, regardless of initial cooling timescale, just as in the 0.8 M case.

The evolution of surface conditions for the case is shown in Fig. 14. For a “cold start” 0.9 M core, we see that and after a Gyr are 530 K and 630 bar respectively (note that the slight continuing evolution of these quantities beyond 100 Myr is predominantly due to ongoing thermal evolution of the planet, and not any appreciable evaporative mass loss). This [, ] combination is conducive to liquid surface water (Fig. 20), and within the range found in deep-ocean hydrothermals vents on Earth.

Fig. 14 also shows the “hot start” core with the same for comparison. The [,] [350 K, 10 bar] achieved at the end of our calculations is of no particular significance, since we expect the tiny remaining atmospheric mass fraction of 10 to be completely removed subsequently by even Jeans escape; we only note that the surface temperature is beginning to flatten out since it cannot fall below the equilibrium temperature (with albedo assumed to be zero) of 300 K at this orbital radius, while the pressure continues to plummet, which will eventually lead to any liquid surface water boiling into vapour. Of more physical interest is the fact that in the “hot start” case is lower than in the “cold start” one, once the rapid initial boil-off stage (which is not plotted in these figures) has ended. This is because hydrodynamic escape requires mass continuity at the sonic point; when a substantial fraction of the atmosphere is blown off from the top, therefore, the underlying layers swiftly expand to take its place, and the upward advection of heat in this process causes the surface to cool (see discussion in Owen & Wu, 2016).

In summary, at K, we expect liquid water to survive on Gyr timescales on the surface of a “cold start” 0.9 M core with ; the partial evaporation of the H/He atmosphere engenders habitable conditions in this case. A “hot start” 0.9 M core with the same , or a core with much smaller and either a “hot” or a “cold” start, however, is unlikely to retain any of its primordial H/He atmosphere; habitability conditions in this case will depend on any secondary atmosphere that arises. Therefore, a 0.9 M core appears to be at the transition between a planet becoming habitable due to evaporation, and remaining uninhabitable due to retention of a significant H/He envelope.

1.0 M: Finally, Fig. 15 shows the the evolution in atmospheric mass fraction and planetary radius for = 1.0 M. There is relatively little evaporation when , with both “hot” and “cold start” cores equilibrating at a mass fraction of 210 (there is some initial boil-off in the “hot start” case, but hardly enough to initiate runaway evaporation; the planet’s gravity is too strong). When 10, on the other hand, runaway evaporation reduces the fraction to 10 in just a few Myr, as in the 0.8 and 0.9 M cases.

In Fig. 16 we plot the evolution of surface conditions for . Both “hot” and “cold” start cores end up with very similar surface temperatures and pressures after a Gyr, with [, ] [800 K, 210 bar]. Fig. 20 demonstrates that these conditions are too extreme for liquid water; any surface water here can only exist as a supercritical fluid.

### 6.3 Evaporation at Outer Edge of HZ ( K)

We now investigate the effects of hydrodynamic evaporation at K. The single case of a 0.8 M core will serve to illustrate the main results here.

0.8 M: We plot the evolution in atmospheric mass fraction and planetary radius for this core mass in Fig. 17. There is hardly any evaporation for either or 4: in both cases, the final mass fraction after a Gyr is very similar to the initial value. Figs. 18 and 19 show the evolution in surface conditions for these two initial mass fractions. For of 10, the final [, ] after a Gyr is [800 K, 410 bar]; as the phase diagram in Fig. 20 indicates, surface water will exist as a supercritical fluid instead of a liquid under these conditions. For of 410, on the other hand, the equilibrium [, ] after a Gyr is [400 K, 280 bar], which is conducive to liquid water.

We can easily extrapolate from these results, and the preceding ones, to deduce the implications for 0.9 and 1.0 M cores at this radial separation. First, surface temperature and pressure always increase with core mass, for a given and orbital radius. Thus, since a 0.8 M core with an initial H/He atmospheric mass fraction of a percent cannot harbour liquid water at K, higher mass cores with the same initial mass fraction cannot either. Second, in the absence of evaporation, and are very similar for 0.8–1.0 M cores at any fixed and orbital radius (e.g., see the no-evaporation cases plotted in Fig. 20 for K: the for 0.8–1.0 M cores are nearly identical, and their vary by 100 K, for both of 10 and 410). We found above that there is hardly any evaporation in the 0.8 M case even with 10 at this orbital radius; therefore, since the evaporative rate decreases with increasing core mass (as the gravitational potential well deepens), we can invoke the no-evaporation results here for the entire range 0.8–1.0 M. This in turn implies that 0.9 and 1.0 M cores with 10 will have [, ] very close to that of the 0.8 M core with the same initial mass fraction, and can thus also harbour liquid surface water on Gyr timescales. These inferences are indeed what our more detailed calculations show.

## 7 Discussion

The central result of the preceding analysis is that solid cores made of 2/3 rock + 1/3 iron, with mass 1 M, cannot lose enough of their envelopes to become habitable in the classical HZ of M dwarfs, if they are born with significant H/He envelopes. The fundamental reason is that hydrodynamic escape is quenched in such planets while a sizeable portion of the H/He envelope still lingers, and subsequent Jeans escape – which can at best extract a fractional mass of 10 over a Gyr – is too feeble to remove this remainder. Only cores 1 M (with the precise limiting mass depending on where in the HZ one is) can undergo hydrodynamic escape all the way down to their surface in spite of having initial H/He mass-fractions of 1%, and may thus be habitable (either due to a tenuous remnant H/He envelope, or a similarly wispy secondary atmosphere) in the classical HZ of M dwarfs.

### 7.1 Implications for Habitable Planets around M Dwarfs

To summarise: Kepler data imply that small (terrestrial-size) planets are ubiquitous at small orbital separations, with nearly every star (statistically speaking) hosting one such “Kepler” planet; around M dwarfs, a sizeable fraction of these close-in planets reside within the classical HZ (see also further below). Concurrently, both observations and theory indicate that such planets are usually born with H/He envelopes with a mass-fraction 1%, and our calculations imply that planets with cores comprising 2/3 rock and 1/3 iron, with core masses 1 M and initial H/He envelope mass-fractions 1%, cannot lose enough of their envelopes to be habitable within the classical HZ of M dwarfs. With these considerations, we propose three possible classes of solid-core habitable planets around M dwarfs:

(1) Sub-Earth mass rock/iron-core planets within the classical HZ: We found above that sub-Earth mass rock/iron cores – with mass 0.9 M near the inner edge of the classical HZ of M dwarfs, and 0.8 M near the outer edge – can be stripped of 1% initial H/He envelopes within a Gyr. Such Venus-mass planets may therefore be habitable within the classical HZ of M dwarfs, if they acquire secondary atmospheres like the solar system terrestrial planets (the HZ results of Kopparapu et al. (2013) can then be applied to such planets).

(2) Sub-Earth to super-Earth mass ice-core planets within the classical HZ: Our discussion so far has focussed on rock/iron cores. Naively, the results from (1) above should apply to ice cores as well (but see below), except for a higher threshold mass (the limiting planet mass below which 1% initial H/He envelopes will be stripped). This is easily seen by examining Fig. 4, where the threshold mass is defined as the intersection of the Kn=1 line and the locus of a core with fixed density. For a given core mass, an ice core, with a density 10 times smaller than a rock/iron one, will have the latter locus shifted to the right (i.e., to a larger radius) by a factor of 102, resulting in a threshold mass larger by a factor of 2 as well from Fig. 4. Thus, to zeroeth order, one expects ice-core planets, with a somewhat (factor of 2) larger limiting mass than rock/iron core ones, to also be habitable within the classical HZ of M dwarfs. However, a correct analysis of evaporation here must also account for the opacity of steam (liberated from the core’s surface) in addition to H/He, which may strongly affect the outcome (since oxygen – in this case within HO molecules – is a strong X-ray absorber). This will be the subject of future work. Note also that icy cores can only form beyond the ice-line, situated at 1.5 AU for a 0.4 M M dwarf at an age of order a Myr^{7}^{7}7Using the formula supplied by Ida & Lin (2008) for the ice-line orbital separation: AU; with a stellar bolometric luminosity of for a 0.4 M star at an age of 1.5 Myr, from the evolutionary tracks of Baraffe et
al. (1998). Thus ice-cores must migrate inwards for the above discussion to be pertinent.

(3) Planets formed after gas-disk dispersal: Finally, our entire analysis in this paper is for planets that form while a significant amount of gas still remains in the surrounding primordial disk, allowing them to accrete relatively massive H/He envelopes; Kepler data, as discussed, imply that such planets are ubiquitous. However, the terrestrial planets in our own solar system evince signatures of having coalesced after the gas disk dispersed (as the radiometricly determined age of the earth – Manhes et al. 1980 – is estimated to be Myr younger than that of the solar-system, e.g. Bouvier & Wadhwa 2010, well after the gas disc would have dispersed). If solid-core planets around M dwarfs can form in a similar fashion, then they might well be habitable in the classical HZ of M dwarfs (and slightly interior / exterior to it depending on planetary mass); the work by Kopparapu et al. (2013), examining solely secondary atmospheres, is more directly applicable in that case than this paper.

We note that there is, prima facie, another interesting possibility. While we find that rock/iron cores 1 M with 1% initial H/He envelopes retain too large a fraction of these envelopes to be habitable within the classical HZ of M dwarfs, they will be stripped of such atmospheres if their orbits are much smaller. Simultaneously, Kopparapu et al. (2013) show that the inner boundary of the HZ (set by either the moist greenhouse or runaway greehouse effect) for such super-Earth mass planets will be closer to the star than for Earth-mass ones (where the latter are used, by definition, to calculate the classical HZ limits). Therefore, one might expect an overlap in space, interior to the classical HZ, between where a super-Earth of given mass can be stripped of a 1% H/He envelope, and where it will be habitable with a secondary atmosphere before greenhouse effects become overwhelming; super-Earths would then be habitable within this overlap region. Unfortunately, such an overlap is implausible. For solar-type stars, Kopparapu et al. (2013) show that the inner boundary of the HZ moves inwards from 0.99 AU for a 1 M planet to 0.94 AU for a 10 M one; i.e., only a very small shift of 0.05 AU for a factor of 10 increase in planetary mass. The effect should be similar around M dwarfs. This tiny decrease in orbital separation will have negligible impact on the incident X-ray flux and thus on the mass-loss rates; as such, we do not expect any HZ for super-Earths interior to the classical HZ around M dwarfs. Conversely, for sub-Earth mass planets, the inner boundary of the HZ moves slightly outwards, as Kopparapu et al. (2013) show; thus, we do not expect any HZ for these planets interior to the classical HZ either.

There are two main ingredients that determine the actual frequency of habitable planets around M dwarfs. The first, as we have shown, is a rigorous treatment of evaporation, using hydrodynamic models that account for radiative cooling as well as the transition to Jeans escape, and including the thermal evolution of the planet. A simplistic “energy-limited” formalism, that moreover assumes that the escape is always hydrodynamic (e.g., Luger et al. 2015a for M dwarfs), leads to a gross overestimation of evaporative mass-loss rates, and thereby an overly optimistic appraisal of the fraction of planets that can be rendered habitable by stripping their primordial H/He envelopes. Our results are in direct conflict with those of Luger et al. (2015a). The latter authors’ calculations suggest that in the majority of cases, M dwarfs can completely strip envelopes with mass fraction 10% (and even larger ones in some cases) at the inner edge of the HZ, and mass fractions 1% at the outer edge of the HZ, from 1–2 M cores. In other words, they conclude that evaporation can generate a plethora of potentially habitable Earth-mass planets around M-dwarfs. Our results indicate that this is not the case: it is the simplifications made by Luger et al. (2015a) thet lead them to conclude otherwise. The main reason is their choice of a constant efficiency energy-limited mass-loss prescription. As we have discussed throughout this work, at late times when the XUV flux declines and the envelope mass fraction drops below 1% (i.e., when the planet’s radius shrinks rapidly with decreasing envelope mass), the efficiency of hydrodynamic mass-loss decreases significantly (c.f. Owen & Wu, 2013); furthermore, the flow transitions to non-hydrodynamic (Jeans escape) at late times. Neglecting these effects, Luger et al. (2015a) find that their planets can undergo runaway mass-loss at ages Gyr; including these effects, however, we find that a 1 M planet can only ever undergo runaway mass-loss at early times ( Gyr), when the flux is high enough to permit hydrodynamic evaporation. In fact, if we adopt an energy-limited mass-loss rate with no hydrodynamic cut-off at low fluxes, we indeed recover the runaway mass-loss at late times found by Luger et al. (2015a). This is demonstrated in Fig. 21, where we compare the results of our calculations (solid line) to those for the energy-limited case with efficiencies () of 5% (dot-dashed), 10% (dashed) & 25% (dotted); note that Luger et al. (2015a)’s default choice is = 30%. The plot shows that, while our planet loses negligible mass at late times due to the hydrodynamic cut-off, the energy-limited cases can undergo unphysical runaway mass-loss at these ages. We note that while an efficiency of produces a good match to the final envelope mass-fraction at 10 Gyr in this particular example, this is not a generally applicable result: indeed, the plot shows that the 10% case (dashed line) is also starting to enter a run-away phase at the end of the calculation. In general, there is no single efficiency value that can mimic the trends in evaporation in our results, especially the hydrodynamic cut-off, which turns out to be extremely important for evaporation and the question of potential habitability. Using an energy-limited formalism to determine evaporation rates when the flow is close to transitioning from hydrodynamic to non-hydrodynamic can lead to rather poor predictions. Finally, we note that our prediction that planets 1 M will retain most of their initial voluminous H/He envelopes, in the classical HZ of M dwarfs, should be directly testable with TESS, which will target a large number of M dwarfs.

The second ingredient that controls the actual prevalence of habitable planets around M dwarfs is the frequency of their planets as a function of planetary mass, at the small end of the planetary spectrum. A recent careful analysis by Morton & Swift (2014) indicates that the occurrence of M dwarf planets continues to increase to planetary radii below 1 R, without any strong evidence for a downturn in the planet radius function that was suggested by previous studies. If this is true, then it bodes well for the frequency of potentially habitable planets orbiting M dwarfs, since the smallest (lowest-mass) planets will be most easily stripped of their initial H/He envelopes. Indeed, applying their improved calculations to previous work by Dressing & Charbonneau (2013) and Kopparapu (2013), Morton & Swift (2014) already show that the frequency of Earth-size planets in the HZ of M dwarfs – calculated by Dressing & Charbonneau (2013) to be 0.15 planets per star – increases to 0.25–0.8, a remarkable fraction. If this frequency continues to increase towards smaller planets, as Morton & Swift (2014) suggest, then evaporation can indeed create a bonanza of potentially habitable planets around these cool red dwarfs. However, the latter results are based on a quite limited sample; they need to be verified and sharpened with more data from the K2 mission, as well as upcoming ones such as TESS.

## 8 Summary

Our key results are as follow:

1) Jeans escape can only remove a very small H/He mass-fraction from low-mass exoplanets in the HZ of M dwarfs (10 over a Gyr): hydrodynamic evaporation is essential for removing significant amounts of H/He.

2) Previously derived mass-loss rates, based on the “energy-limited” (or “locally energy-limited”) formalism, and assuming moreover that escape is always in the hydrodynamic limit, are significant overestimations. Our improved model, accounting for both radiative losses and the transition from hydrodynamic to Jeans escape, and including the thermal evolution of the planet, yields much lower mass loss rates.

3) In particular, for cores made of 2/3 rock + 1/3 iron, we find that only sub-Earth mass cores – 0.9 M near the inner edge of the classical HZ of M dwarfs, and 0.8 M near the outer edge – can be stripped of 1% initial H/He envelopes. Rock/iron cores with mass 1 M, and born with H/He envelope mass-fractions 1%, cannot lose sufficient envelope mass to become habitable within the classical HZ of M dwarfs. Our prediction that 1 M cores in the HZ of M dwarfs should still be shrouded by their voluminous natal H/He envelopes over Gyr timescales will be directly testable with TESS.

4) We propose three classes of potentially habitable planets in the HZ of M dwarfs: (i) planets with sub-Earth mass rock/iron cores: since they can be stripped of 1% H/He envelopes within a Gyr, they may be habitable if they can acquire suitable secondary atmospheres, like solar system terrestrial planets; (ii) planets with ice-cores, with an upper limit to the core mass somewhat higher than the sub-Earth mass for rock/iron cores, may also be able to lose their primordial H/He envelopes and thus become habitable; however, steam opacity must be included in the evaporation calculations to verify this possibility; and (iii) planets (possibly like the terrestrial ones in the solar system) that form after gas disk dispersal, and thus have only tenuous secondary atmospheres, with very little primordial H/He.

## References

- Baraffe et al. (1998) Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403
- Baraffe et al. (2005) Baraffe, I., Chabrier, G., Barman, T. S. et al. 2005, A&A, 436, L47
- Bodenheimer & Lissauer (2014) Bodenheimer, P., & Lissauer, J. J. 2014, ApJ, 791, 103
- Bouvier & Wadhwa (2010) Bouvier, A., & Wadhwa, M. 2010, Nature Geoscience 3, 637
- Chadney et al. (2015) Chadney, J. M., Galand, M., Unruh, Y. C. et al. 2015, Icarus, 250, 357
- Delfosse et al. (1998) Delfosse, X., Forveille, T., Perrier, C., & Mayor, M. 1998, A&A, 331, 581
- Dressing & Charbonneau (2013) Dressing, C. D. & Charbonneau, D. 2013, ApJ, 767, 95
- Dressing et al. (2015) Dressing, C. D., Charbonneau, D., Dumusque, X. et al. 2015, ApJ, 800, 135
- Elkins-Tanton (2008) Elkins-Tanton L. T., 2008, E&PSL, 271, 181
- Elkins-Tanton & Seager (2008) Elkins-Tanton L. T., Seager S., 2008, ApJ, 685, 1237-1246
- Elkins-Tanton (2011) Elkins-Tanton L. T., 2011, Ap&SS, 332, 359
- Ercolano et al. (2003) Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X.-W. 2003, MNRAS, 340, 1136
- Ercolano et al. (2005) Ercolano, B., Barlow, M. J., & Storey, P. J. 2005, MNRAS, 362, 1038
- Ercolano et al. (2008a) Ercolano, B., Young, P. R., Drake, J. J., & Raymond, J. C. 2008, ApJS, 175, 534
- Ercolano et al. (2008b) Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008, ApJ, 688, 398
- Ercolano et al. (2011) Ercolano B., Bastian N., Spezzi L., Owen J., 2011, MNRAS, 416, 439
- Erkaev et al. (2007) Erkaev, N. V., Kulikov, Y. N., Lammer, H. et al. 2007, A&A, 472, 329
- Erkaev et al. (2013) Erkaev, N. V., Lammer, H., Odert, P. et al. 2013, Astrobiology, 13, 1011
- Erkaev et al. (2016) Erkaev, N. V., Lammer, H., Odert, P. et al. 2016, arXiv, arXiv:1601.00452
- Fortney et al. (2007a) Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661
- Fortney et al. (2007b) Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 668, 1267
- Fressin et al. (2013) Fressin, F., Torres, G., Charbonneau, D. et al. 2013, ApJ, 766, 81
- Güdel (2004) Güdel, M. 2004, A&A Rev., 12, 71
- Guillot (2010) Guillot, T. 2010, A&A, 520, AA27
- Hadden & Lithwick (2014) Hadden, S., & Lithwick, Y. 2014, ApJ, 787, 80
- Haworth et al. (2016) Haworth T. J., Clarke C. J., Owen J. E., 2015, arXiv, arXiv:1512.02234
- Hayes et al. (2006) Hayes, J. C., Norman, M. L., Fiedler, R. A. et al. 2006, ApJS, 165, 188
- Heller, Leconte, & Barnes (2011) Heller R., Leconte J., Barnes R., 2011, A&A, 528, A27
- Hernández et al. (2007) Hernández J., et al., 2007, ApJ, 662, 1067
- Howard et al. (2012) Howard, A. W., Marcy, G. W., Bryson, S. T. et al. 2012, ApJS, 201, 15
- Hunten (1973) Hunten, D. M. 1973, Journal of Atmospheric Sciences, 30, 1481
- Hu, Seager, & Yung (2015) Hu R., Seager S., Yung Y. L., 2015, ApJ, 807, 8
- Ida & Lin (2008) Ida, S. & Lin, D. N. C. 2008, ApJ, 685, 595
- Ikoma & Hori (2012) Ikoma M., Hori Y., 2012, ApJ, 753, 66
- Jackson et al. (2012) Jackson, A. P., Davis, T. A. & Wheatley, P. J. 2012, MNRAS, 422, 2024
- Jin et al. (2014) Jin, S., Mordasini, C., Parmentier, V. et al. 2014, ApJ, 795, 65
- Johnstone et al. (2015) Johnstone, C. P., Güdel, M., Stökl, A. et al. 2015, arXiv:1511.03647
- Johnson, Volkov, & Erwin (2013) Johnson R. E., Volkov A. N., Erwin J. T., 2013, ApJ, 768, L4
- Jontof-Hutter et al. (2014) Jontof-Hutter, D., Lissauer, J. J., Rowe, J. F., & Fabrycky, D. C. 2014, ApJ, 785, 15
- Jontof-Hutter et al. (2015) Jontof-Hutter D., et al., 2015, arXiv, arXiv:1512.02003
- Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., Reynolds, R. T. 1993, Icarus, 101, 108
- Kashyap & Drake (2000) Kashyap V., Drake J. J., 2000, BASI, 28, 475
- Kopparapu (2013) Kopparapu, R. K. 2013, ApJ, 767, L8
- Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F. et al. 2013, ApJ, 765, 131
- Lammer et al. (2003) Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121
- Lammer et al. (2009) Lammer, H., Odert, P., Leitzinger, M., et al. 2009, A&A, 506, 399
- Lammer et al. (2013) Lammer, H., Erkaev, N. V., Odert, P., et al. 2013, MNRAS, 430, 1247
- Lammer et al. (2014) Lammer, H., Stökl, A., Erkaev, N. V., et al. 2014, MNRAS, 439, 3225
- Leconte et al. (2015) Leconte, J., Wu, H., Menou, K., & Murray, N. 2015, Science, 347, 632
- Lee et al. (2014) Lee, E. J., Chiang, E., & Ormel, C. W. 2014, ApJ, 797, 95
- Lee & Chiang (2015) Lee, E. J., & Chiang, E. 2015, arXiv:1508.05096
- Liu et al. (2015) Liu S.-F., Hori Y., Lin D. N. C., Asphaug E., 2015, ApJ, 812, 164
- Lopez et al. (2012) Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761, 59
- Lopez & Fortney (2013) Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2
- Lopez & Fortney (2014) Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1
- Luger et al. (2015a) Luger, R., Barnes, R., Lopez, E., et al. 2015a, Astrobiology, 15, 57
- Luger et al. (2015b) Luger, R., Barnes, R., Lopez, E., et al. 2015b, Astrobiology, 15, 57
- Manhes et al. (1980) Manhes G., Allègre C. J., Dupré B., Hamelin B., 1980, E&PSL, 47, 370
- Marcy et al. (2014) Marcy, G. W., Isaacson, H., Howard, A. W., et al. 2014, ApJS, 210, 20
- Marley et al. (2007) Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541
- Morton & Swift (2014) Morton, T. D., & Swift, J. 2014, ApJ, 791, 10
- Mullally et al. (2015) Mullally, F., Coughlin, J. L., Thompson, S. E., et al. 2015, ApJS, 217, 31
- Murray-Clay et al. (2009) Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23
- Owen et al. (2010) Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415
- Owen, Ercolano, & Clarke (2011) Owen J. E., Ercolano B., Clarke C. J., 2011, MNRAS, 412, 13
- Owen & Jackson (2012) Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931
- Owen & Wu (2013) Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105
- Owen & Adams (2014) Owen, J. E., & Adams, F. C. 2014, MNRAS, 444, 3761
- Owen & Alvarez (2016) Owen J. E., Alvarez M. A., 2016, ApJ, 816, 34
- Owen & Wu (2016) Owen J. E., Wu Y., 2016, ApJ, 817, 107
- Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
- Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
- Petigura et al. (2013) Petigura, E. A., Marcy, G. W., & Howard, A. W. 2013, ApJ, 770, 69
- Pierrehumbert & Gaidos (2011) Pierrehumbert, R., & Gaidos, E. 2011, ApJ, 734, L13
- Reiners & Mohanty (2012) Reiners, A. & Mohanty, S. 2012, ApJ, 746, 43
- Rogers et al. (2011) Rogers L. A., Bodenheimer P., Lissauer J. J., Seager S., 2011, ApJ, 738, 59
- Rogers (2015) Rogers, L. A. 2015, ApJ, 801, 41
- Sanz-Forcada et al. (2011) Sanz-Forcada, J., Micela, G., & Ribas, I. et al. 2011, A&A, 532, A6
- Scalo et al. (2007) Scalo J., et al., 2007, AsBio, 7, 85
- Shematovich et al. (2014) Shematovich V. I., Ionov D. E., Lammer H., 2014, A&A, 571, A94
- Shkolnik et al. (2009) Shkolnik, E., Liu, M. C., & Reid, I. N. 2009, ApJ, 699, 649
- Silburt et al. (2015) Silburt, A., Gaidos, E., & Wu, Y. 2015, ApJ, 799, 180
- Spiegel & Burrows (2012) Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174
- Stökl, Dorfi, & Lammer (2015) Stökl A., Dorfi E., Lammer H., 2015, A&A, 576, A87
- Stone & Norman (1992) Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753
- Stone & Proga (2009) Stone, J. M., & Proga, D. 2009, ApJ, 694, 205
- Tian et al. (2005b) Tian F., Toon O. B., Pavlov A. A., De Sterck H., 2005, Sci, 308, 1014
- Tian et al. (2005a) Tian, F., Toon, O. B., Pavlov, A. A., & De Sterck, H. 2005, ApJ, 621, 1049
- Tian et al. (2008) Tian F., Kasting J. F., Liu H.-L., Roble R. G., 2008, JGRE, 113, E05008
- Trammell, Li, & Arras (2014) Trammell G. B., Li Z.-Y., Arras P., 2014, ApJ, 788, 161
- Tripathi et al. (2015) Tripathi A., Kratter K. M., Murray-Clay R. A., Krumholz M. R., 2015, ApJ, 808, 173
- Tu et al. (2015) Tu, L., Johnstone, C. P., Güdel, M., & Lammer, H. 2015, A&A, 577, L3
- Velli (1994) Velli, M. 1994, ApJ, 432, L55
- Watson et al. (1981) Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150
- Weiss & Marcy (2014) Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, L6
- Wolfgang & Lopez (2015) Wolfgang, A., & Lopez, E. 2015, ApJ, 806, 183
- Wolfgang et al. (2015) Wolfgang, A., Rogers, L. A., & Ford, E. B. 2015, arXiv:1504.07557
- Wu & Lithwick (2013) Wu, Y., & Lithwick, Y. 2013, ApJ, 772, 74
- Xu & McCray (1991) Xu, Y., & McCray, R. 1991, ApJ, 375, 190