© ESO 2018

1 Introduction

TRAPPIST-1 planets recently discovered by Gillon et al. (2016, 2017) are the closest known transiting temperate Earth-sized exoplanets. The TRAPPIST-1 system hosts at least seven planets that are similar in size (from ~0.72 to ~1.13 R_⊕) and irradiation (from ~0.14 to ~4.3 S_⊕) to solar system rocky planets. Considering that the parent star TRAPPIST-1 is an ultra-cool (T_eff = 2550 K), low-mass (M_⋆ = 0.09 M_⊙) star, the planets of the system should have potentially followed evolutionary pathways very different from what the solar system planets experienced. They are therefore invaluable probes for comparative planetary science and habitability.

In a first approach, we can speculate on what TRAPPIST-1 planets might look like by comparing their size and irradiation with solar system planets. TRAPPIST-1b (0.64 S_Mercury) and TRAPPIST-1c (1.19 S_Venus) might be airless planets like Mercury, or endowed with a thick atmosphere like Venus. TRAPPIST-1d (1.14 S_⊕) is located nearthe inner edge of the habitable zone, traditionally defined as the range of orbital distances within which a planet can possibly maintain liquid water on its surface (Kasting et al. 1993; Kopparapu et al. 2013). Its capacity to host surface oceans, assuming an Earth-like atmosphere and water content, should depend on 1) its rotation mode and 2) subtle cloud albedo feedbacks (Yang et al. 2013; Kopparapu et al. 2016). TRAPPIST-1e (0.66 S_⊕) lies at the right distance to maintain surface liquid water if enough water is available and the atmosphere suitable. TRAPPIST-1f (0.89 S_Mars) and TRAPPIST-1g (0.60 S_Mars), being slightly bigger than Mars, could have retained a thick CO₂ atmosphere and thus conditions potentially similar to early Mars for a long period of time. Eventually, TRAPPIST-1h (0.30 S_Mars, 12 S_{Titan / Enceladus}) could look like a warm CH₄ /N₂-rich Titan-like planet, or a snowball icy-moon-like planet.

However, TRAPPIST-1 planets formed and evolved in a very different environment from our solar system planets. Each of the seven planets may potentially be airless today, as a result of the star extreme X/EUV irradiation (Wheatley et al. 2017; Bourrier et al. 2017) and stellar wind through history that could have blown away their atmosphere (Airapetian et al. 2017; Dong et al. 2017, 2018; Garcia-Sage et al. 2017). Moreover, during the first several hundreds of million years following their formation, while TRAPPIST-1 was a pre-main-sequence star and its luminosity significantly higher, eachof the seven planets could have faced a runaway phase in which the most condensable volatiles (e.g., water) would have been vaporized, and exposed to atmospheric escape. As much as several times the Earth’s ocean hydrogen content could have been lost in the process (Bolmont et al. 2017; Bourrier et al. 2017). The remaining oxygen could have therefore potentially built up in the atmosphere of the seven TRAPPIST-1 planets (Luger & Barnes 2015).

Transit timing variations (TTVs) measurements of TRAPPIST-1 planets (Gillon et al. 2017; Wang et al. 2017) suggest that planet bulk densities are compatible with terrestrial or volatile-rich composition. The latter possibility, and the fact that TRAPPIST-1 planets are in a near-resonant chain, suggest that the planets could have formed far from their star and migrated afterward to their current position. The planets may thus have been formed near or beyond the snowline and could have remained enriched in volatiles until now, despite a potentially massive early atmospheric escape (Bolmont et al. 2017). Uncertainties on the masses derived from TTVs are still affected by significant uncertainties but TTVs will eventually provide robust constraints on the density and volatile content.

Transit spectroscopy with the Hubble Space Telescope (HST) has been done on the two innermost planets, and suggest that they do not have a cloud/haze free H₂-dominated atmosphere (De Wit et al. 2016). This is somehow consistent with the fact that primordial H₂ enveloppes would have been exposed to efficient atmospheric escape on the small TRAPPIST-1 planets. At any rate, TRAPPIST-1 planets could still harbor a large variety of atmospheres, such as thick H₂O, CO₂, N₂, O₂, or CH₄ dominated atmospheres (see the review by Forget & Leconte 2014). In any case, each of these seven planets should be amenable to further characterization by the James Webb Space Telescope (JWST) as early as 2019 (Barstow & Irwin 2016; Morley et al. 2017).

The goal of the present study is to explore in more detail the possible climates of temperate-to-cold planets orbiting synchronously around cool stars in general, with a focus on the TRAPPIST-1 system (e,f,g,h). The constraints that we derive on their possible atmospheres could serve as a guideline to prepare future observations with JWST. We have explored in this work the conditions required for the coldest TRAPPIST-1 planets to prevent possible volatile species from atmospheric collapse, escape, or photodissociation. TRAPPIST-1 is a particular system where even weakly irradiated planets should likely be tidally locked. On synchronously rotating planets, the surface temperature of the cold points can be extremely low, making the different volatile species (N₂, CH₄, CO₂, etc.) highly sensitive to nightside trapping and potentially atmospheric collapse.

Conversely, we explore the stability of the same volatile species in the condensed phase (either icy or liquid) and on the surface, either on the dayside or the nightside. This condition, widely known for water as the runaway greenhouse limit, is extended here to other molecular species. Because these processes (runaway and collapse) are 3D on a synchronous planet, the most suited tools to explore them are 3D global climate models (GCM).

In Sect. 2, we describe our 3D global climate model, and more generally the physical parameterizations adopted in this work. In Sect. 3 we discuss the effect of tides on the rotation of TRAPPIST-1 planets. In the next sections, we explore the possible climates that can be expected on the four outer TRAPPIST-1 planets assuming that they are tidally locked and endowed with various volatiles: We discuss in Sect. 4 the ability of the three outer TRAPPIST-1 planets to sustain an atmosphere of background gases (N₂, CO, or O₂ ). Then, we explore whether we should expect oxidized CO₂-dominated atmosphere (in Sect. 5) or reduced CH₄-dominated atmosphere (in Sect. 6). Finally, in Sect. 7 we derive all the implications for the habitability of TRAPPIST-1 planets.

2 Method – the LMD generic global climate model (GCM)

The LMD generic model is a full three-dimension GCM which initially derives from the LMDz Earth (Hourdin et al. 2006) and Mars (Forget et al. 1999) GCMs. Since then, it has been extensively used to study a broad range of (exo)planetary atmospheres (Wordsworth et al. 2011, 2013, 2015; Forget et al. 2013; Charnay et al. 2013, 2015a,b; Leconte et al. 2013a,b; Bolmont et al. 2016; Turbet et al. 2016, 2017a,b).

Simulation input parameters include the observed characteristics of TRAPPIST-1 planets (Gillon et al. 2017; Luger et al. 2017), as summarized in Table 1. All simulations were performed assuming a circular orbit, a choice motivated by the small value of the maximum eccentricities derived from the stability of the system (Gillon et al. 2017; Luger et al. 2017). Even for a noncircular orbit, the orbital period is sufficiently small that the eccentricity should probably be quite high to significantly impact the climate of synchronous planets (see Bolmont et al. 2016 for their 10⁻⁴ L_sun case). We assumed that each of the planets is in synchronous rotation with 0° obliquity, as supported by calculations presented in Sect. 3.

The numerical simulations presented in this paper were all carried out at a horizontal resolution of 64 × 48 (e.g., 5.6° × 3.8°) in longitude × latitude. In all the simulations, the dynamical time step was set to 90 s. The physical parameterizations and the radiative transfer were calculated every 15 min and 1 h, respectively. Subgrid-scale dynamical processes (turbulent mixing and convection) were parameterized as in Forget et al. (2013) and Wordsworth et al. (2013). The planetary boundary layer was accounted for by the Mellor & Yamada (1982) and Galperin et al. (1988) time-dependent 2.5-level closure scheme, and complemented by a convective adjustment which rapidly mixes the atmosphere in the case of unstable temperature profiles. A filter was applied at high latitude to deal with the singularity in the grid at the pole (Forget et al. 1999). In the vertical direction, the model is composed of 26 distinct atmospheric layers that were built using hybrid σ coordinates and 18 soil layers. These 18 layers were designed to represent either a rocky ground (thermal inertia I_rock = 1000 J m⁻² K⁻¹ s${}^{-\frac{1}{2}}$), an icy ground (I_ice = 2000 J m⁻² K⁻¹ s${}^{-\frac{1}{2}}$) or an ocean (I_ocean = 20 000 J m⁻² K⁻¹ s${}^{-\frac{1}{2}}$ to take into account the efficient vertical mixing in the first tens of meter of the ocean, as previously done in Leconte et al. 2013a and Charnay et al. 2013) depending on the assumed surface. Since all the simulations were carried out for a synchronous rotation, thermal inertia should only affect the variability of the atmosphere. Oceanic heat transport was not included in this study.

The GCM includes an up-to-date generalized radiative transfer that takes into account the absorption and scattering by the atmosphere, the clouds and the surface from visible to far-infrared wavelengths, as described in Wordsworth et al. (2011). The radiative transfer is performed here for variable gaseous atmospheric compositions made of various cocktails of CO₂, CH₄, N₂, and H₂O, using the correlated-k method (Fu & Liou 1992; Eymet et al. 2016). Molecular absorption lines were taken from HITRAN 2012 (Rothman et al. 2013). Sublorentzian profiles (Perrin & Hartmann 1989; Campargue et al. 2012), collision induced absorptions (Gruszka & Borysow 1997; Baranov et al. 2004; Wordsworth et al. 2010a; Richard et al. 2012) and various other continua (Gruszka & Borysow 1997; Clough et al. 2005; Richard et al. 2012) were properly included in the calculations when needed. For the computation, we used between 32 and 38 spectral bands in the thermal infrared and between 36 and 41 spectral bands in the visible domain, depending on the atmospheric composition considered. Sixteen nonregularly spaced grid points were used for the g-space integration, where g is the cumulative distribution function of the absorption data for each band. We used a two-stream scheme (Toon et al. 1989) to take into account the scattering effects of the atmosphere and the clouds, using the method of Hansen & Travis (1974).

The emission spectrum of TRAPPIST-1 was computed using the synthetic BT-Settl spectrum¹ (Rajpurohit et al. 2013). We assumed a temperature of 2500 K, a surface gravity of 10³ m s⁻² and a metallicity of 0 dex.

The GCM directly computes the wavelength-dependent albedo of water ice and snow from a simplified albedo spectral law of ice and snow, calibrated to get ice and snow bolometric albedo of 0.55 around a Sun-like star, as in Turbet et al. (2016). Around TRAPPIST-1, we calculate that the average bolometric albedo for water ice and snow is ~ 0.21. Around an ultra-cool star like TRAPPIST-1, the bolometric albedo of water ice and snow is drastically reduced (Joshi & Haberle 2012; von Paris et al. 2013b; Shields et al. 2013) due to the shape of its reflectance spectrum (Warren & Wiscombe 1980; Warren 1984).

Melting, freezing, condensation, evaporation, sublimation, and precipitation of H₂O were included in the model. Water vapor was treated as a variable species in most of our simulations. This means that relative water vapor humidity is let free, but is always limited to 100% by our moist convective adjustment scheme. Similarly, we take into account the possible condensation and/or sublimation of CO₂ in the atmosphere (and on the surface) when needed but not the radiative effect of CO₂ ice clouds because their scattering greenhouse effect (Forget & Pierrehumbert 1997) should be low around cool stars like TRAPPIST-1 (Kitzmann 2017) and limited by partial cloud coverage (Forget et al. 2013). The effect of latent heat is properly taken into account when H₂O and/or CO₂ condense, evaporate, or sublimate.

Carbon dioxide and water cloud particle sizes were determined from the amount of condensed material and the number density of cloud condensation nuclei [CCN]. The latter parameter was taken to be constant everywhere in the atmosphere, and equal to 10⁶ kg⁻¹ for liquid water clouds, 10⁴ kg⁻¹ for water ice clouds (Leconte et al. 2013a) and 10⁵ kg⁻¹ for CO₂ ice clouds (Forget et al. 2013). Ice particles and liquid droplets were sedimented following a Stokes law described in Rossow (1978). H₂O precipitation are computed with the scheme from Boucher et al. (1995), with precipitation evaporation also taken into account.

All the numerical climate simulations were run long enough (up to 30 Earth years) to reach equilibrium. Simulations that lead to unstable CO₂ surface collapse were stopped when the rate of CO₂ surface condensation reached a positive constant, as in Turbet et al. (2017b). We note that more detail on the LMD generic model can be found in Forget et al. (1999, 2013); Wordsworth et al. (2011, 2013); Charnay et al. (2013); Leconte et al. (2013a), and Turbet et al. (2016, 2017a).

3 Effect of tides on TRAPPIST-1 planets

All observed TRAPPIST-1 planets are inside an orbital distance of 0.06 au. As a comparison, Mercury orbits at ~ 0.4 au from the Sun. For such close-in planets, tidal interactions are expected to be strong and to influence the orbital and rotational dynamics of the system.

Here we have used a standard equilibrium tide model (Mignard 1979; Hut 1981; Eggleton et al. 1998; Bolmont et al. 2011) to estimate the tidal evolution of the system. We combined an approach based on evolution timescale calculations and N-body simulations of the system using Mercury-T (Bolmont et al. 2015).

Mercury-T is an N-body code, which computes the orbital and rotational evolution of multi-planet systems taking into account tidal forces and their resulting torques (in the equilibrium tide framework), the correction for general relativity and the rotational flattening forces and torques. From this code, the evolution of the orbital parameters can be calculated (such as semi-major axis, eccentricity and inclination) as well as the rotation of the different bodies (i.e., the rotation period and obliquity of the planets). This code has previously been used to study the orbital dynamics of close-in and/or compact and/or near-resonant systems such as 55-Cnc (Bolmont et al. 2013), Kepler-62 (Bolmont et al. 2015) and Kepler-186 (Bolmont et al. 2014), and is now used here for the TRAPPIST-1 system.

3.1 Tidal dissipation and orders of magnitude

Simple order of magnitude calculations have allowed us to determine that the tide raised by the planets in the star is a priori negligible for this system today. Even considering a relatively high dissipation for a purely convective body (i.e., the dissipation of a hot Jupiter as estimated by Hansen 2010), we find semi-major axis and eccentricity evolution timescales of 10⁸ Myr and 5 × 10⁷ Myr respectively. We note that the dissipation is a measure of how fast the system is tidally evolving: the higher the dissipation, the faster the evolution. The age of TRAPPIST-1 has recently been estimated to be between 5 and 10 Gyr (Luger et al. 2017; Burgasser & Mamajek 2017). The evolution timescales for semi-major axis and eccentricity are thus consistent with Bolmont et al. (2011), which showed that the stellar-tide driven evolution around low mass stars and brown dwarfs was negligible for ages superior to ~ 100 Myr due to the decrease of the stellar radius.

The system therefore principally evolves due to the gravitational tide raised by the star in the planets (the planetary tide). The planetary tide mainly acts to decrease the obliquity of the planet, synchronize the rotation and on longer timescales decrease the eccentricity and semi-major axis. The dissipation in the planets depends on their internal structure and thermal state (Henning & Hurford 2014), as well as on the extension of the external fluid envelop (presence of surficial liquids – water ocean, magma ocean – and of a massive atmosphere) (e.g., Dermott 1979; Remus et al. 2015). On Earth, the dissipation is dominated by the dynamical response of the ocean and friction processes along the costline and to a lesser extent to interactions with seafloor topography in deep ocean (Egbert & Ray 2000, 2003; Williams et al. 2014). Dissipation on the Earth is highly dependent on the continent configuration and is therefore expected to significantly change on geological timescale asa consequence of tectonic plate motion (e.g., Poliakow 2005; Williams et al. 2014). The Earth’s dissipation is close to its highest value right now, and could have varied by a factor of almost ten over the last 200 Myr (Poliakow 2005).

In order to take into account the huge uncertainties in the dissipation factors of exoplanets (for which we do not know the internal structure), we consider various dissipation factors for the planets (from 0.1 to 10 times the dissipation of the Earth). The lowest value we consider here is roughly one order of magnitude larger than the dissipation estimated in Saturn (Lainey et al. 2017), representative of planets dominated by tidal response of a massivefluid envelop. The highest value is close to the maximal possible value and would be representative of very hot planets dominated by fluid-solid friction. There is no example in the solar system of such a dissipative object. Even the highly dissipative Jupiter’s moon Io (Lainey et al. 2009) has a dissipation function smaller than this extreme value, which is comparable to the Earth’s value (even if the dissipation process is very different). However, we could envision that Earth-sized bodies with a dissipation process comparable to that of Io could reach such a highly dissipative state. The tidal dissipation is also sensitive to the forcing frequency (e.g., Sotin et al. 2009; Henning & Hurford 2014), and therefore to the distance from the star. For simplicity, we have ignored this effect here and consider constant dissipation functions, independently of the distance from the star and the size of the planet, which is sufficient at first order to provide some typical tendencies.

Considering the dissipation for the planets of the system to be a tenth of the dissipation of the Earth (Neron de Surgy & Laskar 1997; Williams et al. 2014), comparable to the dissipation in Mars for instance (Yoder et al. 2003; Bills et al. 2005), we find evolution timescales for the rotation to range from 10⁻⁴ Myr for TRAPPIST-1b to 7 Myr for TRAPPIST-1h. For the obliquity, the evolution timescales range from 10⁻³ Myr for planet-b to 80 Myr for planet-h. Given the estimated age of the system, all planets are thus expected to have a small obliquity and to be near synchronization.

In the tidal framework we have used here, the rotation of the planets tend to pseudo-synchronization if the orbit is not circular. However, Makarov & Efroimsky (2013) showed that considering a more physical rheology forthe planet rather lead to a succession of spin-orbit resonance captures as the eccentricity of the planet decreases. We discuss the possibility of capture in spin-orbit resonant configuration in the following section.

3.2 Should we expect TRAPPIST-1 planets to be all tidally locked?

With such short period orbits, it is often assumed that bodily tides have spun-down the planets to the spin-orbit synchronous resonance ina relatively short time. However, it is now known that some other processes can sometimes act to avoid the synchronousstate. We will thus briefly review these processes. However, it appears that around such a low mass star, none of them is strong enough to counteract bodily tides so that all TRAPPIST-1 planets are probably in a synchronous-rotation state.

Indeed, one possibility for planets on an eccentric orbit is capture into a higher order spin-orbit resonance (Goldreich & Peale 1966). However, as discussed by Ribas et al. (2016) for the case of Proxima Centauri b, around a low mass star, the question is whether the dissipative tidal torque exerted by the star on the planet is strong enough to avoid capture into resonance (which is permited by the nonaxisymetric deformation of the planet). We used the methods detailed in Sect. 4.6 of Ribas et al. (2016) to calculate the probability of spin-orbit resonance capture of the planets of TRAPPIST-1. This method relies on comparing the tidal torque and the triaxiality torque, which depend strongly on eccentricity. The lower the eccentricity, the lower the spin-orbit resonance capture probability. For capture to be possible, the eccentricity of a given planet in the system would need to be roughly above 0.01. Probability of capture becomes greater than 10% only for an eccentricity greater than 0.03. However, simulations of the system dynamics accounting for tides and planet-planet interactions (see below) seem to show that such eccentricities are on the very high end of the possible scenarios. The spin orbit capture is thus seen as rather improbable in such a compact system.

The other possibility is that thermal tides in the atmosphere can create a strong enough torque to balance the stellar tidal torque on the mantle, as is expected to be the case on Venus (Leconte et al. 2015; Auclair-Desrotour et al. 2017). For this process to be efficient, the planet must be close enough from the star so that tides in general are able to affect the planetary spin, but far enough so that bodily tides are not strong enough to overpower atmospheric tides. In a system around such a low-mass star, this zone rests well beyond the position of the seven discovered planets (see Fig. 3 of Leconte et al. 2015). Atmospheric tides are thus unable to significantly affect the spin of the planet.

Fig. 1 Eccentricity (top panels) and tidal heat flux (bottom panels) for the three inner planets of TRAPPIST-1 for different tidal dissipation factors (different colors): from 0.1 to ten times the Earth’s value (taken from Neron de Surgy & Laskar 1997). Stars indicate the mean values.

Fig. 2 Same as Fig. 1 but for the four outer planets.

3.3 Tidal N-body simulations

We then performed N-body simulations using Mercury-T (Bolmont et al. 2015) to compute the complete evolution of the system, taking into account tides, general relativity and the rotational flattening of the different rotating bodies. We explored the range of dissipation factors discussed above (from 0.1 to ten times Earth dissipation). Figures 1 and 2 show the evolution of the eccentricity and resulting tidal heat flux for the different planets of TRAPPIST-1 and for the different dissipation factors. Table 2 summarizes the mean values of eccentricities and tidal heat fluxes for all the planets. Due to computation time considerations, simulations corresponding to a dissipation of 0.01 Earth dissipation and lower did not reach equilibrium. In the state of equilibrium, one should observe the following trend: the higher the tidal dissipation, the lower the equilibrium eccentricity, but the higher the tidal heat flux (see Bolmont et al. 2013 where this was discussed for 55 Cancri).

The initial state of our simulations corresponds to the orbital state of the system determined in Gillon et al. (2017) for planets b to g, and we used Luger et al. (2017) for the orbital parameters of planet h. As the evolution timescales of rotation and obliquity are small compared to the estimated age of the system, we considered the planets to be initially in synchronization and with a very small obliquity. We considered two sets of initial eccentricities: all eccentricities at 10⁻⁶ and the eccentricities derived from TTVs of Gillon et al. (2017).

All simulations display the same behavior: after a short initial phase of eccentricity excitation, all excentricities decrease on a timescale depending on the dissipation factor, to reach a mean equilibrium value. This equilibrium value is the result of the competition between tidal damping and planet-planet excitations (e.g., Bolmont et al. 2013) and the eccentricity oscillates around it. The eccentricities corresponding to the equilibrium value are all relatively small as they are inferior to 10⁻³ for planets b and c, and inferior to 10⁻² for planets d to h. The equilibrium value depends slightly on the dissipation factor of the planets: the higher the dissipation factor, the smaller the eccentricity. We note that such small eccentricities would have no effect on the climate (e.g., Bolmont et al. 2016); this is why we have assumed circular orbits for all planets in the climate simulations performed in this study.

The same kind of behavior can be seen for the obliquity of the planets: they assume an equilibrium value, result of the competition between tidal damping and planet-planet excitations. The equilibrium values are very small. For instance, the obliquities of the planets are smaller than 1°, which is why here we also assumed a zero obliquity for all planets.

Estimates of internal heat fluxes

Due to planet-planet interactions, the eccentricities and obliquities of the planets are not zero. This means that the planets are constantly submitted to a changing potential and constantly being deformed. This implies that the planets become tidally heated.

Our simulations with Mercury-T allow us to derive a possible state of the system and the corresponding tidal heat flux for all planets (see Fig 3). We find that the equilibrium eccentricity is enough to create a significant heat flux for the inner planets. For instance, assuming the tidal dissipation of the Earth for all the TRAPPIST-1 planets, we find that the eccentricity of planet b varies from ~ 8 × 10⁻⁶ to 1.5 × 10⁻³, with a mean value at 6 × 10⁻⁴. These eccentricities lead to a tidal heat flux which varies from ~ 0.02 W m⁻² to ~ 25 W m⁻², with a mean value at ~ 4 W m⁻². Table 2 shows mean values for the tidal heat flux for all the planets, for different scenarios of dissipation. Figures 1 and 2 show a snapshot of the evolution of the eccentricity and tidal heat flux over 1 Myr for each planet and for three different tidal dissipation factors.

We warn the readers that the mechanism of electromagnetic induction heating recently proposed by Kislyakova et al. (2017) should have a negligible contribution to the surface heat flux. We calculated from Kislyakova et al. (2017; Table 1) that the induction heating should not produce more than 8, 19, 8, 0.6, 0.2, 0.08, and 0.03 mW m⁻² for TRAPPIST-1b, c, d, e, f, g, and h, respectively.

For the four outer planets of the TRAPPIST-1 system, the mean tidal heat fluxes derived from N-body calculations are lower than 0.09 W m⁻², which corresponds roughly to the Earth mean geothermal heat flux (Davies & Davies 2010). Therefore, tidal heat flux is expected to play a minor role on the climate of these planets. It could nonetheless contribute significantly to:

• 1.

  Surface temperature of the cold traps, and hence to the atmospheric condensation of background gas (like N₂).

• 2.

  Maximum amount of the various volatiles that could be trapped on the cold points of the planets.

• 3.

  Maximum depth of a subsurface liquid water ocean (Luger et al. 2017).

• 4.

  More generally, the internal structure and the orbital dynamics of the planets.

The two first effects are explored in the following sections.

Fig. 3 Tidal heat flux map (in W m−2) as a function of semi-major axis (X axis) and eccentricity (Y axis) for planets with the Earth mass and radius. The dissipation efficiency is assumed to be one tenth of the Earth one to account for the dissipation in the mantle only. The tidal flux scales linearly with this parameter. This value can, however, easily change by orders of magnitude with the internal structure of each planet. This map should thus serve as a rough guide only.

Fig. 4 Equilibrium vapor pressures as a function of temperature for nine different species (where experimental data are available) that could be abundant on TRAPPIST-1 planets. Solid black lines where superimposed to indicate the equilibrium temperatures of TRAPPIST-1e,f,g, and h (assuming a surface albedo of 0.2). Stars indicate the positions of the triple points. These curves were adapted from Fray & Schmitt (2009).

4 Could TRAPPIST-1 planets be airless planets?

Leaving aside the case of H₂/He-rich atmospheres, we focus here specifically on the case of the next three most volatile species: N₂, CO, and O₂, because they are the best compromise between volatility (see Fig. 4) and abundance. The cases of CO₂ and CH₄, which are significantly less volatile, are discussed later on in Sects. 5 and 6. The main goal of this section is to assess the necessary conditions for the three outer TRAPPIST-1 planets to sustain a global, background atmosphere (i.e., a rather transparent atmosphere that can ensure the transport of heat and the pressure broadening of absorption lines of greenhouse gases). Background gases are essential because they can prevent the more volatile species such as CO₂ or NH₃ from collapsing on the nightside.

We haveassumed for now that the surface is covered by water (liquid or icy) because it is expected to be the most abundant volatile, as well as the less dense (this is a key element for planetary differentiation) and most condensable (see Fig. 4). We refer the reader to the review by Forget & Leconte (2014; and the references therein) for more information on possible sources and sinks of these volatile species.

4.1 Can a global atmosphere avoid atmospheric collapse?

To begin, we assumed that the three outer TRAPPIST-1 planets initially start with an atmosphere. On synchronously rotating planets, the nightside surface temperature can be so low that the atmosphere itself starts to condense on the surface. We looked for the minimal atmospheric pressure necessary to prevent them from atmospheric collapse, a configuration for which all the volatiles are permanently frozen on the nightside.

For this, we performed several simulations of TRAPPIST-1f, g, and h planets (surface albedo fixed to 0.2 corresponding to a water ice surface around TRAPPIST-1, or coincidentally to a rocky surface) endowed with a pure N₂ atmosphere (with H₂O as a variable species) for various atmospheric pressures (from 1 bar down to 10 millibars). Surface temperature maps corresponding to these experiments are shown in Fig. 5.

We find that a pure N₂ atmosphere (with H₂O as a variablegas) is quite resistant to atmospheric collapse for each of the three TRAPPIST-1 (fgh) outer planets. A collapse would be expected for N₂ partial pressure (pN₂) slightly lower than 10 millibars, and this value should hold for each of the three planets notwithstanding their various levels of irradiation. Our simulations indicate in fact (see Fig. 5) that if TRAPPIST-1h is always globally colder than TRAPPIST-1g (which is globally colder than planet f), it is not necessarily the case for the temperature of their cold points. TRAPPIST-1f, g and h planets have rotation periods of approximately 10¹ Earth days and they lie thus near the transition between slow and fast rotating regimes (Edson et al. 2011; Carone et al. 2015, 2016). They should be in one of these two regimes and could potentially be in both, depending on the initial forcing (Edson et al. 2011). Since the temperature of the cold points is critically dependant on the circulation regime (see Carone et al. 2016, their Figs. 1–3), it is difficult to assess which of these three TRAPPIST-1 planets should be more sensitive to atmospheric collapse.

Similarly to N₂, CO and O₂ are rather transparent in the infrared region of the surface thermal emission (between 10 and 100 microns, here) and have a similar molar mass to that of N₂ (between 28 and 32 g mol⁻¹). We can then safely extend our results for N₂ to CO and O₂-dominated atmospheres. These two gases are slightly more condensable gas and are thus expected to collapse for atmospheric pressure slightly higher than 10 millibars (see the legend of Fig. 5). These results could be tested in future studies with models that would properly take into account the radiative properties of a CO or O₂ -dominated atmosphere, and that would explore the sensitivity of these results to the assumption made on the surface composition (water, rock, etc).

More generally, we note that as the dominant gas becomes less and less volatile, building up an atmosphere becomes more and more complicated due to atmospheric collapse. At any rate, such collapse would trigger a positive feedback, because as the atmosphere condenses, the heat redistribution would become less efficient, leading to even more condensation. This would drive the planets to a complete and irreversible atmospheric collapse.

Fig. 5 Maps of surface temperatures (averaged over 50 Earth days) for TRAPPIST-1f, g and h, assuming initially cold, water ice covered frozen planets, endowed with a pure N2 atmosphere (with H2O as a variable gas) at three different surface pressure (10 millibars, 100 millibars, 1 bar). Blue colored label indicates the mininum temperature reached by the coldest point of the planet throughout the entire simulation. As a reminder, for partial pressures of 10 millibars (100 millibars and 1 bar), N2 is expected to collapse at 53 K (62 K and 79 K), CO at 56 K (66 K and 83 K), and O2 at 61 K (72 K and 87 K). The geothermal heat flux is not taken into account, but the coldest temperatures found in these cases a posteriori show that it can be neglected.

4.2 How much volatile can be trapped on the nightside of an airless planet?

Conversely, we suppose now that the planets initially start without a global atmosphere, which could have been blown away during the early active phase of TRAPPIST-1. In this configuration, all the volatiles (accreted, outgassed, or residual from the initial collapse) are expected to accumulate on the cold side of the planet. We have calculated here the maximum amount of volatiles that could be trapped in ice caps before a global atmosphere would be (re)formed.

The nightside surface temperature T_night on an airless tidally locked planet is determined by the geothermal heat flux F_geo of the planet: $$T_{\text{night}}={\left(\frac{F_{\text{geo}}}{\sigma }\right)}^{\frac{1}{4}}.$$(1)

For geothermal heat flux of 500 mW m⁻² (corresponding to planets with strong tidal dissipation in their interior) (or 50 mW m⁻², respectively), we get a nightside surface temperature of 50 K (30 K, respectively). The temperature at the base of the nightside ice cap (e.g., the temperature below the volatiles) depends on the geothermal heat flux F_geo and the nightside surface temperature T_night. When the ice cap is full (e.g., when ices start to convert into liquids), the temperature at the base of the nightside ice cap should be close (always higher, though) to that of the triple point. The triple point temperature is equal to 63K for N₂, 68K for CO, and 55K for O₂.

At these temperatures that are slightly warmer than those expected at the surface of Pluto (Forget et al. 2017), the viscosity of ices can be rather low. For instance, we estimate from Umurhan et al. (2017, Eq. (7)) that the viscosity of N₂ ice at 45 K (or 52 and 60 K, respectively) should be roughly 1.6 × 10¹⁰ Pa s (8 × 10⁸ and 7 × 10⁷ Pa s, respectively). With this condition in mind, it is not clear whether the maximum size of glaciers – formed by the accumulation of volatiles – should be controlled by the basal melting condition or by the glacial flow. Assessing this question properly would require to compare the efficiency of the glacial flow with the rate at which and the position where gaseous N₂ would condense on the nightside.

We assess below the case of the basal melting condition as it gives us an upper limit on the maximum amount of volatile possibly trapped on the nightside. When the nightside glaciers start to melt at their base, the ice flow should accelerate and expand significantly on the dayside of the planet. Not only basal melting is expected to be a very efficient process to transport ices (e.g., N₂, CO, O₂, etc.) from the nightside to the dayside, but the ices that reach the terminator should get sublimed and transport latent heat from the dayside to the nightside. This positive feedback would drive the planet into a runaway process, resulting in the formation of a new, global atmosphere.

For any species (we arbitrarily chose N₂ here), and for any of the seven TRAPPIST-1 planets, we derive from the basal melting condition the following set of two equations: $$\{\begin{array}{l}{T}_{\text{base,liq}} = T_{\text{ref}} {\text{e}}^{\frac{\left(\frac{1}{{\rho }_{\text{liq}}}-\frac{1}{{\rho }_{\text{ice}}}\right)}{L_{\text{fus}}} (g{\rho }_{\text{ice}}{h}_{\text{max}}+P_{{\text{N}}_2}-P_{\text{ref}})} \text{,}\hfill \\ h_{\text{max}} = \frac{{\lambda }_{{\text{N}}_2} (T_{\text{base,liq}}-T_{\text{night}})}{F_{\text{geo}}} \text{,}\hfill \end{array}$$(2)

with T_base,liq the temperature at the bottom of the glacier, ${\lambda }_{{\text{N}}_2}$ the conductivity of N₂ ice, ρ_liq and ρ_ice the volumetric mass densities of liquid and icy N₂, L_fus the latent heat of N₂ ice melting, P_ref and T_ref the pressure and temperature of the triple point of N₂. $P_{{\text{N}}_2}$ is the partial pressure of N₂ calculated at saturation from the Clausius-Clapeyron relationship (see Fig. 4) at the surface nightside temperature T_night.

This set of two equations corresponds respectively to:

• 1.

  the solid-liquid thermodynamical equilibrium at the base of the glacier. We note that the pressure at the bottom of the glacier is controlled by the weight of the glacier (and marginally, by the atmospheric pressure);

• 2.

  the geothermal gradient. It is assumed that the temperature inside the glacier rises linearly with depth, with a lapse rate fixed by the internal heat flux (conductive regime).

This set of equations can be solved explicitly after several variable changes and using the Lambert W function, defined as the solution of X = e^X, as done in Turbet et al. (2017b).

We calculated and plotted in Fig. 6 the nightside maximum thickness as a function of geothermal heat flux for various ices. For geothermal heat flux ranging from 50 to 500 mW m⁻², the maximum thicknesses range:

• 1.

  from 200 to 5 meters (Global Equivalent Pressure – GEP – from 15 to 0.2 bar) for N₂;

• 2.

  from 300 to 10 m (GEP from 15 to 0.5 bar) for CO;

• 3.

  from 200 to 0 m (GEP from 10 to 0 bar) for O₂. O₂ has the lowest triple point temperature (see Fig. 4).

We note that these values are of the same order of magnitude than in the atmosphere of Venus (~ 3 bars of N₂ ), Earth (0.78 bar of N₂; 0.21 bar of O₂) and Titan (~1.5 bar of N₂), the only solar system rocky bodies that were able to sustain a thick, global atmosphere.

For geothermal heat flux roughly lower than ~ 5 × 10² mW m⁻², there is a strong hysteresis on the initial state (volatiles in the atmosphere; or volatiles condensed at the cold points). Planets that initially lost their atmosphere could stably accumulate quantities of N₂ /O₂/CO up to the equivalent of few bars, in condensed form on the surface of their nightside. If somehow this scenario occured (through massive accretion or outgassing), the volatiles could not be retained on the nightside. The planet would suddenly sublime the entire volatile content of N₂, CO, O₂, CH₄, and so on. forming a brand new, global atmosphere.

Extreme events such as large meteoritic impact events could also have the potential to destabilize the volatiles that have been trapped in condensed form on the surface of the nightside. Only very eccentric bodies orbiting in the TRAPPIST-1 system could hit the planets near the anti-substellar point and potentially sublime the volatiles that should be preferentially trapped there. In the solar system, it has, for instance, been proposed that the observed distribution of impact craters on Mercury could be explained by a large number of very eccentric bodies that would have hit Mercury near the substellar and anti-substellar regions, while the planet was in synchronous rotation (Wieczorek et al. 2012). We note that there should be generally a large proportion of high-eccentricity bodies in the vicinity of the star, favoring subsequently impact events in the anti-substellar region.

For geothermal heat flux roughly higher than ~ 5 × 10² mW m⁻², we find that the planets could easily form an atmosphere even with very low amount of volatiles. In fact, we start to reach a regime here where the geothermal heat flux itself could significantly contribute to limit the atmospheric collapse as discussed in the previous section. We note that, at such high geothermal heat flux, heat could – and should – be transported by convection; this would significantly alter the calculations made here.

As shown in the previous section, eccentricities of TRAPPIST-1 planets are expected to vary with time, and tidal dissipation and surface heating with it, on timescales of around one Earth year (see Luger et al. 2017; their Supplementary Fig. 6). Peaks of tidal surface heating could trigger the destabilization of volatiles trapped on the nightside, although we would expect a delay and smoothing depending on where the tidal dissipation occurs. Actually, the heat itself could alter the internal structure. Detailed calculations of time-dependant tidal-induced surface heat flux (and implications) could be assessed in future studies.

As previously suggested in Turbet et al. (2016), the large-scale gravitational anomalies on tidally locked planets could be aligned with the star-planet axis. This means for instance that it is likely that a large basin (for example impact-induced) could be present at the anti-substellar point of TRAPPIST-1 planets. This could – in the same fashion than nitrogen ice is trapped on Pluto, in Sputnik Planum (Bertrand & Forget 2016) – significantly increase the amount of volatiles possibly trapped at the cold point of the planets. Furthermore, the weight of the ices trapped on the nightside could cause the underlying water ice shell to slump, creating by itself (or amplifying the size of) an anti-substellar basin. Such a process has recently been proposed as one of the possible scenarios to explain the formation of Sputnik Planitia on Pluto (Hamilton et al. 2016).

Fig. 6 Maximum nightside thickness of various types of ice (N2, CO, O2, CH4 and CO2) – assumingthat it is limited by basal melting – as a function of the geothermal heat flux. It is assumed here that the entire atmosphere has collapsed at the cold points of the planet and that the surface temperature at the top of the glacier is controlled by the geothermal heat flux. As a reference, vertical dashed lines indicate the average surface tidal heat flux on TRAPPIST-1 planets derived from Table 2, for a tidal dissipation equal to that of the Earth. We also added (vertical solid gray line) the average geothermal heat flux on Earth, to give the reader a rough sense of the amplitude of the radiogenic heating on TRAPPIST-1 planets. These quantities can be numerically converted in term of global equivalent surface pressure when multiplied by a factor $\frac{\rho g}{2}$. The thermodynamical and rheological properties of the ices were taken from https://encyclopedia.airliquide.com; http://webbook.nist.gov; Roder (1978); Schmitt et al. (1997); Fray & Schmitt (2009); Trowbridge et al. (2016), and Umurhan et al. (2017). Missing rheological data were mimicked on N2.

4.3 Residual atmospheres

Even though the atmosphere may have collapsed on the cold side of the planet, a residual, thin atmosphere could remain. The volatiles trapped on the nightside should be in fact in thermodynamical equilibrium with a residual atmosphere whose thickness depends on the surface temperature of the nightside, and on the type of volatiles trapped (assuming that the reservoir of volatiles is large enough).

For a geothermal heat flux of 100 mW m⁻² (or 200 and 400 mW m⁻², respectively), the temperature of the cold side is ~36 K (42 and 50 K, respectively) and the remnant atmosphere could be as thick as ~ 0.6 Pa of N₂ (18 and 400 Pa, respectively), 7 × 10⁻² Pa of CO (3 and 110 Pa, respectively), and 7 × 10⁻³ Pa of O₂ (0.6 and 30 Pa, respectively). For the other volatiles (CH₄ and CO₂, for example), the thickness (or surface pressure) of a residual atmosphere would be several orders of magnitude lower.

Such residual atmospheres should not be thick enough to significantly increase the global heat redistribution (or possibly the greenhouse effect) and trigger subsequently a N₂, CO or O₂ runaway process. We remind the reader that the minimum atmospheric pressure required to sustain a global atmosphere is ~ 10³ Pa (see Sect. 4.1).

Even though detecting a residual atmosphere of N₂, CO, O₂, etc. might be extremely challenging, as shown above such measurements could tell us a lot about 1) the temperature of the nightside and thus the internal heat flux of the planet and 2) the composition of the nightside reservoir of volatiles.

We also note that volatiles possibly trapped on the nightside of airless close-in planets would form a residual atmosphere that would be exposed to various processes of atmospheric escape (mainly stellar-wind sputtering and X/EUV-driven hydrodynamic escape). This indicates that volatiles trapped on the nightside of geothermally active tidally-locked planets might not be protected from atmospheric escape.

5 CO₂-dominated atmospheres

All the solar system terrestrial planets are either airless bodies (e.g., Mercury) or worlds where CO₂ is – or was – abundant in the atmosphere (e.g., Venus, Mars) and/or in the subsurface (e.g., Earth). We assume in this section that the four TRAPPIST-1 outer planets possess today large quantities of CO₂ either in their atmosphere, on their surface or in their subsurface, and we explore the possible implications.

5.1 Stability of a CO₂-dominated atmosphere

Carbon dioxde is much more condensable than any other species discussed in the previous section, as illustrated in Fig. 4. On synchronously rotating planets, the nightside surface temperature can be extremely low, leading to the condensation of gaseous CO₂ on the surface. This could potentially prevent TRAPPIST-1 planets from building up thick CO₂ atmospheres.

To test this idea, we performed 130 3D climate numerical simulations of the four TRAPPIST-1 outer planets (surface albedo fixed to 0.2) for atmospheres made of various mixtures of N₂ and CO₂. In the same vein as Turbet et al. (2017b), we find that depending on the partial pressure of background gas (N₂, here) and on the partial pressure of CO₂, the gaseous CO₂ might condense or not, as shown in Fig. 7. The shape of the diagrams is controlled by various physical processes:

• 1.

  The higher the background gas content is, the more efficient the heat redistribution is. This tends to increase the temperature of the cold points and limit the CO₂ condensation. High background gas content also favor the pressure broadening of CO₂ absorption lines, which increases the greenhouse effect of the atmosphere.

• 2.

  The higher the CO₂ content is, the higher its greenhouse effect is, but the higher its condensation temperature is. These two processes are in competition with each other, as illustrated in Soto et al. 2015 (in their Fig. 1).

Figure 7 shows a bistability in the CO₂ atmospheric content. If the planet initially starts with a thick CO₂ atmosphere (e.g., ten bars), the greenhouse effect and the heat redistribution are efficient enough for such atmosphere to be stable. Conversely, if the planet initially starts with a low CO₂ atmospheric content or no CO₂ at all and progressively accumulates somehow additional CO₂ in the atmosphere (e.g., by volcanic outgassing), all the extra CO₂ should keep condensing on the nightside. The planet would thus be permanently locked with a cold, thin CO₂ atmosphere.

We can, for instance, see in Fig. 7 that a background atmosphere of ~ 100 millibars of N₂ is not sufficient to build up a CO₂-rich atmosphere from scratch – on any of the four outer planets – due to the nightside surface condensation of CO₂. We can also see that TRAPPIST-1h is unable to sustain a dense CO₂ atmosphere (>100 mbar) even with several bars of N₂. For TRAPPIST-1e, f and g, if the initial CO₂ content is – for a given amount of background gas – below the “unstable” dotted line, then the planets are unable to build up CO₂-rich atmospheres (blue color). However, if the same planets start with an initial CO₂ content higher than this limit, CO₂ thick atmospheres are found stable (red color).

5.2 The fate of surface condensed CO₂

When CO₂ starts to condense on the nightside, we see from Turbet et al. (2017b), there should be two processes that control the maximum amount of CO₂ possibly trapped on the cold side of the planets. Firstly CO₂ ice flow from the nightside to regions of sublimation, and secondly gravitational burial of CO₂ ice beneath water ice cover due to its higher density.

Fig. 7 Climate regimes reached as a function of the partial pressures of N2 and CO2. For each set of (pN2, pCO2 ), it is indicated if the atmosphere is stable (red) or not (blue) to the atmospheric condensation/collapse of CO2. The black arrows indicate how planets that have an unstable atmosphere would evolve on this diagram. Temperatures (in green) correspond to the rough estimate of the temperature of the cold point, at the stable lower boundary (blue is up; red is down). Simulations were performed assuming a surface albedo of 0.2 (corresponding both to a water ice surface around TRAPPIST-1, or a rocky surface). Water vapor is not included in these simulations. On TRAPPIST-1e, the inclusion of water vapor might substantially increase the temperature of the cold points through transport of latent heat from substellar to anti-substellar regions, and through the greenhouse effect of water vapor. On colder planets, the effect should be marginal.

5.2.1 Glacial flow

There are in fact two distinct processes that could limit the growth of CO₂ ice glaciers:

• 1.

  The gravity pushes the glaciers to flow from the nightside to the dayside where CO₂ ice can be sublimated. This limit depends mostly on the gravity of the planet and the rheological properties of CO₂ ice (e.g., viscosity).

• 2.

  The internal heat flux of the planet causes the basal melting of the CO₂ ice glaciers. In such conditions, glaciers would slip and flow to the dayside where they would sublimate. This limit depends mostly on the geothermal heat flux of the planet and the thermodynamical properties of CO₂ ice (e.g., thermal conductivity).

It has, in fact been shown in similar conditions (Turbet et al. 2017b), owing to the low conductivity of CO₂ ice (${\lambda }_{{\text{CO}}_2}$ ~ 0.5 W m⁻¹ K⁻¹; Schmitt et al. 1997 – Part I, thermal conductivity of ices, Fig. 4), that it is mostly the basal melting that controls the maximum size of a CO₂ ice glacier. Using nightside temperatures from GCM simulations (roughly indicated in Fig. 7), we solve a set of two equations similar to Eq. (2) (the only difference being the effect of the partial pressure of the background gas, N₂ here) to derive the maximum thickness of nightside CO₂ ice deposits. These equations are: 1) the solid-liquid thermodynamical equilibrium at the base of the glacier; and 2) the linear relationship between top and bottom glacier temperatures, assuming a fixed lapse rate (conductive regime) forced by the geothermal heat flux.

For 1 bar of background gas (N₂), TRAPPIST-1e and f should be protected from CO₂ atmospheric collapse. However, for TRAPPIST-1g (h), CO₂ could collapse, and as much as 900 / 200 / 80 m of CO₂ (1000 / 250 / 100 m) could be trapped on the nightside for geothermal heat fluxes of 50 / 200 / 500 mW m⁻². This corresponds roughly to Global Equivalent Pressure of 45 / 10 / 4 bars (50 / 12 / 5 bars) of CO₂ that could be trapped.

We note that these quantities are of the same order of magnitude than the amount of CO₂ outgassed in the Venusian atmosphere (~90 bars), or the amount of CO₂ contained in the Earth’s surface, mostly in the form of carbonate rocks on the continents (~ 10² bars; Walker 1985).

Fig. 8 Maps of CO2 ice condensation(+)/sublimation(–) mean day rates (averaged over 50 Earth days) for TRAPPIST-1g and h. Wind vectors at 5 km are presented as black arrows (see the 50 m s−1 arrow for the normalization). The black line contour indicates the horizontal extent of the CO2 ice clouds (at the 1 g m−2 level). It is assumed here that the planets are endowed with a pure CO2 atmosphere and have with a surface that is entirely covered with CO2 ice. CO2 ice albedo is arbitrarily fixed at 0.5. We remind the reader that the radiative effect of CO2 ice clouds is not included here.

5.2.2 When CO₂ ice caps are full

When the CO₂ nightside ice cap becomes "full" (e.g., when CO₂ ice starts to convert into liquid), all the extra CO₂ ice (or liquid) that reaches the irradiated side sublimates (or vaporizes) into the atmosphere. The extra (now gaseous) CO₂ increases the greenhouse effect of the atmosphere. It tends to warm the nightside and thus strengthens the CO₂ ice glacial flow, leading to even more CO₂ ice sublimation. Depending on the level of irradiation, the planet either finds an equilibrium (with stable CO₂ ice and/or liquid deposits) or enters into a CO₂ runaway greenhouse. This scenario has previously been explored for H₂O-covered planets (Leconte et al. 2013b) and is extended here to the case of CO₂.

It has previously been shown (Turbet et al. 2017b; Appendix C) that when full, CO₂ ice caps should be unstable on planets:

• 1.

  that have a low enough geothermal heat flux (typically lower than ~1 W m⁻² for the present study)

• 2.

  and thatabsorb irradiation fluxes equal or larger than TRAPPIST-1f.

CO₂ ice caps are expected to be entirely injected in the atmosphere.

To test this idea, we performed 3D GCM simulations of the four TRAPPIST-1 outer planets (efgh) endowed with a pure CO₂ atmosphere where we artificially entirely covered the surface of the planets with CO₂ ice (with a CO₂ ice cover that is large enough that CO₂ ice is always present everywhere on the surface). CO₂ ice albedo is arbitrarily fixed at 0.5. For TRAPPIST-1e and f, we find that no equilibrium is possible. The planets cannot maintain surface CO₂ ice on their dayside and should always end up in a CO₂ runaway greenhouse. No equilibrium is possible until 1) all the CO₂ ice and/or liquid content has been sublimed or vaporized, or 2) the CO₂ gas content becomes so large that CO₂ greenhouse effect starts to saturate (whereas the CO₂ condensation temperature increases), as discussed in von Paris et al. (2013a).

The scenario is however different for planets that are substantially less irradiated: TRAPPIST-1g and h. For these planets, 3D GCM simulations (see Fig. 8) indicate that an equilibrium where CO₂ ice and gaseous coexist is possible. In this configuration, the surface temperatures are roughly constant over the planet. For a pure CO₂ atmosphere, we find an equilibrium at 150 ± 1 mbar and 174 ± 0.1 K for TRAPPIST-1g (4 ± 0.1 millibars and 145 ± 0.1 K for TRAPPIST-1h). The dayside intense CO₂ ice sublimation is offset by the nightside condensation, as illustrated in Fig. 8. For TRAPPIST-1g, approximately 6 m of CO₂ per Earth year (4 m per Earth year for TRAPPIST-1h) is expected to get sublimed near the substellar point.

This tells us that if TRAPPIST-1g somehow progressively accumulates enough CO₂ on its nightside so that it starts to spill on its dayside and get sublimed, the planet should not be able to accumulate enough CO₂ in the atmosphere to reach the warm state depicted in Fig. 7. Instead, the planet should be permanently trapped in a cold state, with CO₂ ice covering potentially as much as the entire surface of the planet.

Conversely, we note that – for TRAPPIST-1g only – if the planet initially starts with a large content of gaseous CO₂ (so that it lies above the unstable lines, in Fig. 7), then CO₂ ice and/or liquid deposits are unstable.

5.2.3 CO₂ ice gravitational stability

We have assumed here that some CO₂ has condensed on the nightside of TRAPPIST-1 outer planets, above the water ice shell. CO₂ ice is 1.6 times denser than water ice. This difference of density between the two types of ice (CO₂ above, H₂O below) should trigger an instability of Rayleigh-Taylor (Turbet et al. 2017b) that forces CO₂ ice to sink below the water ice cover. At first order, and assuming that both layers of CO₂ and H₂O ices are isoviscous (i.e., have a fixed viscosity), the density contrast should initiate Rayleigh-Taylor instabilities at a timescale τ_R-T given by (Turcotte & Schubert 2001): $${\tau }_{\text{R-T}}=\frac{13 \eta }{\Delta \rho gb},$$(3)

with η the viscosity of the more viscous layer, Δρ the density contrast between the two layers and b the characteristic size of the domain. We discuss in the following paragraph how we estimated the different terms of this equation.

Depending on whether CO₂ is liquid or solid at the interface with the H₂O layer, the density contrast Δρ would range between 240 and 570 kg m⁻³. Depending on the planet and the background gas content, the surface temperature at the cold point is expected to range between 120 K and 160 K (see Fig. 7). However, the basal temperature should rapidly increase with the glacier thickness, given the low conductivity of CO₂ (Schmitt et al. 1997). Assuming that the CO₂ glaciers are nearly full, they should have a thickness ~ 10² m and the basal temperature could be as high as ~218 K (temperature at the CO₂ liquid-solid equilibrium for a pressure of 1.5 MPa). For a stress at the interface between the two layers on the order of 1 MPa and a temperature ~ 218 K, the viscosity of the CO₂ ice layer is estimated ~ 10¹² Pa s, based on available experimental data (Durham et al. 1999). At the same temperature and stress conditions, water ice has a viscosity ≤ 10¹⁶ Pa s, for grain size lower than 1 mm, based on experimental data (Durham et al. 2001; Durham & Stern 2001; Goldsby & Kohlstedt 2001). The water ice layer should thus be the layer controlling the Rayleigh-Taylor timescale. Assuming that a thickness of 10² m of CO₂ deposit is representative of the characteristics size of the domain, the R-T timescale τ_R-T is ≤ 10⁴ Earth years, whichis geologically short.

Once gravitationally destabilized, the CO₂ ice deposit would sink at the base of the water ice shell at a rate that is determined mostly by the viscosity of water ice and the size of the CO₂ ice diapir (e.g., the domed CO₂ ice formation piercing the overlying water ice shell). The time required for a CO₂ ice diapir to cross the water ice layer can be estimated using the Stokes velocity, the terminal velocity of a sphere falling through a constant viscosity medium (Ziethe & Spohn 2007): $$U_s = \frac{2}{9} \Delta \rho g (r^2/\eta ).$$(4)

For a diapir radius r of 100 m (comparable to the thickness of the CO₂ deposit) and a conservative value for water ice viscosity of 10¹⁵ –10¹⁶ Pa s, this leads to a velocity of 0.04–0.4 m per Earth year. As temperature increases as a function of depth (~ 2 F_geo K m⁻¹), the viscosityof water ice is expected to decrease with depth, resulting in an acceleration of the diapir fall. A 100-m diapir of CO₂ ice would thus not need more than ~10⁴ Earth years to reach the bottom of a 1.5-km thick water ice layer, which is the expected depth of a subglacial ocean for a geothermal heat flux ~ 0.1 W m⁻².

These two calculations (Rayleigh-Taylor and diapir fall timescales) tell us that the lifetime of surface CO₂ ice on TRAPPIST-1 planets should be geologically short. In particular, it should be short compared firstly to the volcanic CO₂ outgassing timescale. The present-day Earth CO₂ volcanic outgassing rate is 60 bars/Gy (Brantley & Koepenick 1995; Jarrard 2003). It takes roughly 10⁶ Earth years to outgass a ~60 cm Global Equivalent Layer (GEL) of CO₂ ice (equivalent to a 10² m-thick nightside CO₂ ice cap with a radius of 10³ km).

Secondly, the lifetime of surface CO₂ ice should be short compared to the CO₂ ice flow and sublimation timescale. We have assumed that a 10² m-thick, 10³ km-radius nightside CO₂ ice cap is in dynamic equilibrium with the atmosphere. This means that the CO₂ ice flow – controlled here by the rheological properties of CO₂ ice – has reached a constant, positive rate. This also means that the integrated CO₂ ice sublimation rate at the edges of the glacier is equal to the total gaseous CO₂ condensation rate on the ice cap. We modeled the steady state flow of the CO₂ ice cap using Eqs. (1)–(4) from Menou (2013). We take the flow rate constants of CO₂ ice from Nye et al. (2000), derived from the measurements of Durham et al. (1999). We chose the rheological properties of CO₂ ice for a creep exponent n = 2, and at 218 K. This is the maximum temperature expected at the bottom of the CO₂ ice glacier, before basal melting occurs. These are conservative assumptions in the sense that these are the parameters (creep law and temperature) that maximize the velocity of the CO₂ ice flow. With these assumptions, we estimate that it takes at least 10⁸ Earth years to recycle the entire CO₂ ice cap.

In summary, CO₂ that starts to collapse on the nightside of TRAPPIST-1 planets is expected to form permanent CO₂ ice caps. This condensed CO₂ can be entirely returned into the atmosphere if:

• 1.

  the CO₂ ice caps reach their maximum size, defined by the condition of basal melting; This can be achieved through CO₂ volcanic outgassing;

• 2.

  the CO₂ ice caps are removed efficiently from the nightside by the dry, rheological CO₂ ice flow, and then get sublimated.

However, these two processes operate on much longer timescales than the time required to burry the CO₂ ice caps below the water ice cover. Consequently, CO₂ that condense on the nightside of TRAPPIST-1 planets (and by extent, any water-rich, synchronous planet) can very likely end up buried below the water ice shell.

5.2.4 The fate of buried CO₂

CO₂ ice is expected to completely melt and equilibrate thermally with the surrounding H₂O media when stabilized at the bottom of the water ice shell. The temperature and pressure conditions at the bottom of the water ice layer depend on its thickness and on the geothermal flow. For geothermal heat flux lower than ~ 0.4 W m⁻², the meltingof water ice would be reached for depth larger than ~ 4 × 10² m, and pressure of ~ 3.5 MPa, corresponding to the saturation vapor pressure of CO₂ at ~ 273 K (Lide 2004). Destabilizing the liquid CO₂ would therefore require a geothermal heat flux higher than 0.4 W m⁻². At such large geothermal heat flux, CO₂ ice (or liquid) should indeed get sublimed (or vaporized) within the water ice shell.

However, for geothermal heat flux lower than 0.4 W m⁻², CO₂ ices and/or liquids should be stable during their fall. Even if the density of liquid CO₂ decreases with increasing temperature as it equilibrates with the surrounding water ice media, it remains always denser than water ice (Span & Wagner 1996), and therefore should always accumulate at the bottom of the ice shell. At T = 273 K and pressure between 3.5 and 28 MPa (subglacial pressures estimated for geothermal heat flux between 400 and 50 mW ⁻²), liquid CO₂ has a density very close to that of liquid water (928 and 1048 kg m⁻³, respectively, using the equation of state of Span & Wagner 1996), so that CO₂ should coexist with H₂O at the ice-water interface.

From this point, two processes are expected to occur and compete with each other. Firstly, part of the CO₂ should dissolve in the liquid water. The total amount of CO₂ that could be dissolved in the water layer would depend on the volume (thickness) of the water layer.

Secondly, pressure and temperature conditions expected at the bottom of the ice layer are in the stability field of CO₂ clathrate hydrate (Sloan 1998; Longhi 2005), therefore CO₂ should rapidly interact with H₂O molecules to form clathrate hydrate.

Clathrate hydrates are nonstoichiometric compounds consisting of hydrogen-bonded H₂O molecules forming cage-like structures in which guest gas molecules, such as CO₂, can be trapped (Sloan 1998). Once formed, these clathrates are very stable and can be dissociated only if the temperature is raised about 5–10 K above the melting point of water ice. The storage of CO₂ in the formof clathrate should be particularly efficient as liquid CO₂ and liquid water coexist. As CO₂ clathrate hydrates have a density of about 1150 kg m⁻³ (assuming full cage occupancy, Sloan 1998), they would rapidly sink at the bottom of the water liquid layer, ensuring an almost complete clathration of CO₂. We note that we expect some of the CO₂ to get dissolved in the liquid water during the clathrate sinking. The relative proportion of CO₂ trapped in the form of clathrate hydrate or dissolved in the water layer would depend on the volume of CO₂ that is buried at the base of the ice shell and on the volume (thickness) of the water layer.

In summary, as long as the geothermal heat flux is lower than ~ 0.4 W m⁻², the water ice shell should exceed several hundreds of meters, and CO₂ should remain sequestered either in the form of CO₂ clathrate hydrates or dissolved in the subglacial liquid water ocean. Release of gaseous CO₂ in the atmosphere may occur in particular following local increase of geothermal heat flux resulting in a significant thinning and breaking-up of the water ice shell. The total amount of CO₂ that can be stored in the H₂O layer (by any of the two processes discussed above) depends on the total abundance of H₂O of the planet as well as the CO₂/H₂O ratio. Evaluating the maximum amount of CO₂ that can betrapped underneath the water ice cover require a detailed description of the H₂O layer structure as well as thermodynamic models predicting the partitionning of CO₂ between the different phases.

Fig. 9 Maps of surface temperatures (averaged over 50 Earth days) for TRAPPIST-1e, f, g and h, assuming planets endowed with a 1 bar N2-dominated atmosphere composed of 10% of CH4. Wind vectors at 5 km are presented as black arrows (see the 50 m s−1 arrow for the normalization). Surface albedo is arbitrarily fixed at 0.2.

6 CH₄-dominated worlds

The only solar system terrestrial-size object that possesses a thick atmosphere that deviates from the one discussed as far is Titan. Titan (0.012 S_⊕, 0.4 R_⊕) has a 1.5 bars thick N₂ -dominated atmosphere, with as much as 5 % of methane near the surface (Niemann et al. 2005). We explore in this section the possibility that TRAPPIST-1 outer planets could be hydrocarbon-rich worlds, and the possible implications.

6.1 Warm Titans

What would happen if we suddenly placed Titan at the location of each of the seven TRAPPIST-1 planets and how would each planet evolve? At the equilibrium temperatures of the four TRAPPIST-1 outer planets (e to h), the saturation vapor pressure of CH₄ ranges between 10 and 100 bars (between 5 × 10⁻² and 5 bars for C₂ H₆). Unlike Titan, we should thus expect all the methane and/or ethane content to be vaporized in the atmosphere.

To check this, we performed haze-free 3D numerical climate simulations of N₂ and CH₄ atmospheres, for various CH₄ contents, and for the four TRAPPIST-1 outer planets. Some of these simulations (for a 1 bar N₂ -dominated atmosphere with 0.1 bar of CH₄, similar to Titan) are presented in Fig. 9. Figure 10 shows the mean, maximum and minimum surface temperatures obtained for each of the TRAPPIST-1e, f, g, and h outer planets. The calculated surface temperature of the planets results from a subtle balance between the radiative cooling of stratospheric CH₄ and the greenhouse effect of tropospheric CH₄. Around a star like TRAPPIST-1, absorption of stellar radiation by CH₄ is particularly efficient around 0.89, 1.15, 1.35, 1.65, 2.3, and marginally 3.3 μm bands. Consequently, CH₄ absorption warms the upper atmosphere and also reduces the short wave irradiation flux that reaches the surface and troposphere, contributing to a cooling of the planetary surface. For example, approximately 40% of the incoming stellar radiation is able to reach the surface in the simulations shown in Fig. 9.

Despite the anti-greenhouse effect of CH₄ that tends to cool the surface temperature, it would be extremely hard for TRAPPIST-1 planets to sustain surface liquid and/or icy CH₄. This is illustrated in Fig 10 with the comparison between the saturation pressure curve of CH₄ and the calculated minimum surface temperatures (on the nightside). Additionally, a partial pressure of at least ~ 10⁻² bar of ethane (C₂ H₆) would be required for the coldest planets to start forming nightside surface lakes or seas of ethane, similar to the ones observed on Titan.

Extreme-UV flux should lead to the formation of photochemical hazes in CH₄-rich planetary atmospheres, in the same way as on Titan. Such hazes could potentially have an additional powerful anti-greenhouse (i.e., radiative cooling) effect, by absorbing and reflecting a significant part of the incoming stellar flux (Lorenz et al. 1997). TRAPPIST-1 outer planets (e, f, g, and h) should receive a EUV flux ranging between 600 and 3000 times Titan’s flux (Wheatley et al. 2017). Photochemical hazes could thus form efficiently and accumulate in the atmosphere. This could potentially cause a catastrophic cooling of the planetary surface (McKay et al. 1991), and change the aforementioned conclusions on the likeliness for TRAPPIST-1 outer planets to sustain surface condensed methane.

Nonetheless, using a 1D photochemical-climate model taking properly into account the microphysics and the radiative effect of photochemical hazes (Arney 2016), it has been shown that the thickness (and thus the opacity) of organic hazes should be self-regulated (Arney et al. 2016). In fact, thick hazes should inhibit methane photolysis, which would in turn drastically limit haze production rates. In other words, the rate of methane photolysis should not scale linearly with the incoming EUV flux, but instead should at some point saturate. Moreover, organic hazes are much less opaque at the emission wavelengths of cool stars like TRAPPIST-1 than solar emission ones (Khare et al. 1984; Vinatier et al. 2012; Arney et al. 2017). This indicates that a large part of the incoming stellar flux would reach TRAPPIST-1 planetary surfaces and tropospheres and easily vaporize all the methane in the atmosphere. In other words, even when takinginto account the radiative effect of hazes, all the TRAPPIST-1 planets should be well beyond the CH₄ runaway-greenhouse-like limit. Eventually, it is important to note that, at the high EUV fluxes expected on TRAPPIST-1 planets, CO₂ (if present) could be also photodissociated into oxygen radicals that should seriously limit the build up of the organic hazes (Arney et al. 2017). In particular, if the atmospheric CO₂/CH₄ ratio is high, and if the emission of TRAPPIST-1 is high in the spectral region ~ [120–180] nm, where the UV cross section of CO₂/O₂ is maximum(Arney et al. 2017, Fig. 2c), then the formation of photochemical hazes could be severely halted.

We remind the reader that a potentially thick O₂ atmosphere could have built up abiotically during the early runaway phase (Luger & Barnes 2015; Bolmont et al. 2017) while TRAPPIST-1 was a pre-main-sequence star, playing potentially a strong role here on the haze formation. But above all, the combustion of CH₄ (and more generally, of any reduced compound such as NH₃ or H₂S), following CH₄ + 2O₂ → CO₂ + 2H₂O, should prevent CH₄ from substantially building up in a thick O₂-rich atmosphere. If the build-up of O₂ during the early runaway phase exceeds the total reservoir of CH₄, there might not be enough room left for CH₄ to accumulate in the atmosphere.

Fig. 10 Mean surface temperatures of TRAPPIST-1 outer planets, assuming atmospheres made of N2 and CH4 only. Each simulation was performed with the 3D LMD generic global climate model, for 1 bar of N2, and various CH4 partial pressures (from 10 Pa to 10 bars). The four regions filled in colors show the range of surface temperatures reached by the four outer planets (red, yellow, green and blue for TRAPPIST-1e, f, g, and h, respectively). Solid, dashed, and dash-dotted lines depict the mean, maximum and minimum surface temperatures reached in the simulations, respectively. We note that the highest surface temperature for TRAPPIST-1f, g and h is almost always lower than 273 K (e.g., the melting point of water). The black line indicates the CH4 equilibrium vapor pressure. Dashed red lines show the equilibrium vapor temperature of ethane (C2 H6) for various partial pressures (10−2, 10−4; and 1 bar).The melting point of water is indicated by the black horizontal dashed line. As a reminder, equilibrium temperatures of TRAPPIST-1e, f, g, and h planets are (for a Titan-like bond albedo of 0.3) respectively 230, 198, 182, and 157 K.

6.2 Titan-like world lifetime

Through 1) CH₄ and hydrocarbons photodissociation, 2) organic hazes formation and 3) haze sedimentation, the atmospheric CH₄ and hydrocarbon content of TRAPPIST-1 planets should deplete rapidly. It is, for example, estimated that it should take roughly 10 My for Titan to remove all the methane (0.07 bar) from the atmosphere (Yung et al. 1984), and as much as ~ 30 bars could have been destroyed since the beginning of the solar system.

Therefore, as much as 600–3000 times more methane (averaged over the surface) could potentially be photolyzed on TRAPPIST-1 outer planets. Over the expected age of the TRAPPIST-1 system (between 5 and 10 Gyr, according to Luger et al. 2017 and Burgasser & Mamajek 2017), at least ~ 120 bars of CH₄ (Titan’s limit, including the gravity correction) and as much as 10⁵ bars of CH₄ (when scaling linearly the CH₄ loss with the EUV flux) could have been destroyed by photolysis.

Sustaining continuously a CH₄-rich (and NH₃-rich, by analogy) atmosphere over TRAPPIST-1 lifetime would require an extremely large source of methane. It is in fact widely believed that the CH₄ current level on Titan might be somewhat anormal and produced by an episodic replenishment due to destabilization of methane clathrates in Titan’s subsurface (Tobie et al. 2006).

Similarly, large quantities of N₂ could be photodissociated, forming HCN (Liang et al. 2007; Tian et al. 2011; Krasnopolsky 2009, 2014), and could be lost subsequently in longer carbonated chains that could sedimentate on the surface (see next subsection). This mechanism could remove efficiently N₂ from the atmosphere in the long term. We do, however, acknowledge that the arguments stated in the previous section (especially the haze negative feedback on CH₄ photolysis, as proposed by Arney et al. 2016) could drastically limit the CH₄ photolysis rate and relax the constraint on the methane production rate required to sustain a CH₄-rich atmosphere.

6.3 Surface conditions

Even for large CH₄ contents, and even when neglecting the radiative effect of photochemical hazes, TRAPPIST-1f, g and h should be cold enough (see Fig. 10 and the associated legend) to be covered by a complete layer of water ice. In this case, photolysis of methane would produce organic hazes that should sedimentate and progressively accumulate at the surface in large quantities. On Titan, it is estimated that ~1 m global equivalent layer (GEL) of heavy hydrocarbons – or tholins – are covering the surface (Lorenz et al. 2008). This is in fact two orders of magnitude lower than what we would expect from the direct conversion of current CH₄ photolysis rate through the age of the solar system (Lorenz & Lunine 1996). Possible solutions to this discrepency are discussed in Lorenz et al. (2008).

Similarly, signatures of long carbonated chains have also been detected on many Kuiper Belt Objects (Johnson et al. 2015; see Table 1 and references therein), including Pluto, Triton, Makemake, Sedna, etc. The New Horizon mission has even directly observed and mapped (during its flyby) dark tholins deposits on Pluto, in Cthulhu Regio (Stern et al. 2015).

TRAPPIST-1 planets could thus be covered by a thick surface layer of tholins today. Once the CH₄ atmosphericreservoir would be empty, only condensable hydrocarbons and long-carbonated chains (and potentially gaseous N₂, leftover from NH₃) would remain. In general, because sedimented organic hazes should be rather rich in hydrogen, they should have a density on the same order of magnitude than water ice, and should be in particular more stable than CO₂ to gravitational burial.

Assuming that photochemical hazes have a limited radiative effect in the near-infrared (Khare et al. 1984; Vinatier et al. 2012; Arney et al. 2017) where the emission of TRAPPIST-1 peaks, large quantities of CO₂ (added to CH₄, N₂, and other greenhouse gases) could be sufficient to raise the surface temperature of TRAPPIST-1 outer planets above the melting point of water, although this needs to be tested with coupled photochemical, 3D climate model in the future. In this case, sedimented organic carbonated chains should not accumulate at the surface but instead should get dissolved in the liquid water ocean. This, and the UV shield provided by the photochemical hazes (Wolf & Toon 2010; Arney et al. 2017) could provide TRAPPIST-1 planets with surface conditions favorable for life – as we know it – to emerge and develop. We discuss more generally in the next section the conditions required for TRAPPIST-1 planets to sustain surface habitability.

7 The habitability of TRAPPIST-1 planets

Most of our knowledge of habitability comes from the study of Venus, Mars, and Earth. The system of TRAPPIST-1 displays a fantastic zoology of planets to confront our theories with, and potentially revolutionize all what we know on this domain.

7.1 Habitability of the inner planets TRAPPIST-1bcd

The two inner planets of the system (TRAPPIST-1b and c) are likely to be too hot to sustain global oceans of liquid water (Kopparapu et al. 2013, 2016; Yang et al. 2013). Nonetheless, they could still be desert worlds with limited surface water (Abe et al. 2011)trapped in nightside niches (e.g., land planets) or at the edge of large scale glaciers near the terminator (Leconte et al. 2013b).

TRAPPIST-1d (S_eff ~ 1.14 S_⊕) however is near the inner edge of the habitable zone of synchronously-rotating planets (Yang et al. 2013; Kopparapu et al. 2016). If TRAPPIST-1d is able somehow to sustain a thick, highly reflective water cloud cover near the substellar region, it could sustain surface liquid water global oceans. Detailed 3D modeling of clouds, and more generally of all the possible parameters that could affect the atmospheric circulation, would be required to assess this possibility.

7.2 The remarkable potential of TRAPPIST-1e for habitability

According to our simulations, TRAPPIST-1e is the only planet in its system with the ability to host surface liquid water without the need of greenhouse warming from another gas than H₂O. This requires a sufficient H₂O reservoir covering the whole surface (i.e., that cannot be fully trapped on the nightside). Thanks to the synchronous rotation, the received stellar flux (F ~ 904 W m⁻²) is sufficient to maintain a least a patch of liquid water at the substellar point, even in the absence of a background atmosphere. This configuration is usually known as the eyeball regime (Pierrehumbert 2011).

This situation is similar to that of Proxima Cen b (Anglada-Escudé et al. 2016). This potentially rocky (most probable mass of 1.4M_⊕) planet orbiting the closest star from our Sun receives within uncertainties nearly the same amount of stellar energy (F ~ 890 W m⁻²) as TRAPPIST-1e. Two studies with two different GCMs (Turbet et al. 2016; Boutle et al. 2017) showed that a water-rich and synchronousscenario for Proxima b generates a substellar surface ocean. The reader is referred to Turbet et al. (2016; in particular, their Fig. 1) for a detailed discussion on the possible climate regimes on TRAPPIST-1e, analogous to Proxima Cen b.

We performed several 3D GCM simulations (see Fig. 11), assuming a cold start (T = 210 K everywhere, full water ice coverage). We find that for any atmosphere, TRAPPIST-1e always ends up with surface liquid water, at least in the substellar region. This would hold even with no background atmosphere at all; in this case, the atmosphere would be composed of water vapor. Starting from this point, adding greenhouse gases to the atmosphere would increase the mean surface temperature and increase the size of the patch of liquid water. We not only confirm here that the case of TRAPPIST-1e is analogous to Proxima Cen b, but we also show that, due to the lowered albedo of water ice around TRAPPIST-1, these conclusions do not depend on the initial state.

In summary, if 1) TRAPPIST-1e is in synchronous rotation and 2) is water-rich, then the planet should have a patch of liquid water at its substellar point, whatever its atmosphere (as thin or thick as wanted) and whatever its initial state (fully glaciated or not). This result must hold for any arbitrary atmospheric composition, unless a tremendous anti-greenhouse effect occurs (e.g., absorption by stratospheric methane, or absorption and reflection by photochemical hazes), or unless a tremendous greenhouse effect (by a thick Venus-like atmosphere, for example) raises the mean surface temperatureabove the critical point of water (647 K). These possibilities will be explored in details in future studies.

If low density estimates of TRAPPIST-1e (Gillon et al. 2017; Wang et al. 2017) were to be confirmed in the future, indicating that the planet could have retained large quantities of water, TRAPPIST-1e would thus become a fantastic candidate for habitability outside our solar system. More generally, TRAPPIST-1e together with Proxima Cen b highlight a new type of planets that are very likely to sustain surface liquid water, and that is therefore extremely promising for habitability prospects.

Fig. 11 Four-year average surface temperature maps of TRAPPIST-1e endowed with atmospheres made of N2 and 376 ppm of CO2, and for various atmospheric pressures (10 mbar, 0.1 bar, 1 bar, 4 bars, and 10 bars). Solid line contours correspond to the delimitation between surface liquid water and sea water ice. The figure in the bottom right panel indicates in blue, black and red the minimum, mean, and maximum surface temperatures, respectively. We note that the planets were assumed to be initially cold (T = 210 K everywhere) and completely covered by water ice.

7.3 The habitability of outer planets

Besides TRAPPIST-1e, the three outer planets of the system (TRAPPIST-1f, g and h) are interesting probes to study habitability outside our solar system.

7.3.1 Surface habitability cannot be sustained with background gases only

From the set of simulations described in Sect. 4.1 and extended to atmospheric pressures as thick as 4 bars, we find that none of the three TRAPPIST-1 outer planets (f, g, and h) are able to maintain surface liquid water, assuming a background atmosphere – not able to generate a significant greenhouse effect – that would only be made of N₂, CO, or O₂ (with H₂O included as a variable species). This tells us that TRAPPIST-1f, g, and h need to build up greenhouse gases in their atmosphere (e.g., CO₂, CH₄, NH₃, H₂, etc.) to sustain surface habitability. We explore this possibility in the next sections.

7.3.2 Minimum CO₂ content required for surface habitability

Using 3D and 1D (with our 1D cloud-free climate model (Wordsworth et al. 2010c) that uses the same physical package than the 3DLMD generic GCM described in Sect. 2) simulations of planets endowed with thick CO₂-dominated atmospheres, we find that:

• 1.

  Planet f can maintain surface liquid water for CO₂-dominated atmosphere thicker than ~1 bar. We note that a warm solution is possible for lower CO₂ atmospheric pressures, although CO₂ would condense on the surface, leading to a complete atmospheric collapse.

• 2.

  Planet g can maintain surface liquid water for CO₂-dominated atmospheres thicker than ~5 bars.

• 3.

  Planet h is not suitable for surface liquid water, whatever the thickness of the CO₂ atmosphere is, and even when maximizing the radiative effect of CO₂ ice clouds (parameterized following Forget et al. 2013). In fact, Fig. 7 also tells us that TRAPPIST-1h is unableto build up a CO₂ atmosphere, whatever the background gas content. For a 1-bar N₂ atmosphere, TRAPPIST-1h should not be able to build up more than few tens of ppm of CO₂ in the atmosphere.

Our results are roughly in agreement with the reference papers by Kopparapu et al. (2013, 2014) on habitability. We note however that the recent paper by Wolf (2017), which finds that CO₂-dominated atmospheresas thick as 30 bars cannot warm the surface of TRAPPIST-1f and g above the melting point of water, is at odd with our results, and more generally, with the related literature. We believe that the discrepancy comes potentially from the fact that Wolf (2017)underestimated in his calculations the effect of the lowered albedo of ice and snow around cool stars (mean water ice albedo of 0.21 around TRAPPIST-1) due to the shape of its reflectance spectrum, as supported by experimental data (Warren & Wiscombe 1980; Joshi & Haberle 2012). This would suppress the runaway glaciation positive feedback invoked by Wolf (2017).

The discrepency could also come from differences in the radiative treatment of CO₂-rich atmospheres. Here we have used the parameterization of CO₂ absorption lines following Wordsworth et al. (2010b), using updated line intensities and positions, and half-width at half-maximum from HITRAN-2012 (Rothman et al. 2013), sublorentzian line shape of Perrin & Hartmann (1989) and Collision-Induced Absorptions (CIA) of Gruszka & Borysow (1997) and Baranov et al. (2004). Wolf (2017) used instead CO₂ cross sections from Wolf & Toon (2013), based on HITRAN-2004, and did not include the effect of collision induced absorptions (CIA) and dimer absorptions despite their importance when modeling thick CO₂ atmospheres (Wordsworth et al. 2010a; Turbet & Tran 2017).

More generally, using our 1D and 3D GCM simulations, we find that the outer edge of the classical habitable zone around TRAPPIST-1 (using the TRAPPIST-1 synthetic spectrum, and for an atmosphere of 70 bars of CO₂) lies around 306 W m⁻² (S_eff = 0.225). This value is slightly higher than the one (302 W m⁻²; S_eff = 0.221) given in Kopparapu et al. (2013).

We explored the effect of 1) gravity, 2) rotation mode, 3) changing the stellar spectrum (from synthetic to blackbody), and 4) including the radiative effect of CO₂ ice clouds, and found that their cumulative effect on the limit of the outer edge of the habitable zone should not exceed roughly 30 W m⁻² in the context of the TRAPPIST-1 exoplanetary system.

7.3.3 The case of TRAPPIST-1h

TRAPPIST-1h is a poor candidate for surface habitability, for the following reasons:

• 1

  As explained in Sect. 5.1, TRAPPIST-1h is unable to accumulate a dense CO₂ atmosphere (because of surface condensation) that could warm the surface and favor surface habitability.

• 2

  As shown in Sect. 6, even 1) when considering an unlikely scenario where a CH₄ thick atmosphere would have been built, and even 2) when neglecting the radiative cooling of photochemical hazes, we find that (see Fig. 10) CH₄-dominated atmospheres are unable to raise the mean surface temperature of TRAPPIST-1h above ~160 K for an atmosphere made of 1 bar of N₂ and CH₄ content lower than 1 bar. This result is mostly due to the anti-greenhouse effect of stratospheric methane.

• 3

  Throughcollision induced absorptions, H₂ is an extremely powerful greenhouse gas that could potentially warm the surface of TRAPPIST-1h well above the melting point of water (Stevenson 1999; Pierrehumbert & Gaidos 2011; Luger et al. 2017). However, given the small size of the planet, and given the preliminary results of transit spectroscopy with HST (De Wit et al. 2016), the possibility of an H₂-rich atmosphere around TRAPPIST-1h seems unlikely.

8 Conclusions

In this paper we have used sophisticated numerical models (a N-body code and a Global Climate Model) to better constrain the nature of the TRAPPIST-1 planets. The main conclusions of our paper are summarized below.

8.1 Tidal dynamics constraints

We show that, given the low eccentricities derived from our N-body numerical simulations, the seven planets of the TRAPPIST-1 system are very likely in synchronous rotation today, with one side permanently facing their ultra-cool host star TRAPPIST-1, and one side in the permanent darkness.

Using the same N-body simulations, we also show that tidal heating is expected to be the dominant process of internal heating for the three inner planets of the system (TRAPPIST-1b, c and d). Tidal heating could play a significant role on TRAPPIST-1e (given that TRAPPIST-1 seems to be an old system and that radiogenic heating should have decreased), but should have a much less pronounced effect on the three outer planets (TRAPPIST-1f, g and h).

8.2 Climate diversity constraints

Assuming that the TRAPPIST-1 planets are all in synchronous rotation today, we detail below the main conclusions of our paper regarding the possible climates of TRAPPIST-1 planets.

8.2.1 Airless planets should remain airless

TRAPPIST-1 planets are exposed to X and EUV radiation and stellar wind atmospheric erosion, and could have lost their atmosphere earlier in their history. We show that planets that, at some point, completely lost their atmosphere are more likely to remain airless.

• 1.

  Planets that have a low internal heat flux (e.g., TRAPPIST-1e, f, g and h) have to accumulate very large quantities of volatiles on their nightside before a runaway greenhouse process re-forms a global atmosphere.

• 2.

  Planets that have a large internal heat flux (e.g., TRAPPIST-1b, c and d) would struggle to store and protect volatiles located on their nightside. The warmer temperature of the nightside should be responsible for the formation of a residual atmosphere that would be exposed to atmospheric escape.

However, both TTV analysis of the planets and the compact, resonant architecture of the system suggest that each of the TRAPPIST-1 planets could still be endowed with various volatiles today. Assuming that the four TRAPPIST-1 outer planets (e, f, g, and h) were able to retain various volatiles in their atmosphere, surface or subsurface, we summarize the last part of our results below.

8.2.2 Background atmospheres are stable regarding atmospheric collapse

TRAPPIST-1 planets are all highly resistant to complete atmospheric collapse of N₂ or any other background gas (CO, O₂). Around 10 millibarsof N₂ or any other background gas should suffice to avoid surface condensation on the nightside of the TRAPPIST-1 planets. This is an essential property, because background gases can prevent the more volatile species (e.g., CO₂ and NH₃) from collapsing on the nightside.

8.2.3 CO₂-dominated atmospheresare sensitive to atmospheric collapse and gravitational burial

If TRAPPIST-1e, f, and g outer planets have a CO₂-dominated atmosphere, this atmosphere must be very thick.

• 1.

  Thin CO₂ atmospheres would collapse permanently on the nightside of the planets. For example, a Mars-like atmosphere would be unstable on TRAPPIST-1e, f, g and h.

• 2.

  Thick (multi-bars) CO₂ atmospheres are found stable, thanks to an efficient greenhouse warming and heat redistribution. For example, a Venus-like atmosphere would be stable on TRAPPIST-1e, f and g. We note however that TRAPPIST-1h is beyond the CO₂ condensation limit.

If CO₂ somehow starts to condense on the nightside of TRAPPIST-1 outer planets, it would form CO₂ ice glaciers that would flow toward the substellar region. A complete CO₂ ice cover is not possible for TRAPPIST-1f and the inner planets because they receive an insolation that is greater that the runaway greenhouse threshold for CO₂. A complete CO₂ ice cover is found possible on TRAPPIST-1g and h only, although the CO₂ ice glaciers should be gravitationally unstable and get buried beneath the water ice shell (if present) in geologically short timescales. CO₂ could be permanently sequestred underneath the water ice cover, in the form of CO₂ clathrate hydrates or dissolved in a subglacial water ocean. This makes the presence of surface CO₂ ice deposits rather unlikely on water-rich, synchronous planets.

8.2.4 Sustaining continuously a CH₄-rich atmosphere is challenged by photochemical destruction

Given TRAPPIST-1 planets large EUV irradiation (at least ~ 10³ × Titan’s flux) and the large photodissociation rates that are associated, sustaining continuously a CH₄-rich (and NH₃-rich, by analogy) atmosphere over TRAPPIST-1 lifetime is difficult. Calculations of the surface temperatures of the three TRAPPIST-1 outer planets (f, g and h), under a CH₄-rich atmosphere, indicate that:

• 1.

  Their surface (even on the nightside) should be too warm to sustain oceans of methane and/or ethane.

• 2.

  Their surface should be to cold to sustain surface liquid water. This is mostly due to the anti-greenhouse effect of photochemical hazes and stratospheric methane. The planets could then more likely be covered by water ice.

Photochemical hazes when sedimenting could thus form a surface layer of tholins that would progressively thicken – over the age of the TRAPPIST-1 system – above the surface.

8.3 The habitability of the TRAPPIST-1 system

Remarkably, provided a sufficient H₂O reservoir is present, TRAPPIST-1e is very likely to sustain surface liquid water, at least in the substellar region. Thisstems from the synchronous rotation coupled to an ideal insolation, and is independent of the atmosphericbackground content (from no atmosphere at all, to a thick atmosphere of hundreds of bars). The H₂O reservoir should be large enough to avoid trapping on the nightside.

Conversely, TRAPPIST-1f, g, and h are unable to sustain surface habitability only with background gases (i.e., a rather transparent atmosphere that can ensure the transport of heat and the pressure broadening of absorption lines of greenhouse gases). Approximately 1 bar of CO₂ (~ 5 bars) would be needed to raise the surface temperature above the melting point of water on TRAPPIST-1f (g). A thick CH₄ atmosphere should be unable to sustain surface habitability on TRAPPIST-1f, g, and h.

TRAPPIST-1h is unable to sustain surface habitability with N₂, CO₂, CH₄, and so on. only. This could only be achieved with an unlikely, thick H₂ -dominated atmosphere.

Future atmospheric exploration of the TRAPPIST-1 system with the JWST and other forthcoming astronomical observatories is extremely promising. TRAPPIST-1 planets are about to become invaluable probes for comparative planetary science outside our solar system and possibly habitability. The results of our paper could serve to prepare and then interpret the future observations of the TRAPPIST-1 system and analogous. The various numerical climate simulations presented in this paper will be used in follow-up papers to provide the community with synthetic observables (transit spectra, phase curves, and secondary eclipses), that should be directly comparable with future JWST observations. Finally, we remind the reader that the results of this paper can be applied to any other cool Earth-sized planets orbiting in synchronous rotation around any cool to ultra-cool star.

Acknowledgements

M.T. thanks Tanguy Bertrand, Jan Vatant d’Ollone and Sébastien Lebonnois for fruitful discussions on Pluto and Titan. M.T. also thanks Kevin Olsen for his editorial feedback. M.T. thanks B. Charnay and B. Bézard for insightful discussions on the photochemistry of Titan. On behalf of M. Turbo-King, M.T. thanks warmly his colleagues B.R. Tang (aka Tanguy Bertrand), Z. Habeertable (aka Aymeric Spiga), M.C. Chouffe (aka Apurva Oza), B. Exquisit (aka David Dubois), and L. Keg-Beer (aka Laura Kerber), for their limitless inspiration. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 679030/WHIPLASH).

References

All Tables

All Figures

	Fig. 1 Eccentricity (top panels) and tidal heat flux (bottom panels) for the three inner planets of TRAPPIST-1 for different tidal dissipation factors (different colors): from 0.1 to ten times the Earth’s value (taken from Neron de Surgy & Laskar 1997). Stars indicate the mean values.
In the text

	Fig. 2 Same as Fig. 1 but for the four outer planets.
In the text

Fig. 3 Tidal heat flux map (in W m−2) as a function of semi-major axis (X axis) and eccentricity (Y axis) for planets with the Earth mass and radius. The dissipation efficiency is assumed to be one tenth of the Earth one to account for the dissipation in the mantle only. The tidal flux scales linearly with this parameter. This value can, however, easily change by orders of magnitude with the internal structure of each planet. This map should thus serve as a rough guide only.

In the text

Fig. 4 Equilibrium vapor pressures as a function of temperature for nine different species (where experimental data are available) that could be abundant on TRAPPIST-1 planets. Solid black lines where superimposed to indicate the equilibrium temperatures of TRAPPIST-1e,f,g, and h (assuming a surface albedo of 0.2). Stars indicate the positions of the triple points. These curves were adapted from Fray & Schmitt (2009).

In the text

Fig. 5 Maps of surface temperatures (averaged over 50 Earth days) for TRAPPIST-1f, g and h, assuming initially cold, water ice covered frozen planets, endowed with a pure N2 atmosphere (with H2O as a variable gas) at three different surface pressure (10 millibars, 100 millibars, 1 bar). Blue colored label indicates the mininum temperature reached by the coldest point of the planet throughout the entire simulation. As a reminder, for partial pressures of 10 millibars (100 millibars and 1 bar), N2 is expected to collapse at 53 K (62 K and 79 K), CO at 56 K (66 K and 83 K), and O2 at 61 K (72 K and 87 K). The geothermal heat flux is not taken into account, but the coldest temperatures found in these cases a posteriori show that it can be neglected.

In the text

Fig. 6 Maximum nightside thickness of various types of ice (N2, CO, O2, CH4 and CO2) – assumingthat it is limited by basal melting – as a function of the geothermal heat flux. It is assumed here that the entire atmosphere has collapsed at the cold points of the planet and that the surface temperature at the top of the glacier is controlled by the geothermal heat flux. As a reference, vertical dashed lines indicate the average surface tidal heat flux on TRAPPIST-1 planets derived from Table 2, for a tidal dissipation equal to that of the Earth. We also added (vertical solid gray line) the average geothermal heat flux on Earth, to give the reader a rough sense of the amplitude of the radiogenic heating on TRAPPIST-1 planets. These quantities can be numerically converted in term of global equivalent surface pressure when multiplied by a factor $\frac{\rho g}{2}$. The thermodynamical and rheological properties of the ices were taken from https://encyclopedia.airliquide.com; http://webbook.nist.gov; Roder (1978); Schmitt et al. (1997); Fray & Schmitt (2009); Trowbridge et al. (2016), and Umurhan et al. (2017). Missing rheological data were mimicked on N2.

In the text

Fig. 7 Climate regimes reached as a function of the partial pressures of N2 and CO2. For each set of (pN2, pCO2 ), it is indicated if the atmosphere is stable (red) or not (blue) to the atmospheric condensation/collapse of CO2. The black arrows indicate how planets that have an unstable atmosphere would evolve on this diagram. Temperatures (in green) correspond to the rough estimate of the temperature of the cold point, at the stable lower boundary (blue is up; red is down). Simulations were performed assuming a surface albedo of 0.2 (corresponding both to a water ice surface around TRAPPIST-1, or a rocky surface). Water vapor is not included in these simulations. On TRAPPIST-1e, the inclusion of water vapor might substantially increase the temperature of the cold points through transport of latent heat from substellar to anti-substellar regions, and through the greenhouse effect of water vapor. On colder planets, the effect should be marginal.

In the text

Fig. 8 Maps of CO2 ice condensation(+)/sublimation(–) mean day rates (averaged over 50 Earth days) for TRAPPIST-1g and h. Wind vectors at 5 km are presented as black arrows (see the 50 m s−1 arrow for the normalization). The black line contour indicates the horizontal extent of the CO2 ice clouds (at the 1 g m−2 level). It is assumed here that the planets are endowed with a pure CO2 atmosphere and have with a surface that is entirely covered with CO2 ice. CO2 ice albedo is arbitrarily fixed at 0.5. We remind the reader that the radiative effect of CO2 ice clouds is not included here.

In the text

	Fig. 9 Maps of surface temperatures (averaged over 50 Earth days) for TRAPPIST-1e, f, g and h, assuming planets endowed with a 1 bar N2-dominated atmosphere composed of 10% of CH4. Wind vectors at 5 km are presented as black arrows (see the 50 m s−1 arrow for the normalization). Surface albedo is arbitrarily fixed at 0.2.
In the text

Fig. 10 Mean surface temperatures of TRAPPIST-1 outer planets, assuming atmospheres made of N2 and CH4 only. Each simulation was performed with the 3D LMD generic global climate model, for 1 bar of N2, and various CH4 partial pressures (from 10 Pa to 10 bars). The four regions filled in colors show the range of surface temperatures reached by the four outer planets (red, yellow, green and blue for TRAPPIST-1e, f, g, and h, respectively). Solid, dashed, and dash-dotted lines depict the mean, maximum and minimum surface temperatures reached in the simulations, respectively. We note that the highest surface temperature for TRAPPIST-1f, g and h is almost always lower than 273 K (e.g., the melting point of water). The black line indicates the CH4 equilibrium vapor pressure. Dashed red lines show the equilibrium vapor temperature of ethane (C2 H6) for various partial pressures (10−2, 10−4; and 1 bar).The melting point of water is indicated by the black horizontal dashed line. As a reminder, equilibrium temperatures of TRAPPIST-1e, f, g, and h planets are (for a Titan-like bond albedo of 0.3) respectively 230, 198, 182, and 157 K.

In the text

Fig. 11 Four-year average surface temperature maps of TRAPPIST-1e endowed with atmospheres made of N2 and 376 ppm of CO2, and for various atmospheric pressures (10 mbar, 0.1 bar, 1 bar, 4 bars, and 10 bars). Solid line contours correspond to the delimitation between surface liquid water and sea water ice. The figure in the bottom right panel indicates in blue, black and red the minimum, mean, and maximum surface temperatures, respectively. We note that the planets were assumed to be initially cold (T = 210 K everywhere) and completely covered by water ice.

In the text