Atmospheric escape from the TRAPPIST-1 planets and implications for habitability

TL;DR

This study uses numerical simulations and analytical estimates to assess whether the seven TRAPPIST-1 planets can retain their atmospheres over billions of years, impacting their potential habitability.

Contribution

It provides the first combined numerical and analytical analysis of atmospheric escape for all TRAPPIST-1 planets, highlighting their potential to sustain atmospheres.

Findings

Outer planets can retain atmospheres over billion-year timescales

Numerical and analytical escape rates are in good agreement

Implications for habitability, abiogenesis, and future exoplanet searches

Abstract

The presence of an atmosphere over sufficiently long timescales is widely perceived as one of the most prominent criteria associated with planetary surface habitability. We address the crucial question as to whether the seven Earth-sized planets transiting the recently discovered ultracool dwarf star TRAPPIST-1 are capable of retaining their atmospheres. To this effect, we carry out numerical simulations to characterize the stellar wind of TRAPPIST-1 and the atmospheric ion escape rates for all the seven planets. We also estimate the escape rates analytically and demonstrate that they are in good agreement with the numerical results. We conclude that the outer planets of the TRAPPIST-1 system are capable of retaining their atmospheres over billion-year timescales. The consequences arising from our results are also explored in the context of abiogenesis, biodiversity, and searches for future exoplanets. In light of the many unknowns and assumptions involved, we recommend that these conclusions must be interpreted with due caution.

Figures

Tables

Table 1. Table 1: Ion escape rates in sec -1 .

	O+	O${}_{}{}^{+}{}_2{}^{}$	CO${}_{}{}^{+}{}_2{}^{}$	Total
Maximum total pressure
Trappist-1b	5.56$\times$1027	2.09$\times$1026	1.52$\times$1026	5.92$\times$1027
Trappist-1c	1.54$\times$1027	1.38$\times$1026	1.32$\times$1026	1.81$\times$1027
Trappist-1d	1.29$\times$1027	3.80$\times$1025	1.14$\times$1025	1.34$\times$1027
Trappist-1e	7.01$\times$1026	2.83$\times$1025	1.10$\times$1025	7.40$\times$1026
Trappist-1f	5.23$\times$1026	3.37$\times$1025	1.19$\times$1025	5.68$\times$1026
Trappist-1g	2.17$\times$1026	2.71$\times$1025	1.32$\times$1025	2.58$\times$1026
Trappist-1h	1.06$\times$1026	1.65$\times$1025	6.98$\times$1024	1.29$\times$1026
Minimum total pressure
Trappist-1b	9.33$\times$1026	4.99$\times$1025	2.92$\times$1025	1.01$\times$1027
Trappist-1c	4.23$\times$1026	9.22$\times$1025	2.76$\times$1025	5.42$\times$1026
Trappist-1d	2.81$\times$1026	3.07$\times$1025	1.04$\times$1025	3.23$\times$1026
Trappist-1e	2.20$\times$1026	4.19$\times$1025	1.25$\times$1025	2.74$\times$1026
Trappist-1f	1.88$\times$1026	4.30$\times$1025	1.10$\times$1025	2.42$\times$1026
Trappist-1g	9.33$\times$1025	5.85$\times$1025	1.38$\times$1025	1.66$\times$1026
Trappist-1h	4.52$\times$1025	2.69$\times$1025	4.39$\times$1024	7.66$\times$1025

Table 2. Table 2: Chemical reactions and associated rates adapted from Dong et al. ( 2017b ) . The photoionization frequencies are rescaled to account for the EUV fluxes received at each of the TRAPPIST-1 planets based on the estimates provided in Bolmont et al. ( 2017 ) and Bourrier et al. ( 2017 ) .

Chemical Reaction		Rate Coefficient
Primary Photolysis in s-1
	CO2 + $h\nu$ $\to$ CO${}_{}{}^{+}{}_2{}^{}$ + $e^{-}$	see table caption
	O + $h\nu$ $\to$ O+ + $e^{-}$	see table caption
Ion-neutral Chemistry in cm3 s-1
	CO${}_{}{}^{+}{}_2{}^{}$ + O $\to$ O${}_{}{}^{+}{}_2{}^{}$ + CO	$1.64\times {10}^{-10}$
	CO${}_{}{}^{+}{}_2{}^{}$ + O $\to$ O+ + CO2	$9.60\times {10}^{-11}$
	O+ + CO2 $\to$ O${}_{}{}^{+}{}_2{}^{}$ + CO	$1.1\times {10}^{-9}$ (800/Ti)0.39
	H+ + O $\to$ O+ + H	$5.08\times {10}^{-10}$
Electron Recombination Chemistry in cm3s-1
	O${}_{}{}^{+}{}_2{}^{}$ + $e^{-}$ $\to$ O + O	$7.38\times {10}^{-8}$ (1200/Te)0.56
	CO${}_{}{}^{+}{}_2{}^{}$ + $e^{-}$ $\to$ CO + O	$3.10\times {10}^{-7}$ (300/Te)0.5

Table 3. Table 3: Stellar wind parameters for the TRAPPIST-1 planets

	Nsw(cm-3)	Tsw(K)	Vsw(km/s)	IMF(nT)
Maximum total pressure
Trappist-1b	6.59$\times$104	2.01$\times$106	$(-470,80,-1)$	$(381,81,-147)$
Trappist-1c	2.99$\times$104	1.68$\times$106	$(-527,68,0)$	$(210,42,-111)$
Trappist-1d	1.20$\times$104	1.59$\times$106	$(-566,56,6)$	$(-129,15,-50)$
Trappist-1e	5.79$\times$103	1.26$\times$106	$(-604,50,3)$	$(-149,13,-42)$
Trappist-1f	2.99$\times$103	1.02$\times$106	$(-624,44,3)$	$(-98,7,-34)$
Trappist-1g	1.95$\times$103	8.92$\times$105	$(-637,40,2)$	$(-69,6,-28)$
Trappist-1h	9.52$\times$102	7.17$\times$105	$(-657,37,2)$	$(-44,2,-25)$
Minimum total pressure
Trappist-1b	4.28$\times$103	2.40$\times$106	$(-803,80,10)$	$(-2206,96,-64)$
Trappist-1c	2.09$\times$103	2.35$\times$106	$(-871,68,-3)$	$(-1192,63,-64)$
Trappist-1d	1.06$\times$103	2.00$\times$106	$(-923,56,-5)$	$(-641,41,-59)$
Trappist-1e	5.68$\times$102	1.72$\times$106	$(-972,50,3)$	$(-379,32,-38)$
Trappist-1f	2.93$\times$102	1.44$\times$106	$(-1017,44,19)$	$(-218,23,-25)$
Trappist-1g	2.05$\times$102	1.30$\times$106	$(-1032,40,2)$	$(-159,20,-19)$
Trappist-1h	1.09$\times$102	1.07$\times$106	$(-1049,37,-16)$	$(-89,14,-11)$

Equations

$$H_x=\frac{kT}{m{g}_x},$$

$$N_x\propto {\left(\frac{R_x}{a}\right)}^2{\dot{M}}_{\star },$$

$$\omega \approx {\omega }_c=\frac{eB}{m_{e}c},$$

$$\Phi =\frac{P_{\mathrm{radio}}}{4\pi \Delta f{D}^2},$$

$$P_{\mathrm{radio}}\sim 2\times 10^{-3}\left(\frac{\pi B_{sw}^2{r}_m^2{V}_{sw}}{4\pi }\right),$$

$$r_m=R_p{\left(\frac{B_p^2}{8\pi P_{sw}}\right)}^{1/6},$$

$$\frac{\partial {\rho }_s}{\partial t}+\nabla \cdot \left({\rho }_s\mathbf{u}\right)={\mathcal{S}}_s-{\mathcal{L}}_s$$

$$\frac{\partial (\rho \mathbf{u})}{\partial t}+\nabla \cdot \left(\rho \mathbf{uu}+p\mathbf{I}+\frac{B^2}{2{\mu }_0}\mathbf{I}-\frac{1}{{\mu }_0}\mathbf{BB}\right)=\rho \mathbf{G}-\sum _{s=ions}{\rho }_s\sum _{n=neutrals}{\nu }_{sn}\mathbf{u}-\sum _{s=ions}{\mathcal{L}}_s\mathbf{u}$$

$$\frac{\partial }{\partial t}\left(\frac{\rho u^2}{2}+\frac{p}{\gamma -1}+\frac{B^2}{2{\mu }_0}\right)+\nabla \cdot \left[\left(\frac{\rho u^2}{2}+\frac{\gamma p}{\gamma -1}+\frac{B^2}{{\mu }_0}\right)\mathbf{u}-\frac{(\mathbf{B}\cdot \mathbf{u})\mathbf{B}}{{\mu }_0}+\mathbf{B}\times \frac{\nabla \times \mathbf{B}}{{\mu }_0^2{\sigma }_0}\right]$$

$$=\sum _{s=ions}\sum _{n=neutrals}\frac{{\rho }_s{\nu }_{sn}}{m_s+m_n}\left[3k\left(T_n-T_s\right)-m_s{u}^2\right]+\rho \mathbf{u}\cdot \mathbf{G}$$

$$+\frac{k}{\gamma -1}\sum _{s=ions}\frac{\left({\mathcal{S}}_s{T}_n-{\mathcal{L}}_s{T}_s\right)}{m_s}-\frac{1}{2}\sum _{s=ions}{\mathcal{L}}_s{u}^2+\frac{k}{\gamma -1}\frac{{\mathcal{S}}_e{T}_{n0}-{\mathcal{L}}_e{T}_e}{m_e}$$

$$\frac{\partial B}{\partial t}+\nabla \cdot \left(\mathbf{uB}-\mathbf{Bu}\right)=-\nabla \times \left(\frac{1}{{\mu }_0{\sigma }_0}\nabla \times \mathbf{B}\right)$$

$${\mathcal{S}}_s$$

$${\mathcal{L}}_s$$

$${\mathcal{S}}_e$$

$${\mathcal{L}}_e$$

$${\sigma }_0=\frac{n_e{e}^2}{m_e({\nu }_{ei}+{\nu }_{en})},$$

$${\nu }_{ei}=54.5\frac{n_i}{T_e^{3/2}},$$

$${\nu }_{en}$$

Full text

With the number of detected exoplanets now exceeding 3600 (Winn and Fabrycky, 2015), exoplanetary research has witnessed many remarkable advances recently. One of the most important areas in this field is the hunt for Earth-sized terrestrial planets residing in the habitable zone (HZ) of their host stars - the HZ represents the region within which a planet can support liquid water on its surface (Kopparapu et al., 2013); a probabilistic version of the HZ, encompassing a wide range of planetary and stellar parameters, has also been formulated (Zsom, 2015). The importance of this endeavour stems from the fact that such planets can be subjected to further scrutiny to potentially resolve the question as to whether they may actually harbour life (Cockell et al., 2016).

Most of the recent attention has focused on exoplanets in the HZ of M-dwarfs, i.e. low-mass stars that are much longer lived than the Sun, for the following reasons. Firstly, M-dwarfs are the most common type of stars within the Milky Way (Chabrier, 2003), implying that $\sim 10^{10}$ Earth-sized planets in the HZ of M-dwarfs may exist in our Galaxy (Dressing and Charbonneau, 2015). Second, owing to the HZ being much closer to such stars, it is much easier to detect exoplanets and characterize their atmospheres, if they do exist (Shields et al., 2016). Lastly, this field has witnessed two remarkable advances within the last year: the discovery of Proxima b (Anglada-Escudé et al., 2016) and the seven Earth-sized planets transiting the ultracool dwarf TRAPPIST-1 (Gillon et al., 2016, 2017). The significance of the former stems from the fact that it orbits the star closest to the Solar system, and the latter is important because there exist as many as three planets in the HZ with the possibility of life being seeded by panspermia (Lingam and Loeb, 2017a).

In light of these discoveries, the question of whether terrestrial exoplanets in the HZ of M-dwarfs are habitable is an important one (Shields et al., 2016). Amongst the many criteria identified for a planet to be habitable, the existence of an atmosphere has been posited as being crucial for surficial life-as-we-know-it (Lammer et al., 2009; Cockell et al., 2016). It is therefore evident that the study of atmospheric losses from exoplanets constitutes a crucial line of enquiry. Empirical and theoretical evidence from our own Solar system suggests that the erosion of the atmosphere by the stellar wind plays a crucial role, especially for Earth-sized planets where such losses constitute the dominant mechanism (Lammer, 2013; Brain et al., 2016) and the same could also be true for exoplanets around M-dwarfs (Dong et al., 2017b; Zahnle and Catling, 2017). Recent studies of atmospheric ion escape rates from Proxima b (and other M-dwarf exoplanets) also appear to indicate that the resulting ion losses are significant because of the extreme space weather conditions involved (Garraffo et al., 2016), potentially resulting in the atmosphere being depleted over a span ranging from tens to hundreds of millions of years (Dong et al., 2017b, a; Airapetian et al., 2017; Garcia-Sage et al., 2017).

Hence, in this paper we shall focus primarily on the atmospheric ion escape rates of the seven TRAPPIST-1 planets by adapting a sophisticated multi-species magnetohydrodynamics (MHD) model which self-consistently includes ionospheric chemistry and physics, and electromagnetic forces. In this work, we do not tackle the wide range of hydrodynamic escape mechanics that have been explored for terrestrial planets (Owen and Alvarez, 2016; Zahnle and Catling, 2017) on account of the above reasons.

The stellar wind of TRAPPIST-1

In order to commence our analysis of stellar wind-induced atmospheric loss, the stellar wind parameters of TRAPPIST-1 are required. Since the space weather conditions at the TRAPPIST-1 planets are presently unknown from observations, we must rely upon simulating the stellar wind of TRAPPIST-1. The latter is implemented by means of the Alfvén Wave Solar Model (AWSoM), a data-driven global MHD model that was originally developed for simulating the solar corona and solar wind (Sokolov et al., 2013; van der Holst et al., 2014). The AWSoM has been proven to be successful in reproducing high-fidelity solar corona conditions (van der Holst et al., 2014; Oran et al., 2013), and can readily be adapted to self-consistently model stellar wind profiles for a wide range of stars (Vidotto et al., 2013; Cohen et al., 2014; Garraffo et al., 2016). To adapt the AWSoM for modelling the TRAPPIST-1 stellar wind, we utilize the rotational period, radius, and mass of the star based on the latest estimates (Luger et al., 2017) and a mean magnetic field typical of similar late M-dwarfs (Morin et al., 2010). Further details concerning our approach can be found in Appendix A.

The steady state stellar wind solution is illustrated in Fig. 1. When compared to the normal solar wind solution (Jin et al., 2012), the stellar wind of TRAPPIST-1 is much faster ($\sim 3$ times) when evaluated at the same stellar distance. The critical surface is defined as the region where $v_{s}=v_{f}$, with $v_{s}$ and $v_{f}$ representing the stellar wind and fast magnetosonic speeds respectively. The surface occurs at distances of $\sim 30$ $R_{*}$ and $\sim 20$ $R_{*}$ for the fast and slow stellar wind; the latter duo originate at different regions on the star and have different speeds and densities (van der Holst et al., 2014). Hence, this leads to a very unique feature of the TRAPPIST-1 system, namely that part of the orbit of TRAPPIST-1b (the closest planet) lies within the critical surface, while all the other planets are embedded in the super-magnetosonic stellar wind.

This striking scenario does not exist within our own solar system primarily because of the proximity of TRAPPIST-1b to its host star (in conjunction with a strongly magnetized stellar wind). As TRAPPIST-1b orbits within the critical surface, the planet could magnetically interact with its host star directly. In turn, the star-planet interaction could perhaps (i) regulate the rotational rate (Pont, 2009), (ii) modify the properties of a local dynamo (Cuntz et al., 2000), and (iii) even give rise to a dynamo mechanism (Cébron and Hollerbach, 2014). In this context, we observe that variations in the magnetic field occur during the stellar cycle caused by the dynamo process. Thus, the distance of the critical surface is expected to also vary concomitantly, implying that TRAPPIST-1b could be subject to frequent transitions between sub-magnetosonic and super-magnetosonic stellar wind conditions along the lines of Proxima b (Garraffo et al., 2016).

Another distinguishing feature of the stellar wind from TRAPPIST-1 is its higher density. When combined with the higher wind speed, all of the planets are subjected to a much larger dynamic pressure compared to that experienced by the Earth. At the orbit of TRAPPIST-1b, the dynamic wind pressure is about $10^{3}$-$10^{4}$ times greater than the solar wind dynamic pressure at Earth. Even when we consider the furthermost planet, TRAPPIST-1h, the dynamic pressure is about $100$-$300$ times larger than the near-Earth environment. The existence of such an extreme wind pressure has already been documented for Proxima b (Garraffo et al., 2016), and its effects on the evolution of the planet’s magnetosphere have also been thoroughly investigated (Airapetian et al., 2017; Dong et al., 2017b). The ramifications of these extreme space weather conditions on the atmospheric ion escape rates of the TRAPPIST-1 planets are explored in the subsequent section.

We note that the mass-loss rate from TRAPPIST-1 is $\sim 2.6\times 10^{11}$ g/sec, which is about $10\%$ of the solar mass-loss rate. Although the density and velocity of the stellar wind are higher for TRAPPIST-1, the smaller size of the host star is responsible for yielding a value lower than that of the active young Sun (Airapetian and Usmanov, 2016). The mass-loss rate $\left(\sim 0.1\dot{M}_{\odot}\right)$ obtained for TRAPPIST-1 is broadly consistent with the upper bound of $\sim 0.2\dot{M}_{\odot}$ for the slightly larger star, Proxima Centauri (Wood et al., 2001).

Lastly, the stellar wind parameters provided in this paper are useful in determining the radio auroral emission from the TRAPPIST-1 planets, which can be used to constrain their magnetic fields (Grießmeier, 2015). In Appendix B, we show that the radio emission could potentially peak at $\mathcal{O}(0.1)$ MHz, and result in a radio flux density of $\mathcal{O}(0.1)$ mJy; the latter could be enhanced by $2$-$3$ orders of magnitude during a Coronal Mass Ejection (CME) event.

Ion escape rates for the TRAPPIST-1 planets

To simulate the ion escape rates for the seven planets of the TRAPPIST-1 system, we employ the sophisticated 3-D Block Adaptive Tree Solar-wind Roe Up-wind Scheme (BATS-R-US) multi-species MHD (MS-MHD) model that has been extensively tested and validated in the Solar system for Venus and Mars (Ma et al., 2004; Tóth et al., 2012; Ma et al., 2013; Jakosky et al., 2015), and was recently employed to study the atmospheric losses from Proxima b (Dong et al., 2017b). The reader is referred to these papers and the Appendix for further details concerning the details of the numerical implementation, the model equations, and the physical and chemical processes encoded within the model. Note that the neutral atmosphere is the source of the produced ions through, e.g., photoionization and charge exchange, and that only a small fraction of them will escape into space. Hence, the atmosphere will, in addition to being eroded, also undergo changes in the chemical composition (Dong et al., 2017b).

The input parameters of the model

We shall, for the most part, concern ourselves with describing and motivating our choice of the different input parameters required for the BATS-R-US model. Before proceeding further, we wish to caution the readers that many of the relevant planetary and stellar wind parameters of the TRAPPIST-1 system are unknown or poorly constrained. Hence, it is important to recognize that, on account of the many uncertainties involved, the ensuing escape rates may not necessarily be representative of the TRAPPIST-1 system.

The BATS-R-US model relies upon an atmospheric composition akin to Venus and Mars, implying that the TRAPPIST-1 planets are also assumed to possess a similar composition. There are several factors that must be noted in this context. Firstly, as seen from our Solar system, the ion escape rates for Venus, Mars and Earth are similar despite their compositions, sizes and magnetic field strengths being wildly dissimilar (Lammer, 2013; Brain et al., 2016), thereby indicating that the ion escape rates may be relatively sensitive to stellar wind parameters compared to planetary properties111The difference in stellar wind parameters at Venus, Earth and Mars is only up to a $\mathcal{O}(1)$ factor.; this is also partly borne out by the atmospheric ion escape rate calculations for Proxima b (Airapetian et al., 2017; Dong et al., 2017b). In addition, it has been shown recently that the ion escape rates are only weakly dependent on the surface pressure (Dong et al., 2017b). We observe that the inner planets of the TRAPPIST-1 system could have experienced significant losses of $H_{2}$ and water over fast timescales (Bolmont et al., 2017; Bourrier et al., 2017), leaving behind other atmospheric components. Lastly, a Venus-like atmosphere for the inner planets cannot be ruled out empirically, as noted in (de Wit et al., 2016).222In broader terms, gaining a thorough understanding of Venus-type exoplanets is highly relevant, because it allows us to compare and contrast their properties against exo-Earths (Angelo et al., 2017).

The next two input parameters to be specified are the surface pressure and the scale height for each of the planets. The former remains unknown at this stage, and we work with the fiducial value of $1$ atm at this stage. We anticipate that the surface pressure does not significantly alter the escape rates, at least for extreme stellar wind conditions, as demonstrated in (Dong et al., 2017b). The scale height $H_{x}$ is defined as

[TABLE]

where $g_{x}$ is the acceleration due to gravity for planet $X$. The latter quantity can be easily computed for all the planets since their masses and radii are known (Gillon et al., 2017; Luger et al., 2017).

The stellar wind parameters are obtained from the model described in the previous section. We must also prescribe the extreme ultraviolet (EUV) fluxes received at each of these planets, since the EUV flux plays an important role in regulating the extent of photoionization and the resultant stellar heating. This is accomplished by utilizing the values for the TRAPPIST-1 planets computed in (Bolmont et al., 2017; Wheatley et al., 2017; Bourrier et al., 2017).

The planetary magnetic field is a potentially important factor in regulating the ion escape rates. We consider the scenario where the planets are unmagnetized because this case yields an upper limit on the allowed escape rates (Dong et al., 2017b). Hence, if a planet was characterized by “low” escape rates in the unmagnetized limit, it would also typically possess low escape rates in the presence of a magnetic field. It must also be borne in mind that the planets orbiting TRAPPIST-1 are likely to be tidally locked, and it has been argued that such planets are likely to possess weak magnetic fields (Khodachenko et al., 2007). If a planet is weakly magnetized, it is likely that the total ion escape rate will be comparable to (but slightly lower than) the unmagnetized case (Dong et al., 2017b).

Results from the model

For each planet, we consider two limiting cases. The first corresponds to the scenario with maximum dynamic (and total) pressure over one orbit of the planet. The second corresponds to the case with minimum total pressure, but with the maximum magnetic pressure. The corresponding stellar wind parameters have been provided in the Appendix, and the escape rates in Table 1.

For all seven planets, the case with maximum total pressure yields a total atmospheric ion escape rate that is a few times higher than the corresponding case with minimum total pressure. The innermost trio of planets (‘b’, ‘c’ and ‘d’) have escape rates higher than $10^{27}$ sec*-1*, while the last four have rates lower than this value when the case with maximum total pressure is considered. In comparison, the escape rates for Mars, Venus and Earth are $\sim 10^{24}-10^{25}$ sec*-1* (Lammer, 2013), while that of Proxima b is $\sim 2\times 10^{27}$ (Dong et al., 2017b). In our subsequent analysis, we shall focus on this case (maximum total pressure) since it leads us towards determining the upper bounds on the escape rates.

Using the mixing-length formalism of (Canto and Raga, 1991) in conjunction with the definitions of the stellar and planetary mass-loss rates, it can be shown that

[TABLE]

where $N_{x}$ is the atmospheric escape rate arising from stellar wind stripping, $\dot{M}_{\star}$ is the stellar mass-loss rate, $R_{x}$ and $a$ are the radius and semi-major axis of the planet $X$ respectively (Zendejas et al., 2010; Lingam and Loeb, 2017b). We have normalized the escape rates for the TRAPPIST-1 planets in terms of the escape rate for TRAPPIST-1b, and compared the numerical and analytical predictions in Fig. 2.

An inspection of Fig. 2 reveals that the analytical formula is in good agreement with the numerical simulations, especially for the last four planets (‘e’ to ‘h’) which are regarded as being potentially capable of retaining atmospheres over Gyr timescales. Even for the inner planets, we find that the analytical results fall within the numerical values by a factor of $\lesssim 2$. Hence, Eq. 2 may facilitate a quick estimation of the atmospheric escape rates from unmagnetized planets. However, for both the simulations and the analytic results, it is important to note that the atmospheric escape rates were higher in the past, although quantitative estimates are difficult since the time-dependent stellar mass-loss rates are poorly constrained.

Fig. 3 depicts how the oxygen ions escape from two planets of the TRAPPIST-1 system (‘b’ and ‘g’) along with the interplanetary magnetic field lines. From the color contours, it is evident that the ion loss rate at TRAPPIST-1b is larger than TRAPPIST-1g. The solar wind magnetosonic Mach number of TRAPPIST-1b for the case with minimum total pressure (Ptot) equals $0.94$, therefore there is no shock formed in front of the planet (because of its sub-magnetosonic nature). This is particularly noteworthy since this condition is not prevalent in our Solar system for any of the eight planets. On the other hand, such a condition has been observed for moons located inside the planetary global magnetosphere (e.g., Ganymede and Titan). In contrast to the above case, the bow shock is present in all other cases and planets of the TRAPPIST-1 system.

Based on the background colors of the second row, minimum Ptot corresponds to the fast stellar wind and maximum Ptot signifies the slow stellar wind, which is fully consistent with the classical solar/stellar wind picture. The high Ptot manifested in the ‘slow’ stellar wind is mostly a consequence of the high stellar wind density. The second row of Fig. 3 also indicates that the stellar wind accelerates away from the star; the stellar wind velocity at TRAPPIST-1g is obviously larger than that of TRAPPIST-1b. These conclusions can also be verified through an inspection of Table 3 in the Appendix.

From Fig. 3, we observe that the plasma boundaries are compressed because of the fact that the total pressure is much higher at these planets compared to that experienced by Venus. The escape of ions is primarily driven by the pressure gradient and ${\bf J}\times{\bf B}$ forces in the momentum equation (Ma et al., 2013). Lastly, we refer the reader to Fig. 4 in the Appendix that depicts the ionospheric profiles of TRAPPIST-1g for the cases with minimum and maximum total pressure.

Discussion and Conclusions

We have arrived at an upper bound on the ion escape rates by considering the scenario where the planets are unmagnetized and subject to the maximum total (or dynamic) pressure. It was concluded that the innermost planet TRAPPIST-1b has an upper bound of $5.92\times 10^{27}$ sec*-1* while the corresponding value for the outermost planet TRAPPIST-1h is $1.29\times 10^{26}$ sec*-1*. As the planets are approximately Earth-sized and assumed to have a surface pressure of $1$ atm, we can estimate the timescales over which these planets can retain their atmospheres. The values range from $\mathcal{O}\left(10^{8}\right)$ years for TRAPPIST-1b to $\mathcal{O}\left(10^{10}\right)$ years for TRAPPIST-1h. Moreover, we also see from Fig. 2 that the overall escape rate declines monotonically as one moves outwards, from TRAPPIST-1b to TRAPPIST-1h. Hence, taken collectively, this may suggest that TRAPPIST-1h ought to be most “stable” planet amongst them, when viewed purely from the perspective of atmospheric ion loss. Along the same lines, it seems likely that TRAPPIST-1g will represent the best chance for a planet in the HZ of this planetary system to support a stable atmosphere over long periods.3333D climate simulations appear to suggest that TRAPPIST-1f and TRAPPIST-1g may not be amenable to surficial life as they enter a snowball state (Wolf, 2017), but the effects of tidal heating, which are expected to be considerable (Luger et al., 2017), were not included in the model.

At this stage, we must reiterate the caveats discussed earlier. Firstly, most of the planetary and stellar parameters are partly or wholly unknown since the appropriate observations are not currently existent. The presence of (i) a more massive atmosphere, (ii) a different atmospheric composition, and (iii) the planet’s magnetic field, are likely to alter the extent of atmospheric loss to some degree. Nonetheless, it can be surmised that the unmagnetized cases (with maximum total pressure) considered herein do yield robust upper bounds on the atmospheric ion escape rates. Apart from atmospheric loss, it is possible that outgassing processes could very well replenish the atmosphere (Kasting and Catling, 2003). Hence, resolving the existence of an atmosphere over Gyr timescales necessitates an in-depth understanding of the interplay between source and loss mechanisms. Lastly, it is important to note that stellar properties evolve over time, implying that the escape rates are also likely to change accordingly. In the case of M-dwarfs such as TRAPPIST-1, the pre-main-sequence phase is particularly long and intense, and expected to have an adverse impact on atmospheric losses (Luger and Barnes, 2015; Shields et al., 2016). The ensuing effect of XUV on hydrodynamic and ion escape rates during this phase has not been investigated in the paper.

Bearing these limitations in mind, let us now turn to a discussion of the implications. We have argued that TRAPPIST-1h and TRAPPIST-1g represent the most promising candidates in terms of retaining atmospheres over Gyr timescales. Instead, if the atmosphere were to be depleted over $\mathcal{O}\left(10^{8}\right)$ years, this could prove to be problematic for the origin of life (abiogenesis) on the planet although it must be acknowledged that the actual timescale for abiogenesis on Earth and other planets remains unknown (Spiegel and Turner, 2012). Abiogenesis has been argued to be accompanied by an increase in biological (e.g. genomic) complexity over time (Adami et al., 2000; Purvis and Hector, 2000) although this growth is not uniform, and may be contingent on environmental fluctuations. Hence, ceteris paribus, a planet capable of sustaining a stable atmosphere over long time periods (along with retaining a stable climate) might have a greater chance of hosting complex surficial organisms. The outer planets of the TRAPPIST-1 system may therefore lead to more diverse biospheres eventually.

We have also shown that the ionospheric profiles for the TRAPPIST-1 planets in the HZ are not sensitive to the stellar wind conditions at altitudes $\lesssim 200$ km (see Fig. 4 of the Appendix). This is an important result in light of the considerable variability and intensity of the stellar wind, since it suggests that the lower regions (such as the planetary surface) may remain mostly unaffected from under normal space weather conditions.444Stellar flares may have either a deleterious or beneficial effect on prebiotic chemistry that is dependent on a complex and interconnected set of factors (Dartnell, 2011; Ranjan et al., 2017).

Let us turn our attention to Table 1 for the seven Earth-sized exoplanets of the TRAPPIST-1 system. It is seen that the ion escape rate reduces as one moves outwards. Hence, for similar multi-planetary systems around low-mass stars, it may be more prudent to focus on the outward planet(s) in the HZ for detecting atmospheres since their escape rates could be lower. Similarly, when confronted with two planets with similar values of $R_{x}$ and $a$ (see Eq. 2), we propose that searches should focus on stars with lower mass-loss rates and stellar magnetic activity. Further details concerning the implications and generalizations of Eq. 2 will be reported in a forthcoming publication. Lastly, our results and implications are also broadly applicable to future planetary systems detected around M- and K-dwarfs endowed with similar features (Lingam and Loeb, 2017a).

To summarize, we have studied the atmospheric ion escape rates from the seven planets of the TRAPPIST-1 system by assuming a Venus-like composition. This was done by utilizing numerical models to compute the properties of the stellar wind and the escape rates, and the latter were shown to match the analytical predictions. We demonstrated that the outer planets of the TRAPPIST-1 system (most notable ‘h’ and ‘g’) are capable of retaining their atmospheres over Gyr timescales. However, as many factors remain unresolved at this stage, future missions such as the James Webb Space Telescope (JWST) will play a crucial role in constraining the atmospheres of the TRAPPIST-1 planets (Barstow and Irwin, 2016). In particular, a recent study concluded that spectral features for six of the seven TRAPPIST-1 planets could be detected with $<20$ transits with 5 $\sigma$ accuracy (Morley et al., 2017). Such observations would help constrain theoretical predictions, pave the way towards looking for biosignatures, and empirically estimate the putative habitability of these planets.

Acknowledgements

The authors thank A. Burrows, J.-F. Donati, S. Kane, A. Knoll, A. Loeb and J. Winn for the helpful discussions and comments. The support provided by the NASA Living With a Star Jack Eddy Postdoctoral Fellowship Program, administered by the University Corporation for Atmospheric Research is acknowledged. M.J. is supported by NASA’s SDO/AIA contract (NNG04EA00C) to LMSAL. Resources for this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The Space Weather Modeling Framework that contains the BATS-R-US code used in this study is publicly available from http://csem.engin.umich.edu/tools/swmf. For distribution of the model results used in this study, please contact the corresponding author.

Appendix A Modelling the stellar wind of TRAPPIST-1

In this Section, we describe the methodology involved in modelling the stellar wind of TRAPPIST-1 in more detail.

A.1 Description of the Alfvén Wave Solar Model

In this study, the TRAPPIST-1 stellar wind is simulated by means of the Alfvén Wave Solar Model (AWSoM; van der Holst et al., 2014), a data-driven global MHD model that had been initially developed for simulating the solar atmosphere and solar wind. The inner boundary condition of the magnetic field can be specified by different magnetic maps from available observations. The initial conditions for the stellar wind plasma are determined through the Parker solution, while the initial magnetic field is based on the Potential Field Source Surface (PFSS) model with the Finite Difference Iterative Potential Solver (FDIPS) described in Tóth et al. (2011).

Alfvén waves are driven at the inner boundary with a Poynting flux that scales with the surface magnetic field. The stellar wind is heated by Alfvén wave dissipation and accelerated by thermal and Alfvén wave pressure. Electron heat conduction (that includes both collisional and collisionless contributions) and radiative cooling are also included in the model. In the AWSoM, the electron and proton temperatures are treated separately, while the electrons and protons are assumed to have the same bulk velocity. However, heat conduction is applied only to the electrons, owing to their much higher thermal velocity. The system of governing equations is solved numerically using the Block Adaptive Tree Solar Wind Roe-type Upwind Scheme (BATS-R-US) code within the Space Weather Modeling Framework (Tóth et al., 2012).

By using a physically consistent treatment of wave reflection, dissipation, and heat partitioning between the electrons and protons, the AWSoM has successfully reproduced solar coronal conditions to a high degree of precision.

A.2 Application of the AWSoM to TRAPPIST-1

To adapt the AWSoM for the TRAPPIST-1, we modify the rotational mass, radius, and period of the star in accordance with the latest observational data (Luger et al., 2017). Hence, we specify $M_{\star}=0.08M_{\odot}$, $R_{\star}=0.11R_{\odot}$, and $P_{\star}=3.3$ days.

Due to the lack of direct surface magnetic field observations of TRAPPIST-1, we utilize a solar magnetogram under the solar minimum condition (GONG magnetogram of Carrington Rotation 2077; Jin et al., 2012) and scale the mean and radial magnetic field strength based on the magnetic field observations of similar late M-dwarfs (Morin et al., 2010). Based on the fact that the X-ray luminosity of TRAPPIST-1 is similar to that of the quiet Sun (Bourrier et al., 2017), we modify the Poynting flux parameter of the model such that the same amount of Poynting flux is generated as in the solar case. The simulation domain is extended to 250 $R_{\star}$ to ensure that the orbits of all seven planets in the TRAPPIST-1 system are duly encompassed.

Appendix B Deducing the magnetic fields of the TRAPPIST-1 planets

Since the past 40 years, it has been well known that radio emission from exoplanets can be used to detect them, and thereby determine their magnetic fields (Yantis et al., 1977). The basic idea behind detecting the planetary magnetic field ($B_{p}$) is that the cyclotron maser instability drives the emission of radio waves at a frequency approximately equal to the gyrofrequency (Zarka, 1998). The maximum emission frequency $\omega$ is given by

[TABLE]

where $m_{e}$ and $e$ are the electron mass and charge respectively. Hence, we propose that the radio auroral emission from the TRAPPIST-1 planets may be detectable by ground-based observatories when the frequency is above $\sim 10$ MHz; the Earth’s ionosphere reflects radio waves below this value thereby requiring space- or lunar-based telescopes (Grießmeier, 2015). As an illustrative example, let us suppose that $B_{p}\approx 0.1$ G for the TRAPPIST-1 planets, which leads us to $f=\omega/2\pi\approx 0.3$ MHz. This fiducial value is motivated by planetary dynamo scaling laws (Christensen, 2010), as well as the lower limit of the hypothesized magnetic field for Proxima b (Burkhart and Loeb, 2017).

Next, it is necessary to obtain a heuristic estimate of the radio flux density $\Phi$ at Earth to know whether the emission would be detectable. It can be estimated via

[TABLE]

where bandwidth $\Delta{f}\approx f/2$, $D$ is the distance to Earth and $P_{\mathrm{radio}}$ represents the planetary radio power (Zarka, 2007). Thus, a knowledge of $P_{\mathrm{radio}}$ and $\omega$ would suffice to determine $\Phi$. Although there exist several ways of computing $P_{\mathrm{radio}}$ (Grießmeier, 2015), it has been noted that the planetary radio power is dominated by the dissipation of the magnetic power carried by the stellar wind (Vidotto and Donati, 2017):

[TABLE]

where $B_{sw}$ is the interplanetary magnetic field (IMF) and $V_{sw}$ is the stellar wind velocity, and both of them have been listed in Table 3. Here, $r_{m}$ denotes the planetary magnetospheric radius and is given by

[TABLE]

where $P_{sw}\approx\rho_{sw}v_{sw}^{2}$ is the dynamic pressure and $\rho_{sw}$ is the stellar wind density provided in Table 3. Hence, it is possible to combine Eqs. (4), (5) and (6) to arrive at an estimate of $\Phi$. We find that the radio flux density, for the above choice of values, is $\mathcal{O}\left(10^{-4}\right)$ Jy. However, we wish to note that $\Phi$ can be boosted by a factor of $\sim 10^{2}-10^{3}$ during a strong Coronal Mass Ejection (CME) event, which would imply that it can attain a value of $\sim 10-100$ mJy. This follows from the fact that the stellar wind parameters are enhanced during a CME event, and the corresponding radio power is also increased accordingly (Grießmeier et al., 2007). We note that a similar approach has also been employed in the context of Proxima b, where it was shown that the radio flux density attained a peak value of $\sim 1-10$ Jy during a large Carrington-type CME event.

Thus, to summarize, future space-based (or lunar) low-frequency observations may be able to constrain the planetary magnetic fields of the TRAPPIST-1 planets by measuring the radio flux density and extrapolating backwards to determine the value of $B_{p}$.

Appendix C Atmospheric ion escape rates for the TRAPPIST-1 planets

In this Section, we briefly describe the workings of the code and provide additional results pertaining to the atmospheric ion escape rates.

C.1 Physical model and computational methodology

The 3-D Block Adaptive Tree Solar-wind Roe Upwind Scheme (BATS-R-US) multi-species MHD (MS-MHD) model was initially developed in the context of our Solar system i.e. for studying Mars (Ma et al., 2004) and Venus (Ma et al., 2013). We rely upon the code developed for Venus, and the neutral atmospheric profiles are based on the solar maximum conditions. The MS-MHD comprises of a separate continuity equation for each ion species, in conjunction with one momentum equation and one energy equation for the four ion fluids H+, O+, O${}_{2}^{+}$, CO${}_{2}^{+}$ (Ma et al., 2004, 2013). Unlike most global (Earth) magnetosphere models that commence from $2$-$3$ Earth radii, the Mars/Venus MS-MHD model contains a self-consistent ionosphere, and thus the lower boundary extends down to an altitude of 100 km above the planetary surface. The MS-MHD model, which serves as the basis of our paper, accounts for a diverse array of chemical processes, such as charge exchange, photoionization and electron recombination.

The various chemical reactions and their corresponding rate coefficients have been delineated in Table 2, and the reader may also consult Dong et al. (2017b) for further details. The densities of O+, O${}_{2}^{+}$, CO${}_{2}^{+}$ at the lower boundary satisfy the photochemical equilibrium condition as described in Dong et al. (2017b). The model also assumes that the plasma temperature (sum of ion and electron temperatures) is approximately double that of the neutral temperature at the lower boundary because of the high collision frequency.

The grid is also taken to be non-uniform and spherical in nature to accurately capture the multi-scale physics operating in different regions. Hence, the radial resolution ranges from around half the scale height at the lower boundary to several thousands of kilometers at the outer boundary. The horizontal resolution is chosen to be $3.0^{\circ}$ (in both longitude and latitude), while the simulation domain ranges from $-$45 R${}_{x}\leq$ X $\leq$ 15 Rx and $-$30 R${}_{x}\leq$ Y, Z $\leq$ 30 Rx, where Rx denotes the radius of planet $X$. The code is run in the Planet-Star-Orbital (PSO) coordinate system, where the $x$-axis is directed from the planet towards TRAPPIST-1, the $z$-axis is perpendicular to the planet’s orbital plane, and the $y$-axis completes the right-hand system.

The multi-species MHD equations are summarized as follows:

[TABLE]

where Eqs. (7), (8) and (9) represent the conservation of mass (of each species), momentum and energy respectively, whilst Eq. (10) is the magnetic induction equation. The total mass density $\rho=\sum_{s=ions}\rho_{s}$ and $\nu_{sn}$ is the elastic ion-neutral collision frequency (Ma et al., 2004, 2013). The ratio of the specific heats $\gamma$ is taken to be 5/3. The subscripts $s$, $n$ and $e$ indicate ion species $s$, neutral species $n$ and electron $e$. The other symbols have their usual definitions (Ma et al., 2004, 2013). In order to account for photoionization, we calculate the optical depth of the neutral atmosphere by applying the Chapman functions based on the numerical evaluation given by Smith and Smith (1972). The photoelectron gains an excess energy through the photoionization process, as indicated by the presence of $T_{n0}$ in Eq. (9). Therefore, we include the stellar heating via photoionization as those hydrodynamic models (e.g., Tripathi et al., 2015). In the above set of equations, note that the source ($\mathcal{S}$) and loss ($\mathcal{L}$) terms of species $s$ associated with photoionization ($\nu_{ph,s^{\prime}}$), charge exchange ($k_{is^{\prime}}$) and recombination ($\alpha_{R,s}$) are shown as follows:

[TABLE]

In Eq. (10), $\sigma_{0}$ is the electrical conductivity. In the model, the electrical conductivity in the planetary ionosphere is calculated using:

[TABLE]

where me is the electron mass, $\nu_{ei}$ and $\nu_{en}$ are the electron-ion and electron-neutral collision frequencies, respectively. The collision frequency $\nu_{ei}$ and $\nu_{en}$ are given by Ma et al. (2013)

[TABLE]

[TABLE]

The above set of equations is solved using an upwind finite-volume scheme based on an approximate Riemann solver, to ensure the appropriate conservation of plasma variables. To determine the steady-state solution, the simulation starts with a 2-stage local-time stepping scheme that enables different grid cells to select different advance time step thus accelerating convergence to the steady-state to save computational resources. Because of the stiffness of the source terms, a point implicit scheme is used for handling them.

C.2 Auxiliary results

In the main paper, the atmospheric ion escape rates were provided for two cases: (i) maximum dynamic and total pressure, and (ii) minimum total pressure, but maximum magnetic pressure. The corresponding stellar wind parameters for these two cases have been delineated in Table 3.

In Table 3, note that $N_{sw}$, $T_{sw}$ and $V_{sw}$ denote the number density, temperature and velocity of the stellar wind respectively at the locations of the seven planets, whereas IMF denotes the Interplanetary Magnetic Field. The case with minimum total pressure corresponds to the fast solar wind, which features a higher velocity albeit with a lower density, whereas the converse is true for the maximum total pressure that can be associated with the slow solar wind.

In Fig. 4, the ionospheric profiles for TRAPPIST-1g have been provided for the cases with minimum and maximum total pressure. An inspection of Table 3 reveals that the cases with minimum and maximum total pressure have very different stellar wind parameters, for e.g. density, velocity and interplanetary magnetic field. Despite the considerable variability in the stellar wind parameters, it is evident that the ionospheric profiles remain mostly unaffected. This would appear to indicate that the lower regions are effectively immune to the effects of the stellar wind, which lends some credibility to the fact that surface biological processes (if present) may not be significantly affected.