The Carbon-deficient Evolution of TRAPPIST-1c

It is currently unknown whether small M dwarfs like TRAPPIST-1 are hospitable to atmospheres at any orbital distance. M dwarfs can be highly magnetically active for up to 6–7 billion years, significantly longer than Sun-like stars, which diminish in activity by 1 billion years (West et al. 2006). Magnetic activity presents itself as high-energy X-ray and ultraviolet (XUV) radiation, flares, winds, and coronal mass ejections, all of which could be harmful to planetary atmospheres.

In the solar system, there is evidence that the terrestrial planets formed with primordial atmospheres of hydrogen and helium, which were mostly lost as the gas of the protoplanetary disk dissipated (Lammer et al. 2018; Young et al. 2023). Secondary atmospheres made up largely of volatiles such as carbon dioxide (CO₂) and water (H₂O) then formed by planetary outgassing as magma oceans solidified. Earth was able to maintain at least part of its initial water reservoir during this hot magma ocean phase, and the geological carbon cycle (i.e., carbon–silicate cycle) has subsequently recycled volatiles into and out of the atmosphere, regulating the planet's climate.

Conversely, Venus may be representative of planets that become desiccated during the magma ocean phase, because they receive too much stellar flux for their steam atmospheres to cool efficiently, triggering a runaway greenhouse (Hamano et al. 2013). Water is then lost to space, by its rapid photodissociation and subsequent hydrodynamic escape of hydrogen (Kasting 1988). This desiccation left Venus with the CO₂-dominated atmosphere that we see today. Mars also has a CO₂-dominated atmosphere, albeit much thinner, but had a warmer, wetter atmosphere in the past, most of which has been lost. These atmospheres are not static, even today, as modern probes detect the steady loss of atmospheric material from Earth, Venus, and Mars, mostly in the form of ion escape by solar wind (Edberg et al. 2011; Jakosky et al. 2018). In fact, escape has affected their atmospheres over gigayear timescales (Lammer et al. 2008).

The history of the solar system's terrestrial planets can inform our understanding of how the atmospheres of terrestrial exoplanets evolve, specifically close-in planets such as TRAPPIST-1c, which is the subject of this paper. For planets that orbit close enough to their host star to evaporate water into the atmosphere, a period of H₂O photodissociation and intense early hydrodynamic loss of hydrogen (which results in net water loss) is expected. Hydrodynamic escape occurs when the light constituent of an atmosphere reaches a high enough temperature to leave in an evaporative wind. During this hydrodynamic phase, it is also possible for heavy constituents to be dragged along by the hydrodynamic wind. Therefore, a net loss of heavier constituents such as CO₂ and O₂ can occur. Once the hydrodynamic phase is over, we expect the remaining atmosphere to be dominated by heavier molecules such as CO₂, as in the case of Venus.

Now, over the rest of the planet's lifetime, CO₂ can be lost slowly through other mechanisms. Jeans escape, another type of thermal escape, is a mechanism that occurs when the high-velocity, high-altitude particles of a gas in thermodynamic equilibrium have velocities that exceed the escape velocity of the planet. Jeans escape is not dominant in a CO₂-dominated atmosphere, because CO₂ is too heavy and cools efficiently through the 15 μm band (Gordiets & Kulikov 1985). Therefore, we expect that long-term escape on these types of planets is dominated by nonthermal escape, particularly in the form of stellar wind stripping as seen in the solar system.

Even as stellar phenomena strip material from planetary atmospheres, terrestrial planet atmospheres are replenished by geological processes, such as volcanic outgassing, as seen in recent geological history on Earth, Venus, and Mars. This outgassing may not last for the entire lifetime of the planet, as Unterborn et al. (2022) suggest for the TRAPPIST-1 planets. Nevertheless, we can model both the escape and outgassing of terrestrial planet atmospheres as coupled time-dependent processes that evolve an atmosphere through time.

In this paper, we use coupled, time-dependent simulations of outgassing and escape to model the terrestrial planet TRAPPIST-1c. Because TRAPPIST-1c orbits close to its host star, receiving slightly more flux than Venus, we assume here that it has lost effectively all H₂O and possesses a Venus-like CO₂-dominated atmosphere, as is common in the solar system. Recent observations have been taken for TRAPPIST-1c, using JWST's MIRI instrument to obtain secondary eclipse photometry in the 15 μm CO₂ absorption band. The data from these observations (Zieba et al. 2023) point to the absence of a thick CO₂ atmosphere, ruling out pure CO₂ atmospheres with surface pressures ≥0.1 bar. Zieba et al. (2023) find that a bare rock surface or a thin O₂-dominated atmosphere could also be consistent with their data. Lincowski et al. (2023) further explored the potential compositions of TRAPPIST-1c's present atmosphere. They acknowledge that CO₂ is likely to be a principal component of outgassing for Earth-sized or larger planets (Gaillard & Scaillet 2014), but test many different atmospheric species in their models. They confirm the results of Zieba et al. (2023) but additionally find that a maximum of 10% H₂O abundance is consistent with the observations.

Assuming CO₂ has dominated TRAPPIST-1c's atmosphere for most of its history, the constraint on CO₂ surface pressure from Zieba et al. (2023) allows us to constrain the history of its atmosphere with simulations and assess its early volatile inventory. We base our simulation framework and study upon the work of Kane et al. (2020), which found volatile-poor formation conditions for LHS 3844b. Our simulations differ from Krissansen-Totton & Fortney (2022), who made predictions for the atmospheres of TRAPPIST-1 planets. Krissansen-Totton & Fortney (2022) include a range of atmospheric species, including H₂O, CO₂, and O₂, and runs coupled simulations starting in the magma ocean phase all the way to present day. We first calculate the amount of CO₂ potentially lost during the hydrodynamic escape phase, and then run coupled time-dependent simulations of outgassing and escape of CO₂ starting after solidification of the magma ocean and atmospheric desiccation. Our simulations therefore track the amount of CO₂ lost in the long-term stellar wind stripping phase of TRAPPIST-1c's evolution. This allows us to determine whether long-term stellar wind stripping is efficient enough to chisel the atmosphere down to what we see today. If not, formation conditions or significant early escape must be responsible for the current shortage of CO₂ in TRAPPIST-1c's atmosphere.

We organize this manuscript as follows. In Section 2, we describe the components of our simulations of TRAPPIST-1c, including stellar properties, planet properties, planetary outgassing, and atmospheric escape. We detail how we calculate outgassing rates from interior thermal evolution and escape rates from stellar-wind-induced ion escape. We explain how we couple these processes in time-dependent simulations and vary key parameters to create a grid of potential atmospheric histories of TRAPPIST-1c. In Section 3, we first provide calculations of CO₂ loss during the early hydrodynamic escape phase, and then we analyze the results of our long-term simulations, identifying the parameters of those that fit the observational constraint. In Section 4, we discuss the implications of these results for the TRAPPIST-1 system as a whole. We address caveats of our models, diagnosing future developments that will help us better characterize the atmosphere of TRAPPIST-1c and its neighbors. In Section 5, we summarize our findings and conclude.

We initialize our simulations with the properties of the M8V star TRAPPIST-1 (Costa et al. 2006). According to Agol et al. (2021), TRAPPIST-1 has mass M_⋆ = 0.0898 ± 0.0023 M_⊙, radius R_⋆ = 0.1192 ± 0.0013 R_⊙, effective temperature T_eff = 2566 ± 26 K, and bolometric luminosity L_bol = 0.000553 ±0.000018 L_⊙.

We use the age of TRAPPIST-1 presented in Burgasser & Mamajek (2017). They conclude that it has an age of 7.6 ± 2.2 Gyr. The rotation period of TRAPPIST-1 is 3.30 ± 0.14 days as found photometrically by Luger et al. (2017). Age and rotation period are theoretically related to each other and to magnetic activity, which we address further in Sections 2.3 and 4.4.

In this manuscript, we focus on the second-closest planet to the host star, TRAPPIST-1c. TRAPPIST-1c has mass M_p = 1.308 ± 0.056 M_⊕, radius ${R}_{p}={1.097}_{-0.012}^{+0.014}\ {R}_{\oplus }$, and semimajor axis a = 0.01580 ± 0.00013 au (Agol et al. 2021). These planetary properties are displayed with the stellar properties in Table 1. Constraints on mass and radius provide an estimate for the density of TRAPPIST-1c, approximately ρ = 5.45 g cm⁻³, which is almost exactly Earth-like. The surface gravity of TRAPPIST-1c is 10.65 ± 0.42 m s⁻², slightly larger than Earth's. Given the semimajor axis and stellar luminosity, we calculate the incident bolometric flux on TRAPPIST-1c to be approximately F_inc = 2.22F_inc,⊕, which is near the value for Venus (1.91F_inc,⊕).

Table 1. Relevant Properties of TRAPPIST-1 and TRAPPIST-1c as Cited in the Text

Parameter	Value
Star—TRAPPIST-1
Mass M⋆	0.0898 ± 0.0023 M⊙
Radius R⋆	0.1192 ± 0.0013 R⊙
Effective temperature Teff	2566 ± 26 K
Bolometric luminosity Lbol	0.000553 ± 0.000018 L⊙
Age	7.6 ± 2.2 Gyr
Rotation period Prot	3.30 ± 0.14 days
Planet—TRAPPIST-1c
Mass Mp	1.308 ± 0.056 M⊕
Radius Rp	${1.097}_{-0.012}^{+0.014}\ {R}_{\oplus }$
Semimajor axis a	0.01580 ± 0.00013 au

Download table as:  ASCIITypeset image

The aforementioned properties of TRAPPIST-1c are compared to those of its neighbors, solar system planets, and other likely rocky exoplanets in Figure 1. The left panel displays the incident flux, escape velocity, and stellar temperature for the TRAPPIST-1 planets, Earth, Venus, and Mars. The right panel displays the same information for the TRAPPIST-1 planets and other exoplanets whose atmospheres have been observed, or will be observed in JWST Cycle 1 and Cycle 2 General Observers and Guaranteed Time Observations programs. Black Xs are placed on the exoplanets LHS 3844b, GJ 1252b, K2-141b, and TRAPPIST-1b, which have been observed to have no thick atmosphere (Kreidberg et al. 2019; Crossfield et al. 2022; Zieba et al. 2022; Greene et al. 2023). It is plausible that the existence of atmospheres has some dependence on all three variables. The closer a planet is to its host star, the more high-energy radiation and bombardment of stellar material it is expected to receive. Higher escape velocity makes it less likely that atmospheric material will escape at a given temperature and mean molecular mass. Finally, atmospheres may be favored in cases of higher stellar temperature where high-activity pre-main-sequence phases are shorter. The gray line in the right panel is the "cosmic shoreline" (${F}_{\mathrm{inc}}\propto {v}_{\mathrm{esc}}^{4}$) presented in Zahnle & Catling (2017), and is an empirical division of planets known with and without atmospheres. This dividing line might suggest that the TRAPPIST-1 planets should all have atmospheres, but that appears not to be the case. The exact effect that incident flux, escape velocity, and stellar temperature have on the presence of atmospheres is thus still a lingering question and must be explored through observations and simulations.

Zoom In Zoom Out Reset image size

In our simulations, the CO₂ atmosphere of TRAPPIST-1c evolves through two mechanisms: planetary outgassing and atmospheric escape. Outgassing is a result of processes in the planet's interior that are driven by mantle convection, thereby cooling the mantle over time. Tracking the mantle temperature over time along with several other important quantities allows outgassing rates to be calculated. Atmospheric escape is a result of many complex physical and chemical processes that we do not model directly here, but we instead employ the results of MHD models that trace ion escape driven by the stellar wind. We discuss this choice further in Section 4.

2.3.1. Planetary Outgassing

We employ a new version of the thermal evolution model presented in Foley & Smye (2018), translated to Python and available at https://github.com/katieteixeira/atmospheric_evolution and doi:10.5281/zenodo.10052104 (Teixeira & Foley 2023), to calculate outgassing rates. This thermal evolution code models Earth-like planets in the stagnant-lid regime, where lithospheres are rigid and immobile, but volcanism still occurs. As far as we know, Earth is the only solar system planet that has plate tectonics. Both Venus and Mars have stagnant lids, so we expect this tectonic mode to be common among terrestrial planets (Foley & Driscoll 2016). The likely lack of liquid water on TRAPPIST-1c also disfavors plate tectonics, as water has been shown to be important for plate tectonics (Hirth & Kohlstedt 1996). Water lowers the strength of the lithosphere, weakening plate boundaries and allowing subduction to occur (Regenauer-Lieb et al. 2001; Korenaga 2007).

The original code presented in Foley & Smye (2018) couples the thermal evolution of the interiors of stagnant-lid planets to melting, crustal growth, weathering, and CO₂ outgassing. Convection transports heat and carbonated material to the base of the stagnant lid, where melting and eruption allow them to escape to the surface. On stagnant-lid planets, this volcanism itself is responsible for crustal growth and the creation of weatherable rock. If carbon is weathered back into the crust, it can be buried and recycled either into the mantle, by carbon foundering, or into the atmosphere, by metamorphic outgassing.

As in Foley & Smye (2018), the temperature of the mantle evolves according to the following equation:

$$\begin{eqnarray}&&{V}_{\mathrm{man}}\rho {c}_{p}\displaystyle \frac{{{dT}}_{p}}{{dt}}={Q}_{\mathrm{man}}-{A}_{\mathrm{man}}{F}_{\mathrm{man}}-{f}_{m}{\rho }_{m}({c}_{p}{\rm{\Delta }}{T}_{m}+{L}_{m})\end{eqnarray} \tag{ 1 }$$

where V_man is the volume of convecting mantle beneath the stagnant lid, ρ is the bulk density of the mantle, c_p is the heat capacity, T_p is the potential temperature of the upper mantle, t is time, Q_man is the radiogenic heating rate in the mantle, A_man is the surface area of the top of the convecting mantle, F_man is the heat flux from the convecting mantle at the base of the stagnant lid, f_m is the volumetric melt production, ρ_m is the density of the mantle melt, L_m is the latent heat of fusion of the mantle, and ΔT_m is the temperature difference between the melt erupted at the surface and the surface temperature. This assumes that heat from the core does not play a role in the thermal evolution of the mantle.

Coupled to Equation (1) and several other differential equations, the following equation calculates the rate at which carbon is outgassed from the mantle to the atmosphere:

$$\begin{eqnarray}&&\displaystyle \frac{{{dR}}_{\mathrm{man}}}{{dt}}=-\displaystyle \frac{{f}_{m}{R}_{\mathrm{man}}[1-{(1-\phi )}^{1/{D}_{{\mathrm{CO}}_{2}}}]}{\phi ({V}_{\mathrm{man}}+{V}_{\mathrm{lid}})}\end{eqnarray} \tag{ 2 }$$

where R_man is the carbon reservoir of the mantle, ϕ is the melt fraction, ${D}_{{\mathrm{CO}}_{2}}$ is the distribution coefficient of CO₂, and V_lid is the volume of the stagnant lid. This assumes that a certain fraction of CO₂ enters the melt phase and a fraction of the melt degasses its CO₂ to the atmosphere. We refer the reader to Foley & Smye (2018) for more details.

By default, the thermal evolution code assumes that water is present, and thus weathering occurs. Since we assume TRAPPIST-1c has no liquid water on its surface, we "turn off" weathering in the code. No CO₂ is returned to the crust, and therefore no CO₂ is recycled into the mantle, so metamorphic outgassing and carbon foundering do not operate. This implies that all CO₂ that is outgassed will simply build up in the atmosphere. This is captured in Figure 2 with the major processes that potentially move CO₂ into and out of TRAPPIST-1c's atmosphere.

Zoom In Zoom Out Reset image size

2.3.2. Atmospheric Escape

We parameterize the atmospheric escape component of our simulations based on Dong et al. (2018), which used the BATS-R-US MHD code (van der Holst et al. 2014) to model the stellar wind of TRAPPIST-1 and its interaction with the seven TRAPPIST-1 planets. BATS-R-US has been tested and validated on solar system planets and is based on the composition of Mars and Venus (e.g., Ma et al. 2004, 2013). To obtain an upper bound on escape rates, the TRAPPIST-1 planets are assumed to have no intrinsic magnetic field that could protect their atmospheres. The stellar magnetic field induces small planetary magnetic fields, but the upper atmospheres of the planets are still subject to erosion. The model presented in Dong et al. (2018) captures photoionization, ion–neutral chemistry, and electron recombination chemistry that occurs in the upper atmosphere, above the neutral atmosphere dominated by CO₂. The model then solves the MHD equations for the four ion fluids that exist in this case: H⁺, O⁺, ${{\rm{O}}}_{2}^{+}$, and ${\mathrm{CO}}_{2}^{+}$. The physical processes that are captured by solving these equations include charge exchange in the stellar wind, ion pickup, and ion sputtering, all of which act to remove atmospheric material by accelerating it away from the planet. Therefore, Dong et al. (2018) are able to calculate the ion escape rates from the TRAPPIST-1 planets, which we add and use here as a total mass-loss rate in kg s⁻¹. Ion escape rates are calculated at the point where the stellar wind produces the maximum total pressure at the planets.

Given a stellar mass-loss rate from TRAPPIST-1 of ${\dot{M}}_{\star }\,=2.6\times {10}^{8}$ kg s⁻¹, the atmospheric loss rate from TRAPPIST-1c is calculated to be ${\dot{M}}_{\mathrm{atm}}=57\,\mathrm{kg}$ s⁻¹. Dong et al. (2018) find that the atmospheric loss rate numerically calculated in their simulations scales with planet radius, semimajor axis, and stellar mass-loss rate as follows:

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{atm}}\propto {\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{2}{\dot{M}}_{\star }.\end{eqnarray} \tag{ 3 }$$

This scaling relation allows us to test different stellar mass-loss rates from the literature. For example, Garraffo et al. (2017) cite the stellar mass-loss rate of TRAPPIST-1 as ${\dot{M}}_{\star }=1.89\,\times {10}^{9}$ kg s⁻¹, almost an order of magnitude higher than the value used by Dong et al. (2018). This produces an atmospheric loss rate of ${\dot{M}}_{\mathrm{atm}}=421\,\mathrm{kg}$ s⁻¹.

These instantaneous mass-loss rates do not take into account the change in the stellar wind over time. We estimate this by connecting several empirical relationships for M dwarfs. As previously described, the age, rotation, high-energy radiation, and mass-loss rates of stars are thought to be intimately connected. In fact, theoretical and observational work has been done to model the stellar mass-loss rate of M dwarfs as a function of age based on their X-ray luminosities and rotation periods. We employ scaling relations from Engle & Guinan (2018) and Magaudda et al. (2020) to calculate rotation period and X-ray luminosity, respectively:

$$\begin{eqnarray}&&{P}_{\mathrm{rot}}[\mathrm{days}]=\displaystyle \frac{\mathrm{age}[\mathrm{Gyr}]-0.012\pm 0.221}{0.061\pm 0.002}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{L}_{{\rm{X}}[\mathrm{erg}\,{{\rm{s}}}^{-1}]}=\left\{\begin{array}{ll}{C}_{\mathrm{sat}}{P}_{\mathrm{rot}}^{{\beta }_{\mathrm{sat}}}, & {P}_{\mathrm{rot}}\leqslant {P}_{\mathrm{rot},\mathrm{sat}}\\ {C}_{\mathrm{unsat}}{P}_{\mathrm{rot}}^{{\beta }_{\mathrm{unsat}}}, & {P}_{\mathrm{rot}}\gt {P}_{\mathrm{rot},\mathrm{sat}}\end{array}\right.\end{eqnarray} \tag{ 5 }$$

where ${C}_{n}=({L}_{{\rm{X}},n}/{P}_{\mathrm{rot}}^{{\beta }_{n}})$ with n = (sat, unsat), β_sat = −0.19 ± 0.11, β_unsat = −3.52 ± 0.02, P_rot,sat = 33.7 ± 4.5 [days], and log(L_X,sat(P_rot = 1 [day])) = 28.54 ± 0.20 [erg s⁻¹]. The stellar mass-loss rate can then be calculated from the X-ray surface flux as follows (Wood et al. 2021):

$$\begin{eqnarray}&&\displaystyle \frac{\dot{{M}_{\star }}}{\mathrm{SA}}\left[\displaystyle \frac{{\dot{M}}_{\odot }}{{\mathrm{SA}}_{\odot }}\right]={10}^{(0.77\pm 0.04\cdot {\mathrm{log}}_{10}({F}_{{\rm{X}},\mathrm{surface}})-3.42)}\end{eqnarray} \tag{ 6 }$$

where SA is the surface area of the star. Note that Equation (4) does not correctly predict the known rotation period of TRAPPIST-1 at its current age. These scaling relations only provide estimates for the evolution of these quantities for a representative sample of M dwarfs, and there is significant dispersion for individual stars.

In this work, we use three different atmospheric loss rate prescriptions: one based on the stellar mass-loss rate presented in Dong et al. (2018; Low), another based on the stellar mass-loss rate presented in Garraffo et al. (2017; High), and a variable stellar mass rate based on the scaling relations for M dwarfs (Variable). These three different cases are shown in Figure 3.

Zoom In Zoom Out Reset image size

We ran 2,430,000 total simulations, varying atmospheric loss rate prescription (Low, High, and Variable) and four key geological parameters, three of which are initial conditions for TRAPPIST-1c after the hydrodynamic escape phase of its evolution:

• 1.  
  C_tot: Initial carbon budget of the planet. This is the amount of CO₂ in moles that the planet initially possesses after solidification of the magma ocean and atmospheric desiccation, all located in the mantle unless otherwise specified. Earth's estimated budget is 10²² mol according to Foley & Smye (2018). Venus' C_tot must be at least 10²² mol, given that is approximately how much exists in the atmosphere at present.
• 2.  
  Heat-producing elements (HPEs): Initial budget of radiogenic heat-producing elements (²³⁸U, ²³⁵U, K, and Th) in the mantle. These elements act to heat the mantle and may prolong outgassing.
• 3.  
  μ_ref: Mantle reference viscosity. This measures the convecting mantle material's resistance to flow, specifically at a mantle temperature equal to that of present-day Earth.
• 4.  
  T_init: Initial mantle temperature. This is the starting temperature of the mantle before significant cooling occurs, and it depends on the planet's initial thermal budget, which may come from heat from formation, impacts, etc.

Our parameter space grid comprised 30 different values of each of the four geological parameters, evenly chosen from the following distributions:

$$\begin{eqnarray}&&{{\rm{log}}}_{10}\left(\displaystyle \frac{{C}_{{\rm{tot}}}}{{C}_{{\rm{tot}},\oplus }}\right)\sim {\mathscr{U}}\left({{\rm{log}}}_{10}\left(0.025),{{\rm{log}}}_{10}(2.5\right)\right)\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{HPE}}{{\mathrm{HPE}}_{\oplus }}\sim {\mathscr{U}}\left(0.5,2\right)\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\left(\displaystyle \frac{{\mu }_{\mathrm{ref}}}{\mathrm{Pa}\,{\rm{s}}}\right)\sim {\mathscr{U}}\left(20,22\right)\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{init}}}{{\rm{K}}}\sim {\mathscr{U}}\left(1700,2000\right).\end{eqnarray} \tag{ 10 }$$

Here C_tot,⊕ = 10²² mol according to an estimate of Earth's mantle and surface reservoirs from Sleep & Zahnle (2001). Our chosen range for C_tot represents a range of 10⁻⁴–10⁻² wt% concentration of CO₂ in the mantle. We test a range of initial budgets for radiogenic heat-producing elements that spans 50%–200% of Earth's budget, about 80–320 TW, estimated from stellar abundance measurements of radionuclides (Unterborn et al. 2015; Botelho et al. 2019). Our range of mantle reference viscosity encompasses the range of typical estimates for Earth's value. Finally, we consider a wide range of initial mantle temperature, which is relatively unconstrained for Earth.

During the hydrodynamic hydrogen loss phase, we emphasize that it is possible for CO₂ and other components to be dragged along with the hydrodynamic wind due to the high XUV flux of the star (Fleming et al. 2020). Here we estimate the length of the hydrodynamic H loss phase and the amount of CO₂ lost during it according to Odert et al. (2018). For a hydrogen-dominated atmosphere with a CO₂ component, the XUV-driven hydrodynamic flux of hydrogen in kg s⁻¹ m⁻² is given by

$$\begin{eqnarray}&&{F}_{{\rm{H}}}=\displaystyle \frac{{\beta }^{2}\eta {F}_{\mathrm{XUV}}}{4{\rm{\Delta }}{\rm{\Phi }}({m}_{{\rm{H}}}+{m}_{{\mathrm{CO}}_{2}}{f}_{{\mathrm{CO}}_{2}}{x}_{{\mathrm{CO}}_{2}})}.\end{eqnarray} \tag{ 11 }$$

Here β = R_XUV/R_p is the ratio of the effective radius out to which the bulk of the XUV radiation is absorbed to the radius of the planet, η is the efficiency of heating, F_XUV is the XUV flux, ΔΦ = GM_p /R_p , m_H is the mass of a hydrogen atom, ${m}_{{\mathrm{CO}}_{2}}$ is the mass of a CO₂ molecule, ${f}_{{\mathrm{CO}}_{2}}$ is the mixing ratio of CO₂ to H, and ${x}_{{\mathrm{CO}}_{2}}$ is the fractionation factor of CO₂. The fractionation factor ${x}_{{\mathrm{CO}}_{2}}$ is calculated from the following equation:

$$\begin{eqnarray}&&{x}_{{\mathrm{CO}}_{2}}=1-\displaystyle \frac{g({m}_{{\rm{H}}}-{m}_{{\mathrm{CO}}_{2}})b(T)}{{F}_{{\rm{H}}}{k}_{{\rm{B}}}T(1+{f}_{{\mathrm{CO}}_{2}})}\end{eqnarray} \tag{ 12 }$$

where b is the binary diffusion parameter and T is the temperature of the upper atmosphere. b = 8.4 × 10¹⁹ T^0.6 m⁻¹ s⁻¹ for H and CO₂. The flux of CO₂ hydrodynamically dragged is given by

$$\begin{eqnarray}&&{F}_{{\mathrm{CO}}_{2}}={F}_{{\rm{H}}}{f}_{{\mathrm{CO}}_{2}}{x}_{{\mathrm{CO}}_{2}}.\end{eqnarray} \tag{ 13 }$$

In order to track the amount of H and CO₂ over time, Equations (11) and (12) must first be solved for F_H. Then, this equation, along with Equation (13), must be integrated using a differential equation solver. When performing these calculations, we assumed β = 1, η = 0.15 (Salz et al. 2015; Kubyshkina et al. 2018), and F_XUV is the sum of EUV and X-ray flux:

$$\begin{eqnarray}&&{F}_{\mathrm{XUV}}={F}_{\mathrm{EUV}}+{F}_{{\rm{X}}}.\end{eqnarray} \tag{ 14 }$$

We calculate X-ray flux as in Equation (5) and EUV flux according to Sreejith et al. (2020) and Boudreaux et al. (2022). We varied T between 300 and 8000 K (Kulikov et al. 2006) and found that it had minimal effect on the outcome. Finally, we calculated the amount of hydrogen that could be derived from a certain number of Earth oceans, and varied the number of Earth oceans (EOs) between 0.01 and 100. This initial amount of water had the greatest effect on our calculations. If TRAPPIST-1c started with 0.01 EOs, its hydrodynamic loss phase could have lasted for ∼1000 yr and it could have lost about 0.1 bar of CO₂. Starting with 1 EO, the hydrodynamic phase could last for about 100,000 years, and about 10 bar of CO₂ could be lost. With 100 EOs, the hydrodynamic phase could last for ∼10 Myr, and about 1000 bar of CO₂ could be lost. We see that the hydrodynamic phase is extremely short compared to the lifetime of TRAPPIST-1c in all cases, but the amount of CO₂ lost during this phase can vary widely.

After calculating the rapid hydrodynamic CO₂ loss, we focus on the evolution of CO₂ over the rest of TRAPPIST-1c's lifetime. In our long-term, coupled simulations, we track several different quantities that vary as a function of time, including mantle temperature, total amount of CO₂ in the mantle, thickness of the lithosphere, and heat production from the radiogenic elements in the crust and mantle. Importantly, the observable quantity that we track is the amount of CO₂ present in the atmosphere, which can be represented by mass in kg or surface pressure (pCO₂) in bar.

The modeled evolution of pCO₂ on TRAPPIST-1c is shown for the edge cases of our varied parameters in Figure 4. The baseline parameters are C_tot[⊕] = 0.25, μ_ref = 1.0 × 10²¹ Pa s, T_init = 2000 K, and HPE[⊕] = 1 in each panel. All simulations are characterized by a sharp increase in pCO₂ in the first 1 Gyr of evolution due to early vigorous outgassing, as shown in Figure 5. Some simulations result in pCO₂ that is relatively constant over the rest of TRAPPIST-1c's lifetime; others result in order-of-magnitude losses in pCO₂. Specifically, the Low atmospheric loss rate prescription results in a minimal decrease in pCO₂ in most cases. The High prescription leads to significant decreases in pCO₂ that deplete the atmosphere by TRAPPIST-1's age in most cases. The Variable prescription produces pCO₂ curves that decrease slightly in the 1-6 Gyr range but then level out as stellar activity decreases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The four geological parameters each have some noticeable effect on the evolution of pCO₂. Most clearly, the initial carbon budget C_tot has a significant effect on the value of pCO₂ that results from early outgassing. For example, C_tot = 2.5C_tot,⊕ leads to pCO₂ almost as high as 200 bar, but with C_tot = 0.025C_tot,⊕, pCO₂ does not even reach 2 bar. We see that in the C_tot = 2.5C_tot,⊕ case, none of the atmospheric loss rate prescriptions lead to a depleted atmosphere, whereas for C_tot = 0.025C_tot,⊕, the atmosphere is depleted before 6 Gyr for all atmospheric loss rate prescriptions. A mixture of outcomes occurs for C_tot = 0.25C_tot,⊕ depending on atmospheric loss rate prescription.

HPE has a moderate effect on evolution of pCO₂. Earth-like and twice-Earth HPE budgets produce similar outcomes, while a half-Earth HPE budget leads to slightly delayed outgassing and lower pCO₂ values throughout time. μ_ref and T_init each have minimal effects on the final pCO₂ value at TRAPPIST-1's current age, but high μ_ref values and low T_init values do noticeably delay the onset of outgassing.

Zieba et al. (2023) find that TRAPPIST-1c possesses pCO₂ < 0.1 bar. This helps us constrain its atmospheric history over its 7.6 ± 2.2 Gyr lifetime, particularly after it lost its H and other constituents through hydrodynamic escape. Two broad possibilities exist: (1) an atmosphere never accumulated or (2) a non-negligible atmosphere existed at some point but was later lost. In Figure 5, we see that intense outgassing occurs early in the planet's lifetime and diminishes by ∼2.5 Gyr as the mantle cools. This outgassing rate depends on geological properties of the planet, specifically the initial carbon budget C_tot. Therefore, there exists a degeneracy between the atmospheric mass-loss rate and the outgassing rate when attempting to differentiate between these two possibilities. If the atmospheric mass-loss rate is low, as in our Low models, then the outgassing rate must have been low. If the atmospheric mass-loss rate is high, then the outgassing rate could have been high or low. In Figure 6, we present pCO₂ as a function of time for the simulations, in each atmospheric loss rate prescription, that fit the observational constraint and have the highest maximum pCO₂. In the High models, we see that the highest pCO₂ value that may exist in our simulations of TRAPPIST-1c is ∼16 bar. We note that this limit is relevant for the period of TRAPPIST-1c's lifetime after it became desiccated. TRAPPIST-1c could have lost a non-negligible amount of CO₂ through drag during the hydrodynamic H loss phase as shown in Section 3.1.

Zoom In Zoom Out Reset image size

We now examine the values of each parameter that produced nominally "successful" models, or models that ended with pCO₂ < 0.1 bar at TRAPPIST-1's age, 7.6 ± 2.2 Gyr. First, note that each atmospheric loss rate prescription produced a different total number of successful models as visualized in Figure 7. 0.095% of Low models were successful, 0.495% of High models were successful, and 0.408% of Variable models were successful. Assuming our parameter space encompasses all possible values of the four geological parameters and that each combination of the four parameters is equally likely in reality, this provides evidence that a stellar mass-loss-rate function similar to that modeled in Garraffo et al. (2017) or in our variable function is more likely for TRAPPIST-1 than that from Dong et al. (2018). Within the C_tot parameter space, lower values are preferred in all atmospheric loss rate prescriptions. For Low, High, and Variable, respectively, the median C_tot values are 10^20.47, 10^20.88, and 10^20.81 mol. These are all less than Earth's estimated carbon budget C_tot,⊕ ∼ 10²² mol (which is also Venus' lower limit). In fact, while 100% of Low simulations with C_tot of 10^20.47 are successful, this percentage drops to 18% for C_tot of 10^20.54, and down to 8% for C_tot of 10^20.60. Similar drop-offs in likelihood occur in High and Variable simulations at C_tot of ∼10^21.36 and C_tot of ∼10^21.23, respectively. The median HPE values are 0.914, 1.22, and 1.22 HPE_⊕ for Low, High, and Variable, respectively. The distributions of successful models in μ_ref and T_init are relatively uniform, with slightly lower values preferred for T_init in the Low atmospheric loss rate prescription.

Zoom In Zoom Out Reset image size

Given that atmospheric loss rate, carbon budget, and pCO₂ of TRAPPIST-1c all have uncertain values at present, we ran additional simulations varying both the atmospheric loss rate and carbon budget. We chose 10 values of C_tot uniformly log-spaced over the previously used range and 10 values of constant atmospheric loss rate uniformly log-spaced between 1% and 10,000% of the constant Low value. The resulting pCO₂ values at TRAPPIST-1's age for these 100 combinations are shown in Figure 8. We repeated this for variable atmospheric loss rate of the form in Variable but with different reference loss rates at 7.6 Gyr distributed as in the previous case. This possibly more realistic set of simulations are shown in Figure 9. Given that any two of atmospheric loss rate, carbon budget, and pCO₂ at 7.6 Gyr are well known, the third quantity could be constrained using these plots.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Constraining the atmospheric evolution of TRAPPIST-1c allows us to predict the atmospheres of the other TRAPPIST-1 planets as well. Equation (3) allows us to scale the atmospheric mass-loss rate as a function of semimajor axis. The outer planets have semimajor axes a = 0.02227, 0.02925, 0.03849, 0.04683, and 0.06189 au, respectively. We assume that each planet has the same size as TRAPPIST-1c, with uniform partitioning of volatiles and radiogenic heat-producing elements across the TRAPPIST-1 system, to estimate pCO₂ for the other TRAPPIST-1 planets.

Assuming they "started" out with the same atmospheric size as desiccated TRAPPIST-1c (resulting in TRAPPIST-1c having pCO₂ < 0.1 bar at 7.6 Gyr), and scaling the High atmospheric mass-loss rate, the outer planets, d, e, f, g, and h, should have less than approximately 9, 12, 15, 16, and 16 bar, respectively. The assumption that the outer planets are of identical size to TRAPPIST-1c is a very good approximation for TRAPPIST-1g, a fair approximation for TRAPPIST-1f, and less so for the remaining smaller outer planets. These calculations suggest that the five outer planets in the TRAPPIST-1 system may possess atmospheres even if TRAPPIST-1c has little atmospheric material. Further, the assumption that the outer planets had the same amount of atmospheric material as TRAPPIST-1c after it became desiccated is also likely not a good approximation, because the process of runaway greenhouse and hydrodynamic escape would not occur to the same extent for the outer planets. This implies that they would have even larger atmospheres than calculated here. This agrees qualitatively with the recent paper by Krissansen-Totton (2023), who found that TRAPPIST-1e and f retained atmospheres in 98% of their model runs.

Our simulations assume a pure CO₂ atmosphere, which is characteristic of present-day Venus and Mars, but may not be representative of all terrestrial planet atmospheres. In fact, non-CO₂ atmospheres remain a viable possibility consistent with the observations for TRAPPIST-1c (Lincowski et al. 2023; Zieba et al. 2023). In the solar system, Titan possesses a N₂-dominated atmosphere with a non-negligible amount of CH₄. We chose a pure CO₂ atmosphere because it is common in the solar system, and is the highest mean molecular weight atmosphere possible. If any atmospheric material remains after intense escape, it should be CO₂. In this sense, assuming a pure CO₂ atmosphere represents an upper limit on the size of terrestrial planet atmospheres. A large abiotic O₂ component may be left after the early hydrodynamic escape phase, but this will escape more easily than CO₂. Additional species could exist, depending on available material during formation and the redox state of the planet's mantle. If we were to add more components, like N₂ and CH₄, we would expect these components to be lost more easily, and more complex photochemistry would operate.

While a large amount of water likely escapes from TRAPPIST-1c early, it is also reasonable to assume that some non-negligible amount of water (along with other trace gases) is outgassed with CO₂. This would replenish the atmosphere with water vapor and other gases over the lifetime of TRAPPIST-1c, which would then be removed simultaneously with CO₂. Given the same total mass-loss rate, and the observed constraint on pCO₂, the initial carbon budget after desiccation C_tot might be even lower. Thus, the initial carbon budget after desiccation that we find is an upper bound.

Another escape mechanism that could, in theory, affect the composition and size of TRAPPIST-1c's atmosphere is impact erosion, which we do not model in this manuscript. However, Raymond et al. (2022) run N-body simulations and find that high rates of impact erosion are quite unlikely in the TRAPPIST-1 system after formation because perturbations from additional objects could break the multiresonant orbital configuration of the system.

We have assumed that the planet operates under a stagnant-lid tectonic regime due to its lack of water (Tikoo & Elkins-Tanton 2017) and the fact that plate tectonics are rare in the solar system. However, TRAPPIST-1c has higher gravity than terrestrial planets in the solar system, which may increase the likelihood of subduction and therefore plate tectonics. If plate tectonics does operate on TRAPPIST-1c, previous work has suggested that outgassing rates may be higher for these "active lids" as opposed to "sluggish" or stagnant lids (Fuentes et al. 2019). Without water, recycling would still not occur, so it is plausible that the mantle's carbon reservoir would be depleted more quickly. Depending on the atmospheric mass-loss rates, this could result in a shorter atmospheric retention time.

Another geophysical consideration that is relevant to TRAPPIST-1c is the effect of tidal heating. Thus far, we have assumed that the only heating source in the mantle is from radiogenic elements. However, Dobos et al. (2019) found that the internal heat flux from tidal heating for TRAPPIST-1c is 0.62 W m⁻², comparable to that observed on Io (1–2 W m⁻²). We added this constant heating rate to our thermal evolution, which has the effect of increasing the equilibrium mantle temperature and pushing vigorous outgassing to even earlier times. We ran the edge cases of our four geological parameters again as in Figure 4, to see how tidal heating affected the evolution of pCO₂. The results are shown in Figure 10. Comparing to the results without tidal heating in Figure 4, the final pCO₂ value at TRAPPIST-1's age is relatively unchanged in each case. The largest noticeable difference is that pCO₂ increases at earlier times in all cases, due to an earlier onset of vigorous outgassing. The specific HPE value becomes unimportant in determining the evolution of pCO₂, because tidal heating becomes the dominant heating term. Overall, tidal heating does not change our conclusion that low initial carbon budgets are significantly preferred for TRAPPIST-1c to match observations.

Zoom In Zoom Out Reset image size

One more geophysical process that occurs during the magma ocean phase of a planet and subsequent solidification of the magma ocean is volatile partitioning. We assume in this manuscript that the initial carbon budget of the planet (after solidification of the magma ocean and planetary desiccation) is all located in the mantle, and it subsequently outgasses from there. However, studies of magma ocean evolution have shown that the majority of CO₂ may partition into the atmosphere (e.g., Hier-Majumder & Hirschmann 2017; Bower et al. 2022). We performed additional simulations varying the initial fraction of CO₂ f_atm in the atmosphere for different initial carbon budgets C_tot. In this case, C_tot still refers to the total amount that the planet possesses, but it is divided between the atmosphere and the mantle. Figure 11 shows the evolution of pCO₂ from these simulations. It is clear that different f_atm values slightly change the evolution of pCO₂ before 1 Gyr, but do not affect the outcome at 7.6 Gyr.

Zoom In Zoom Out Reset image size

A number of assumptions introduced by the MHD models (Dong et al. 2018) whose results we use in our simulations may impact our results. The ion escape rates are calculated in the case where the TRAPPIST-1 planets are unmagnetized and subject to maximum total wind pressure, so as to obtain upper limits on escape rates. In particles per second, Dong et al. (2018) calculate the total ion escape rates to be of the order of 10²⁷ for the TRAPPIST-1 planets. Dong et al. (2019) did similar simulations for the TRAPPIST-1 planets, varying the planetary obliquity, but only included H⁺ and O⁺ ions, so the ion escape rates are higher (10²⁸). Other recent work has calculated ion escape rates for CO₂-containing atmospheres in the case of early Earth, and found a value of ∼10²⁷ s⁻¹ at 4.0 Gyr ago (Grasser et al. 2023). This is similar to that found in Dong et al. (2018) for a much higher stellar wind pressure. This discrepancy can be explained by the fact that Grasser et al. (2023) modeled a CO₂/N₂ atmosphere, which is more susceptible to escape and included XUV-induced expansion of the upper atmosphere, which Dong et al. (2018) do not.

Dong et al. (2018) use a solar magnetogram as input to their models and scale the magnetic field strength based on observations of similar late M dwarfs (Morin et al. 2010). (Dong et al. (2018) do not report the exact value they use, but the average of the targets measured in Morin et al. (2010) is ∼500 G.) However, Garraffo et al. (2017) use the magnetogram of M6.5 dwarf GJ 3622 and 600 G for the magnetic field strength as found from Zeeman broadening of TRAPPIST-1 (Reiners & Basri 2010). These differences are likely the source of the discrepancy in stellar mass-loss rates that Dong et al. (2018) and Garraffo et al. (2017) report. Future observations should aim to better constrain these inputs. Further, the MHD models require an input of XUV luminosity to calculate photoionization and resultant stellar heating, which Dong et al. (2018) draw from Bourrier et al. (2017) to be ${F}_{\mathrm{XUV}}\,={801}_{-217}^{+436}{\rm{erg}}\,{{\rm{s}}}^{-1}{{\rm{cm}}}^{-2}$. The XUV flux of TRAPPIST-1 should also vary as a function of stellar age, which is not included in the MHD models currently.