# The Carbon-Deficient Evolution of TRAPPIST-1c

Katie E. Teixeira  Caroline V. Morley  Department of Astronomy, University of
Texas at Austin, Austin, TX, USA  Center for Planetary Systems Habitability,
University of Texas, Austin, TX, USA  Bradford J. Foley  Department of
Geosciences, Pennsylvania State University, University Park, PA, USA  Cayman
T. Unterborn  Space Instruments and Payload Department, Southwest Research
Institute, San Antonio, TX, USA

###### Abstract

Transiting planets orbiting M dwarfs provide the best opportunity to study the
atmospheres of rocky planets with current facilities. As JWST enters its
second year of science operations, an important initial endeavor is to
determine whether these rocky planets have atmospheres at all. M dwarf host
stars are thought to pose a major threat to planetary atmospheres due to their
high magnetic activity over several billion-year timescales, and might
completely strip atmospheres. Several Cycle 1 and 2 GO and GTO programs are
testing this hypothesis, observing a series of rocky planets to determine
whether they retained their atmospheres. A key case-study is TRAPPIST-1c,
which receives almost the same bolometric flux as Venus. We might, therefore,
expect TRAPPIST-1c to possess a thick, $\mathrm{CO}_{2}$-dominated atmosphere.
Instead, Zieba et al. (2023) show that TRAPPIST-1c has little to no CO2 in its
atmosphere. To interpret these results, we run coupled time-dependent
simulations of planetary outgassing and atmospheric escape to model the
evolution of TRAPPIST-1c’s atmosphere. We find that the stellar wind stripping
that is expected to occur on TRAPPIST-1c over its lifetime can only remove up
to $\sim 16$ bar of CO2, less than the modern $\mathrm{CO}_{2}$ inventory of
either Earth or Venus. Therefore, we infer that TRAPPIST-1c either formed
volatile-poor, as compared to Earth and Venus, or lost a substantial amount of
$\mathrm{CO}_{2}$ during an early phase of hydrodynamic hydrogen escape.
Finally, we scale our results for the other TRAPPIST-1 planets, finding that
the more distant TRAPPIST-1 planets may readily retain atmospheres.

exoplanets, atmospheres, habitability, TRAPPIST-1c

##  1 Introduction

Over 5000 exoplanets have now been discovered, the majority of them found
through transits with NASA’s Kepler and TESS missions. About 30% of known
exoplanets have radii less than 2 $R_{\oplus}$, classifying them as possibly
rocky (Kaltenegger, 2017; Fulton et al., 2017). While it has been possible to
constrain the radii and sometimes the masses of these potentially rocky
worlds, it has not been possible, until recently, to characterize their
atmospheres. While the atmospheres of rocky planets around the smallest stars
may be detectable with JWST (Morley et al., 2017; Lustig-Yaeger et al., 2019),
there remains great uncertainty in whether planets orbiting these small stars
will retain their atmospheres for billions of years (Lammer et al., 2007;
Tian, 2009; Cohen et al., 2015; Bourrier et al., 2017; Garraffo et al., 2017;
Becker et al., 2020; France et al., 2020; Krissansen-Totton & Fortney, 2022;
do Amaral et al., 2022).

Prior observations of rocky planets orbiting late M dwarfs have allowed us to
empirically rule out the presence of thick atmospheres around three such
planets. Spitzer, along with data from TESS and K2, allowed us to obtain
phases curves for LHS 3844b (Kreidberg et al., 2019) and K2-141b (Zieba et
al., 2022) and secondary eclipse photometry for GJ 1251b (Crossfield et al.,
2022). Each of these planets were found to have dayside temperatures
consistent with that of a bare rock with no thick atmosphere redistributing
heat to the nightside. As of mid-2023, JWST has completed its first year of
science operations, and we are now observing a series of terrestrial planets
orbiting nearby M dwarfs. Among the first observations published is a
transmission spectrum of LHS 475b, which is flat, consistent with no
atmosphere or a high-altitude cloud (Lustig-Yaeger et al., 2023). Another such
observation of the transmission spectrum of GJ 486b (Moran et al., 2023) has
found a feature in the spectrum, which is consistent with either a water-rich
atmosphere or stellar contamination on a planet orbiting an M3.5 dwarf.

Of particular interest to the exoplanet community is the M8V star TRAPPIST-1
(Gillon et al., 2017). Its seven transiting planets occupy locations which are
interior to, within, and outside of the habitable zone. This makes it a case
study in which to probe the distance from the host star where atmospheres are
retained. This is important because atmospheres are taken to be a prerequisite
for life. All of the TRAPPIST-1 planets will be observed by JWST in its first
two Cycles of operations, with a mix of transit observations (all planets) and
eclipse observations (planets b and c). The first JWST eclipse observations of
the TRAPPIST-1 system have already shown that the innermost planet,
TRAPPIST-1b, has no thick $\mathrm{CO}_{2}$ atmosphere (Greene et al., 2023),
with models from Ih et al. (2023) ruling out pure CO2 atmospheres thicker than
6.5 mbar and atmospheres containing 100 ppm CO2 thicker than 0.3 bar.

###  1.1 Atmospheric Escape from Terrestrial Planets

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
that 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 ($\mathrm{CO}_{2}$) and water ($\mathrm{H}_{2}\mathrm{O}$) then formed
by planetary outgassing as magma oceans solidified. Earth was able 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 which 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 $\mathrm{CO}_{2}$-dominated
atmosphere that we see today. Mars also has a $\mathrm{CO}_{2}$-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 (Jakosky
et al., 2018; Edberg et al., 2011). In fact, escape has affected their
atmospheres over Gyr timescales (Lammer et al., 2008).

The history of the Solar System terrestrial planets can inform our
understanding of how the atmospheres of terrestrial exoplanets evolve,
specifically close-in planets like 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 $\mathrm{H}_{2}\mathrm{O}$ photo-
dissociation 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 via the hydrodynamic wind. Therefore, a
net loss of heavier constituents like $\mathrm{CO}_{2}$ and $\mathrm{O}_{2}$
can occur. Once the hydrodynamic phase is over, we expect the remaining
atmosphere to be dominated by heavier molecules such as $\mathrm{CO}_{2}$, as
in the case of Venus.

Now, over the rest of the planet’s lifetime, $\mathrm{CO}_{2}$ 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 which exceed the escape
velocity of the planet. Jeans escape is not dominant in a
$\mathrm{CO}_{2}$-dominated atmosphere, because $\mathrm{CO}_{2}$ is too heavy
and cools efficiently through the 15 $\mu m$ band (Gordiets & Kulikov, 1985).
Therefore, we expect that long-term escape on these types of planets is
dominated by non-thermal escape, particularly in the form of stellar wind
stripping as seen in the Solar System.

Even as stellar phenemona 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) suggests 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.

###  1.2 This Work in Context

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 $\mathrm{H}_{2}\mathrm{O}$ and possesses
a Venus-like $\mathrm{CO}_{2}$-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 $\mu m$
$\mathrm{CO}_{2}$ absorption band. The data from these observations (Zieba et
al., 2023) point to the absence of a thick $\mathrm{CO}_{2}$ atmosphere,
ruling out pure $\mathrm{CO}_{2}$ atmospheres with surface pressures $\geq
0.1$ bar. Zieba et al. (2023) find that a bare rock surface or a thin
$\mathrm{O}_{2}$-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 $\mathrm{CO}_{2}$ 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$\%$ $\mathrm{H}_{2}\mathrm{O}$
abundance is consistent with the observations.

Assuming $\mathrm{CO}_{2}$ has dominated TRAPPIST-1c’s atmosphere for most of
its history, the constraint on $\mathrm{CO}_{2}$ 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) which made predictions for the atmospheres
of TRAPPIST-1 planets. Krissansen-Totton & Fortney (2022) includes a range of
atmospheric species, including $\mathrm{H}_{2}\mathrm{O}$, $\mathrm{CO}_{2}$,
and $\mathrm{O}_{2}$, and runs coupled simulations starting in the magma ocean
phase all the way to present day. We first calculate the amount of
$\mathrm{CO}_{2}$ potentially lost during the hydrodynamic escape phase, and
then run coupled time-dependent simulations of outgassing and escape of
$\mathrm{CO}_{2}$ starting after magma ocean solidification and atmospheric
desiccation. Our simulations, therefore, track the amount of $\mathrm{CO}_{2}$
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 $\mathrm{CO}_{2}$ 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 $\mathrm{CO}_{2}$ 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.

Figure 1: Comparing TRAPPIST-1 planets to Solar System planets (left) and
likely rocky exoplanets already observed or planned to be observed by JWST in
Cycle 1 and 2 (right). Incident flux is plotted against escape velocity, so
that planets in the top left are less likely to retain atmospheres, and
planets in the bottom right are more likely to retain atmospheres. Stellar
temperature, represented by color, may also play a role in atmosphere
retention. Black Xs are placed on the labeled exoplanets which have already
been observed to have no thick atmosphere. The gray line in the right panel is
the “cosmic shoreline” ($\mathrm{F_{inc}}\propto\mathrm{v_{esc}}^{4}$)
presented in Zahnle & Catling (2017).

##  2 Methods

We model the atmosphere of TRAPPIST-1c as a coupled system that includes the
interior of the planet and the stellar environment. We encode observationally-
constrained properties of TRAPPIST-1c and its host star into simulations which
evolve the atmosphere through the planet’s lifetime. In our simulations,
planetary outgassing delivers $\mathrm{CO}_{2}$ to the atmosphere and is based
on geophysical phenomena observed on Earth, except without plate tectonics, as
on present-day Venus. Atmospheric escape removes $\mathrm{CO}_{2}$ from the
atmosphere and is based on the stellar-wind-induced ion escape of Venus and
Mars observed and modeled by magnetohydrodynamic (MHD) simulations, accounting
for the differing properties of M dwarfs like TRAPPIST-1. Coupling these
processes then allows us to track the amount of $\mathrm{CO}_{2}$ in the
atmosphere at any given time.

###  2.1 Stellar Properties

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_{\star}=0.0898\pm 0.0023\ M_{\odot}$, radius $R_{\star}=0.1192\pm 0.0013\
R_{\odot}$, effective temperature $T_{\mathrm{eff}}=2,566\pm 26$ K, and
bolometric luminosity $L_{\mathrm{bol}}=0.000553\pm 0.000018\ L_{\odot}$.

We use the age of TRAPPIST-1 presented in Burgasser & Mamajek (2017). They
conclude that TRAPPIST-1 has an age of $7.6\pm 2.2$ Gyr. The rotation period
of TRAPPIST-1 is $3.30\pm 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 Section 2.3 and Section 4.4.

###  2.2 Planet Properties

In this manuscript, we focus on the second-closest planet to the host star,
TRAPPIST-1c. TRAPPIST-1c has mass $M_{p}=1.308\pm 0.056\ M_{\oplus}$, radius
$R_{p}=1.097\begin{subarray}{c}+0.014\\\ -0.012\end{subarray}\ R_{\oplus}$,
and semi-major axis $a=0.01580\pm 0.00013$ AU (Agol et al., 2021). These
planetary properties are displayed with the stellar properties in Table 1.
Mass and radius constraints provide an estimate for the density of
TRAPPIST-1c, approximately $\rho=5.45\ \mathrm{g}/\mathrm{cm}^{3}$, which is
almost exactly Earth-like. The surface gravity of TRAPPIST-1c is $10.65\pm
0.42\ \mathrm{m}/\mathrm{s}^{2}$, slightly larger than Earth’s. Given the
semi-major axis and stellar luminosity, we calculate the incident bolometric
flux on TRAPPIST-1c to be approximately
$F_{\mathrm{inc}}=2.22F_{\mathrm{inc},\oplus}$, which is near the value for
Venus ($1.91F_{\mathrm{inc},\oplus}$).

$\mathrm{Parameter}$ | $\mathrm{Value}$
---|---
Star - TRAPPIST-1 |
Mass $M_{\star}$ | $0.0898\pm 0.0023\ M_{\odot}$
Radius $R_{\star}$ | $0.1192\pm 0.0013\ R_{\odot}$
Effective temperature $T_{\mathrm{eff}}$ |  $2566\pm 26$ K
Bolometric luminosity $L_{\mathrm{bol}}$ | $0.000553\pm 0.000018\ L_{\odot}$
Age |  $7.6\pm 2.2$ Gyr
Rotation period $P_{\mathrm{rot}}$ |  $3.30\pm 0.14$ days
Planet - TRAPPIST-1c |
Mass $M_{p}$ | $1.308\pm 0.056\ M_{\oplus}$
Radius $R_{p}$ | $1.097\begin{subarray}{c}+0.014\\\ -0.012\end{subarray}\ R_{\oplus}$
Semi-major axis $a$ |  $0.01580\pm 0.00013$ AU
Table 1: Relevant properties of TRAPPIST-1 and TRAPPIST-1c as cited in the
text.

The aforementioned properties of TRAPPIST-1c are compared to that 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 GO and GTO 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
stellar material bombardment 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
higher stellar temperature cases where high-activity pre-main sequence phases
are shorter. The gray line in the right panel is the “cosmic shoreline”
($\mathrm{F_{inc}}\propto\mathrm{v_{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.

###  2.3 Atmospheric Evolution

In our simulations, the $\mathrm{CO}_{2}$ atmosphere of TRAPPIST-1c evolves
through two mechanisms: planetary outgassing and atmospheric escape.
Outgassing is a result of planet interior processes 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 magnetohydrodynamic (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 (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 stagnant lid planet interiors to melting, crustal growth,
weathering, and $\mathrm{CO}_{2}$ 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 the atmosphere,
by metamorphic outgassing.

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

$V_{man}\rho
c_{p}\frac{dT_{p}}{dt}=Q_{man}-A_{man}F_{man}-f_{m}\rho_{m}(c_{p}\Delta
T_{m}+L_{m})$ (1)

where $V_{man}$ is the volume of convecting mantle beneath the stagnant lid,
$\rho$ 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}$ the heat flux from
the convecting mantle at the base of the stagnant lid, $f_{m}$ is the
volumetric melt production, $\rho_{m}$ is the density of the mantle melt,
$L_{m}$ is the latent heat of fusion of the mantle, and $\Delta 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
mantle thermal evolution.

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 atmopshere:

$\frac{dR_{man}}{dt}=-\frac{f_{m}R_{man}[1-(1-\phi)^{1/D_{CO_{2}}}]}{\phi(V_{man}+V_{lid})}$
(2)

where $R_{man}$ is the mantle carbon reservoir, $\phi$ is the melt fraction,
$D_{CO_{2}}$ is the distribution coefficient of $\mathrm{CO}_{2}$, and
$V_{lid}$ is the volume of the stagnant lid. This assumes that a certain
fraction of $\mathrm{CO}_{2}$ enters the melt phase and a fraction of the melt
degasses its $\mathrm{CO}_{2}$ 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 $\mathrm{CO}_{2}$ is
returned to the crust, and, therefore, no $\mathrm{CO}_{2}$ is recycled into
the mantle, so metamorphic outgassing and carbon foundering do not operate.
This implies that all $\mathrm{CO}_{2}$ which is outgassed will simply build
up in the atmosphere. This is captured in Figure 2 with the major processes
that potentially move $\mathrm{CO}_{2}$ into and out of TRAPPIST-1c’s
atmosphere.

Figure 2: Processes that move atmospheric species into and out of a stagnant
lid planet’s atmosphere. Without liquid water, metamorphic outgassing and
carbon foundering do not occur.

####  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 which could protect
their atmospheres. The stellar magnetic field produces small, induced
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
$\mathrm{CO}_{2}$. The model then solves the MHD equations for the 4 ion
fluids that exist in this case: $\mathrm{H}^{+}$, $\mathrm{O}^{+}$,
$\mathrm{O}_{2}^{+}$, and $\mathrm{CO}_{2}^{+}$. The physical processes that
are captured by solving these equations include stellar wind charge exchange,
ion pickup, and ion sputtering, all which act to remove atmospheric material
by accelerating it away from the planet. Therefore, Dong et al. (2018) is 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\cdot
10^{8}$ kg/s, the atmospheric loss rate from TRAPPIST-1c is calculated to be
$\dot{M}_{\mathrm{atm}}=57$ kg/s. Dong et al. (2018) finds that the
atmospheric loss rate numerically calculated in their simulations scales with
planet radius, semi-major axis, and stellar mass-loss rate as follows:

$\dot{M}_{\mathrm{atm}}\propto\left({\frac{R_{p}}{a}}\right)^{2}\dot{M}_{\star}$
(3)

This scaling relation allows us to test different stellar mass-loss rates from
the literature. For example, Garraffo et al. (2017) cites the stellar mass-
loss rate of TRAPPIST-1 as $\dot{M}_{\star}=1.89\cdot 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$ 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 M dwarf
empirical relationships. 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:

$P_{\mathrm{rot}}[\mathrm{d}]=\frac{\mathrm{age[Gyr]}-0.012\pm 0.221}{0.061\pm
0.002}$ (4)
$L\mathrm{{}_{X}[erg/s]}=\begin{cases}C_{\mathrm{sat}}P_{\mathrm{rot}}^{\beta_{\mathrm{sat}}},&{P_{\mathrm{rot}}}\leq{P_{\mathrm{rot,sat}}}\\\
C_{\mathrm{unsat}}P_{\mathrm{rot}}^{\beta_{\mathrm{unsat}}},&{P_{\mathrm{rot}}}>{P_{\mathrm{rot,sat}}}\\\
\end{cases}$ (5)

where $C_{n}=(L_{\mathrm{X},n}/P_{\mathrm{rot}}^{\beta_{n}})$ with
n=(sat,unsat), $\beta_{\mathrm{sat}}=-0.19\pm 0.11$,
$\beta_{\mathrm{unsat}}=-3.52\pm 0.02$, $P_{\mathrm{rot,sat}}=33.7\pm 4.5$
[d], and log$(L_{\mathrm{X,sat}}(P_{\mathrm{rot}}=1[\mathrm{d}]))=28.54\pm
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):

$\frac{\dot{M_{\star}}}{\mathrm{SA}}\left[\frac{\dot{M}_{\odot}}{\mathrm{SA_{\odot}}}\right]=10^{(0.77\pm
0.04\cdot\mathrm{log}_{10}(\mathrm{F_{X,surface}})-3.42)}$ (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 M dwarf scaling
relations (Variable). These three different cases are shown in Figure 3.

Figure 3: Three different planetary atmospheric loss rate prescriptions as a
function of time used in this work. The Low rate is taken from Dong et al.
(2018), the High rate is derived from Garraffo et al. (2017) using Equation 3,
and the Variable rate is calculated from M dwarf scaling relations referenced
in Section 2.3.2.

###  2.4 Simulations

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

  * •

$\mathrm{C_{tot}}$: Initial carbon budget of the planet. This is the amount of
$\mathrm{CO}_{2}$ in moles that the planet initially possesses after magma
ocean solidification and atmospheric desiccation, all located in the mantle
unless otherwise specified. Earth’s estimated budget is $10^{22}$ mol per
Foley & Smye (2018). Venus’ $\mathrm{C_{tot}}$ must be at least $10^{22}$ mol,
given that is approximately how much exists in the atmosphere at present.

  * •

HPE: Initial radiogenic heat-producing element ($\mathrm{U}^{238}$,
$\mathrm{U}^{235}$, K, and Th) budget in the mantle. These elements act to
heat the mantle and may prolong outgassing.

  * •

$\mu_{\mathrm{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.

  * •

$T_{\mathrm{init}}$: Initial mantle temperature. This is the starting
temperature of the mantle before significant cooling occurs and 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:

$\log_{10}\left(\frac{\mathrm{C}_{\mathrm{tot}}}{\mathrm{C}_{\mathrm{tot},\oplus}}\right)\sim\mathscr{U}\left(\log_{10}(0.025),\log_{10}(2.5)\right)$
(7)
$\frac{\mathrm{HPE}}{\mathrm{HPE_{\oplus}}}\sim\mathscr{U}\left(0.5,2\right)$
(8)
$\log_{10}\left(\frac{\mu_{\mathrm{ref}}}{\mathrm{Pa}\cdot\mathrm{s}}\right)\sim\mathscr{U}\left(20,22\right)$
(9)
$\frac{T_{\mathrm{init}}}{\mathrm{K}}\sim\mathscr{U}\left(1700,2000\right)$
(10)

Here $\mathrm{C}_{\mathrm{tot},\oplus}=10^{22}$ per an estimate of Earth’s
mantle and surface reservoirs from Sleep & Zahnle (2001). Our chosen range for
$\mathrm{C}_{\mathrm{tot}}$ represents a range of $10^{-4}$ to $10^{-2}$ wt%
concentration of $\mathrm{CO}_{2}$ in the mantle. We test a range of initial
radiogenic heat-producing element budget that spans 50% to 200% of Earth’s
budget, about 80 to 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.

##  3 Results

###  3.1 Early Hydrodynamic $\mathrm{CO}_{2}$ Loss

During the hydrodynamic hydrogen loss phase, we emphasize that it is possible
for $\mathrm{CO}_{2}$ 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
$\mathrm{CO}_{2}$ lost during it per Odert et al. (2018). For a hydrogen-
dominated atmosphere with a $\mathrm{CO}_{2}$ component, the XUV-driven
hydrodynamic flux of hydrogen in $\mathrm{kg}/\mathrm{s/m}^{2}$ is given by:

$\mathrm{F_{H}}=\frac{\beta^{2}\eta\mathrm{F_{XUV}}}{4\Delta\Phi(m_{\mathrm{H}}+m_{\mathrm{CO}_{2}}f_{\mathrm{CO}_{2}}x_{\mathrm{CO}_{2}})}$
(11)

Here $\beta=\mathrm{R_{XUV}}/\mathrm{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, $\eta$ is the efficiency of heating, $\mathrm{F_{XUV}}$ is the XUV
flux, $\Delta\Phi=\mathrm{GM_{p}}/\mathrm{R_{p}}$, $m_{\mathrm{H}}$ is the
mass of a hydrogen atom, $m_{\mathrm{CO}_{2}}$ is the mass of a
$\mathrm{CO}_{2}$ molecule, $f_{\mathrm{CO}_{2}}$ is the mixing ratio of
$\mathrm{CO}_{2}$ to H, and $x_{\mathrm{CO}_{2}}$ is the fractionation factor
of $\mathrm{CO}_{2}$. The fractionation factor $x_{\mathrm{CO}_{2}}$ is
calculated by the following equation:

$x_{\mathrm{CO}_{2}}=1-\frac{g(m_{\mathrm{H}}-m_{\mathrm{CO}_{2}})b(T)}{\mathrm{F_{H}}k_{B}T(1+f_{\mathrm{CO}_{2}})}$
(12)

where $b$ is the binary diffusion parameter and $T$ is the upper atmosphere
temperature. $b=8.4\times 10^{19}T^{0.6}\ \mathrm{m}^{-1}\ \mathrm{s}^{-1}$
for H and $\mathrm{CO}_{2}$. The flux of $\mathrm{CO}_{2}$ hydrodynamically
dragged is given by:

$\mathrm{F_{CO_{2}}}=\mathrm{F_{H}}f_{\mathrm{CO}_{2}}x_{\mathrm{CO}_{2}}$
(13)

In order to track the amount of H and $\mathrm{CO}_{2}$ over time, Equation 11
and 12 must first be solved for $\mathrm{F_{H}}$. Then, this equation, along
with Equation 13, must be integrated using a differential equation solver.
When performing these calculations, we assumed $\beta=1$, $\eta=0.15$ (Salz et
al., 2015; Kubyshkina et al., 2018), and $\mathrm{F_{XUV}}$ is the sum of EUV
and X-ray flux:

$\mathrm{F_{XUV}}=\mathrm{F_{EUV}}+\mathrm{F_{X}}$ (14)

We calculate X-ray flux as in Equation 5 and EUV flux per 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 to 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 $\sim$ 1000 years and it could have lost about 0.1
bar of $\mathrm{CO}_{2}$. Starting with 1 EO, the hydrodynamic phase could
last for about 100 thousand years, and about 10 bar of $\mathrm{CO}_{2}$ could
be lost. With 100 EOs, the hydrodynamic phase could last for $\sim$ 10 Myr,
and about 1000 bar of $\mathrm{CO}_{2}$ could be lost. We see that the length
of the hydrodynamic phase is extremely short compared to the lifetime of
TRAPPIST-1c in all cases, but the amount of $\mathrm{CO}_{2}$ lost during this
phase can vary widely.

###  3.2 Evolution Over Parameter Space

After calculating the rapid hydrodynamic $\mathrm{CO}_{2}$ loss, we focus on
the $\mathrm{CO}_{2}$ evolution 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
$\mathrm{CO}_{2}$ 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 $\mathrm{CO}_{2}$
present in the atmosphere, which can be represented by mass in kg or surface
pressure (p$\mathrm{CO}_{2}$) in bar.

The modeled evolution of p$\mathrm{CO}_{2}$ on TRAPPIST-1c is shown for the
edge cases of our varied parameters in Figure 4. The baseline parameters are
$\mathrm{C_{tot}}[\oplus]=0.25$, $\mu_{\mathrm{ref}}=1.0\times 10^{21}\
\mathrm{Pa}\cdot\mathrm{s}$, $\mathrm{T_{init}}=2000\ \mathrm{K}$ and
$\mathrm{HPE}[\oplus]=1$ in each panel. All simulations are characterized by a
sharp increase in p$\mathrm{CO}_{2}$ in the first 1 Gyr of evolution due to
early vigorous outgassing, as shown in Figure 5. Some simulations result in
p$\mathrm{CO}_{2}$ that is relatively constant over the rest of TRAPPIST-1c’s
lifetime; others result in order of magnitude losses in p$\mathrm{CO}_{2}$.
Specifically, the Low atmospheric loss rate prescription results in a minimal
decrease in p$\mathrm{CO}_{2}$ in most cases. The High prescription leads to
significant decreases in p$\mathrm{CO}_{2}$ that deplete the atmosphere by
TRAPPIST-1’s age in most cases. The Variable prescription produces
p$\mathrm{CO}_{2}$ curves that decrease slightly in the 1-6 Gyr range but then
level out as stellar activity decreases.

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

HPE has a moderate effect on p$\mathrm{CO}_{2}$ evolution. Earth-like and
twice Earth HPE budgets produce similar outcomes, while a half Earth HPE
budget leads to slightly delayed outgassing and lower p$\mathrm{CO}_{2}$
values throughout time. $\mu_{\mathrm{ref}}$ and $T_{\mathrm{init}}$ each have
minimal effects on the final p$\mathrm{CO}_{2}$ value at TRAPPIST-1’s current
age but high $\mu_{\mathrm{ref}}$ values and low $T_{\mathrm{init}}$ values do
noticeably delay the onset of outgassing.

Figure 4: Evolution of atmospheric surface pressure with time in edge cases of
geological parameter space with three different atmospheric loss rate
prescriptions. The baseline parameters in each panel are
$\mathrm{C_{tot}}[\oplus]=0.25$, $\mu_{\mathrm{ref}}=1.0\times 10^{21}\
\mathrm{Pa}\cdot\mathrm{s}$, $\mathrm{T_{init}}=2000\ \mathrm{K}$ and
$\mathrm{HPE}[\oplus]=1$.

###  3.3 Observational Constraints on Evolution

Zieba et al. (2023) find that TRAPPIST-1c possesses p$\mathrm{CO}_{2}$ $<$ 0.1
bar. This helps us constrain its atmospheric history over its $7.6\pm 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 $\sim$ 2.5 Gyr as the mantle cools. This
outgassing rate depends on geological properties of the planet, specifically
the initial carbon budget $\mathrm{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
p$\mathrm{CO}_{2}$ as a function of time for the simulations, in each
atmospheric loss rate prescription, that fit the observational constraint and
have the highest maximum p$\mathrm{CO}_{2}$. In the High models, we see that
the highest p$\mathrm{CO}_{2}$ value that may exist in our simulations of
TRAPPIST-1c is $\sim$ 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 $\mathrm{CO}_{2}$ through drag during the
hydrodynamic H loss phase as shown in Section 3.1.

Figure 5: Outgassing rate as a function of age for three different carbon
budgets $\mathrm{C_{tot}}$. Other parameters are set at
$\mu_{\mathrm{ref}}=1.0\times 10^{21}\ \mathrm{Pa}\cdot\mathrm{s}$,
$\mathrm{T_{init}}=2000\ \mathrm{K}$ and $\mathrm{HPE}[\oplus]=1$. Figure 6:
p$\mathrm{CO}_{2}$ as a function of time for the simulations, in each
atmospheric loss rate prescription, that fit the observational constraint
(p$\mathrm{CO}_{2}<0.1$ bar at 7.6 Gyr) and have the highest maximum
p$\mathrm{CO}_{2}$.

We now examine the values of each parameter that produced nominally
“successful” models, or models that ended with p$\mathrm{CO}_{2}$ $<$ 0.1 bar
at TRAPPIST-1’s age, $7.6\pm 2.2$ Gyr. First, note that each atmospheric loss
rate prescription produced different total numbers 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 are
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 $\mathrm{C_{tot}}$ parameter space, lower values are
preferred in all atmospheric loss rate prescriptions. For Low, High, and
Variable, respectively, the median $\mathrm{C_{tot}}$ values are $10^{20.47}$,
$10^{20.88}$ and $10^{20.81}$ moles. These are all less than Earth’s estimated
carbon budget $\mathrm{C_{tot,\oplus}}\sim 10^{22}$ moles (which is also
Venus’ lower limit). In fact, while 100$\%$ of Low simulations with
$\mathrm{C_{tot}}$ of $10^{20.47}$ are successful, this percentage drops to
18$\%$ for $\mathrm{C_{tot}}$ of $10^{20.54}$, and down to 8$\%$ for
$\mathrm{C_{tot}}$ of $10^{20.60}$. Similar drop-offs in likelihood occur in
High and Variable simulations at $\mathrm{C_{tot}}$ of $\sim 10^{21.36}$ and
$\mathrm{C_{tot}}$ of $\sim 10^{21.23}$, respectively. The median HPE values
are 0.914, 1.22, and 1.22 $\mathrm{HPE}_{\oplus}$ for Low, High, and Variable,
respectively. The distributions of successful models in $\mu_{\mathrm{ref}}$
and $T_{\mathrm{init}}$ are relatively uniform, with slightly lower values
preferred for $T_{\mathrm{init}}$ in the Low atmospheric loss rate
prescription.

Figure 7: Histograms of parameter values of simulations (in each atmospheric
loss rate prescription) that result in p$\mathrm{CO}_{2}$ $<$ 0.1 bar at
TRAPPIST-1’s measured age, 7.6 Gyr, shown in blue. The error in TRAPPIST-1’s
age, $\pm 2.2$ Gyr, translates to some variation in these histograms which is
shown in yellow and red. The median of each parameter distribution that
results in p$\mathrm{CO}_{2}$ $<$ 0.1 bar at 7.6 Gyr is plotted as a vertical
dashed line.

Given that atmospheric loss rate, carbon budget, and p$\mathrm{CO}_{2}$ 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 $\mathrm{C_{tot}}$ uniformly log-spaced over the previously used
range and 10 values of constant atmospheric loss rate uniformly log-spaced
between 1% and 10000% of the constant Low value. The resulting
p$\mathrm{CO}_{2}$ 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 p$\mathrm{CO}_{2}$ at 7.6 Gyr are well known, the
third quantity could be constrained using these plots.

Figure 8: Atmospheric surface pressure as a function of carbon budget and loss
rate, assuming loss rate is constant over time. Figure 9: Atmospheric surface
pressure as a function of carbon budget and reference loss rate at 7.6 Gyr,
assuming loss rate is variable (left). Atmospheric mass-loss rate versus time
is shown for the 10 variable cases (right). Note that the variable function is
multiplied by a single factor rather than translated.

##  4 Discussion

###  4.1 Scaling Our Results for TRAPPIST-1d, e, f, g, and h

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 function of semi-major axis. The outer
planets have semi-major 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 volatile and radiogenic heat-producing element
partitioning across the TRAPPIST-1 system, to estimate p$\mathrm{CO}_{2}$ for
the other TRAPPIST-1 planets.

Assuming they “started” out with the same atmospheric size as desiccated
TRAPPIST-1c (resulting in TRAPPIST-1c having p$\mathrm{CO}_{2}<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 identically-sized 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 5 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.

###  4.2 Atmospheric Processes

Our simulations assume a pure $\mathrm{CO}_{2}$ atmosphere, which is
characteristic of present-day Venus and Mars, but may not be representative of
all terrestrial planet atmospheres. In fact, non-CO2 atmospheres remain a
viable possibility consistent with the observations for TRAPPIST-1c (Zieba et
al., 2023; Lincowski et al., 2023). In the solar system, Titan possesses a
$\mathrm{N}_{2}$-dominated atmosphere with a non-negligible amount of
$\mathrm{CH}_{4}$. We chose a pure $\mathrm{CO}_{2}$ 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 $\mathrm{CO}_{2}$. In this sense, assuming a pure
$\mathrm{CO}_{2}$ atmosphere represents an upper limit on the size of
terrestrial planet atmospheres. A large abiotic $\mathrm{O}_{2}$ component may
be left after the early hydrodynamic escape phase, but this will escape more
easily than $\mathrm{CO}_{2}$. 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 $\mathrm{N}_{2}$ and
$\mathrm{CH}_{4}$, 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 $\mathrm{CO}_{2}$. This would replenish
the atmosphere with water vapor and other gases over the lifetime of
TRAPPIST-1c, which would then be removed simultaneously with
$\mathrm{CO}_{2}$. Given the same total mass-loss rate, and the observed
constraint on p$\mathrm{CO}_{2}$, the initial carbon budget after desiccation
$\mathrm{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 multi-resonant orbital configuration of the system.

###  4.3 Geophysical Processes

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 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
$\mathrm{W/m}^{2}$, comparable to that observed on Io (1-2
$\mathrm{W/m}^{2}$). We added this a 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 $\mathrm{pCO}_{2}$ evolution. The results are shown
in Figure 10. Comparing to the results without tidal heating in Figure 4, the
final $\mathrm{pCO}_{2}$ value at TRAPPIST-1’s age is relatively unchanged in
each case. The largest noticeable difference is that $\mathrm{pCO}_{2}$
increases at earlier times in all cases, due to an earlier onset of vigorous
outgassing. The specific HPE value becomes unimportant in determining the
$\mathrm{pCO}_{2}$ evolution, 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.

Figure 10: Evolution of atmospheric surface pressure with time in edge cases
of geological parameter space with three different atmospheric loss rate
prescriptions with constant tidal heating rate added. The baseline parameters
in each panel are $\mathrm{C_{tot}}[\oplus]=0.25$,
$\mu_{\mathrm{ref}}=1.0\times 10^{21}\ \mathrm{Pa}\cdot\mathrm{s}$,
$\mathrm{T_{init}}=2000\ \mathrm{K}$ and $\mathrm{HPE}[\oplus]=1$.

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 magma ocean solidification 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
$\mathrm{CO}_{2}$ may partition into the atmosphere (e.g. Hier-Majumder &
Hirschmann 2017; Bower et al. 2022). We performed additional simulations
varying the fraction of $\mathrm{CO}_{2}$ $f_{\mathrm{atm}}$ starting out in
the atmosphere for different initial carbon budgets $\mathrm{C_{tot}}$. In
this case, $\mathrm{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 $\mathrm{pCO}_{2}$ evolution from these simulations. It is clear
that different $f_{\mathrm{atm}}$ values slightly change the evolution of
$\mathrm{pCO}_{2}$ before 1 Gyr, but do not affect the outcome at 7.6 Gyr.

Figure 11: Evolution of atmospheric surface pressure with time, varying
$\mathrm{C_{tot}}$ and $f_{\mathrm{atm}}$, the fraction of $\mathrm{C_{tot}}$
placed in the atmosphere. The baseline parameters are
$\mu_{\mathrm{ref}}=1.0\times 10^{21}\ \mathrm{Pa}\cdot\mathrm{s}$,
$\mathrm{T_{init}}=2000\ \mathrm{K}$ and $\mathrm{HPE}[\oplus]=1$.

###  4.4 Star-Planet Interaction

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, as to obtain upper limits on escape
rates. In particles per second, Dong et al. (2018) calculates the total ion
escape rates to be on the order of $10^{27}$ for the TRAPPIST-1 planets. Dong
et al. (2019) did similar simulations for the TRAPPIST-1 planets, varying the
planetary obliquity, but only included $\mathrm{H}^{+}$ and $\mathrm{O}^{+}$
ions, so the ion escape rates are higher ($10^{28}$). Other recent work has
calculated ion escape rates for $\mathrm{CO}_{2}$-containing atmospheres in
the case of early Earth, and found a value of $\sim 10^{27}$/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
$\mathrm{CO}_{2}$/$\mathrm{N}_{2}$ atmosphere which is more susceptible to
escape and included XUV-induced expansion of the upper atmosphere which Dong
et al. (2018) does not.

Dong et al. (2018) uses a solar magnetogram as input to their models and
scales the magnetic field strength based on observations of similar late M
dwarfs (Morin et al., 2010) (Dong et al. (2018) don’t report the exact value
they use, but the average of the targets measured in Morin et al. (2010) is
$\sim 500$ G). However, Garraffo et al. (2017) uses the magnetogram of M6.5
dwarf GJ3622 and 600 G for the magnetic field strength as found by 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) draws from Bourrier et al. (2017) to be
$F_{\mathrm{XUV}}=801\begin{subarray}{c}+436\\\ -217\end{subarray}$
ergs/s/$\mathrm{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.

##  5 Conclusions

In this manuscript, we have presented simulations of outgassing and escape
that evolve the atmosphere of TRAPPIST-1c, assuming it is
$\mathrm{CO}_{2}$-dominated. Given our assumptions about atmospheric escape
and geological phenomena, and based on the constraint of p$\mathrm{CO}_{2}$
$<$ 0.1 bar from observations, we find that:

  1. 1.

Long-term stellar wind stripping is not efficient enough to remove a large
$\mathrm{CO}_{2}$ atmosphere from TRAPPIST-1c. Specifically,

     * •

The median initial carbon budget after desiccation ($\mathrm{C_{tot}}$) is
significantly less than modern Earth’s carbon budget and Venus’ lower limit,
$10^{22}$ mol, in all atmospheric loss rate prescriptions.

     * •

The surface pressure of $\mathrm{CO}_{2}$ in TRAPPIST-1c’s past, following
early hydrodynamic loss, could have only been as high as $\sim 16$ bar.

  2. 2.

TRAPPIST-1c must have either formed volatile-poor as compared to Earth and
Venus, or lost a substantial amount of $\mathrm{CO}_{2}$ during the early
hydrodynamic escape phase which requires a high initial water inventory – at
least $\sim$ 10 Earth oceans.

  3. 3.

The outer TRAPPIST-1 planets may retain substantial atmospheres which are
$\sim 2$ orders of magnitude larger than that of TRAPPIST-1c.

Our results show that TRAPPIST-1c’s low constraint on p$\mathrm{CO}_{2}$ today
traces back to a lower carbon budget than Earth and Venus after desiccation.
This divergence in evolution could be due to a difference in total volatile
inventory (Scenario 1) which may be due to variations in protoplanetary disk
composition and size, or it might relate to orbital dynamics, such as the
scattering of volatile-rich planetesimals due to giant planets during
formation (Raymond et al., 2004). An equally valid scenario (2) exists where
TRAPPIST-1c’s volatile inventory was similar to that of Earth and Venus, and
the high-energy radiation environment of TRAPPIST-1 was harsh enough to
deplete it in its first $10^{5}$–$10^{7}$ years.

More JWST observations in various wavelength bands will give us further
insight into the TRAPPIST-1 planets’ atmospheric sizes and compositions. While
we expect $\mathrm{CO}_{2}$ to be the most resistant species to intense
escape, other species could represent signatures of atmospheric evolution. It
is also plausible that TRAPPIST-1c and its neighbors do, in fact, lack
atmospheres; the observations presented here from Zieba et al. (2023) only
provide an upper limit. Regardless, evolutionary models, which couple stellar,
atmospheric, and geological processes, like that which we have presented here,
will be key to understanding how these planets became what they are today.

##  6 Code Availability

We make our code publicly available. The software is available on GitHub
(https://github.com/katieteixeira/atmospheric_evolution) under a MIT License
and version 1.0.0 is archived in Zenodo (Teixeira & Foley, 2023). It can be
used to run simulations, save and load data, and make figures.

We thank Brendan Bowler, Adam Kraus, Amber Medina, David Wilson, and Cynthia
Froning for fruitful discussions that aided the development of this project
and manuscript. This work was supported through a Student Research Award in
Planetary Habitability from the UT Center for Planetary Systems Habitability.
Support for programs HST-AR-17025.001-A, HST-AR-15805.001-A, and JWST-
GO-02304.002-A were provided by NASA through a grant from the Space Telescope
Science Institute, which is operated by the Associations of Universities for
Research in Astronomy, Incorporated, under NASA contract NAS5-26555. C.V.M.
acknowledges support from the Alfred P. Sloan Foundation under grant number
FG-2021-16592.

## References

  * Agol et al. (2021) Agol, E., Dorn, C., Grimm, S. L., et al. 2021, , 2, 1, doi: 10.3847/PSJ/abd022
  * Becker et al. (2020) Becker, J., Gallo, E., Hodges-Kluck, E., Adams, F. C., & Barnes, R. 2020, AJ, 159, 275, doi: 10.3847/1538-3881/ab8fb0
  * Botelho et al. (2019) Botelho, R. B., Milone, A. d. C., Meléndez, J., et al. 2019, MNRAS, 482, 1690, doi: 10.1093/mnras/sty2791
  * Boudreaux et al. (2022) Boudreaux, E. M., Newton, E. R., Mondrik, N., Charbonneau, D., & Irwin, J. 2022, ApJ, 929, 80, doi: 10.3847/1538-4357/ac5cbf
  * Bourrier et al. (2017) Bourrier, V., Ehrenreich, D., Wheatley, P. J., et al. 2017, A&A, 599, L3, doi: 10.1051/0004-6361/201630238
  * Bower et al. (2022) Bower, D. J., Hakim, K., Sossi, P. A., & Sanan, P. 2022, , 3, 93, doi: 10.3847/PSJ/ac5fb1
  * Burgasser & Mamajek (2017) Burgasser, A. J., & Mamajek, E. E. 2017, The Astrophysical Journal, 845, 110, doi: 10.3847/1538-4357/aa7fea
  * Cohen et al. (2015) Cohen, O., Ma, Y., Drake, J. J., et al. 2015, ApJ, 806, 41, doi: 10.1088/0004-637X/806/1/41
  * Costa et al. (2006) Costa, E., Mendez, R. A., Jao, W. C., et al. 2006, The Astronomical Journal, 132, 1234
  * Crossfield et al. (2022) Crossfield, I. J. M., Malik, M., Hill, M. L., et al. 2022, ApJ, 937, L17, doi: 10.3847/2041-8213/ac886b
  * do Amaral et al. (2022) do Amaral, L. N. R., Barnes, R., Segura, A., & Luger, R. 2022, ApJ, 928, 12, doi: 10.3847/1538-4357/ac53af
  * Dobos et al. (2019) Dobos, V., Barr, A. C., & Kiss, L. L. 2019, A&A, 624, A2, doi: 10.1051/0004-6361/201834254
  * Dong et al. (2019) Dong, C., Huang, Z., & Lingam, M. 2019, ApJ, 882, L16, doi: 10.3847/2041-8213/ab372c
  * Dong et al. (2018) Dong, C., Jin, M., Lingam, M., et al. 2018, Proceedings of the National Academy of Science, 115, 260, doi: 10.1073/pnas.1708010115
  * Edberg et al. (2011) Edberg, N. J. T., Nilsson, H., Futaana, Y., et al. 2011, Journal of Geophysical Research (Space Physics), 116, A09308, doi: 10.1029/2011JA016749
  * Engle & Guinan (2018) Engle, S. G., & Guinan, E. F. 2018, Research Notes of the American Astronomical Society, 2, 34, doi: 10.3847/2515-5172/aab1f8
  * Fleming et al. (2020) Fleming, D. P., Barnes, R., Luger, R., & VanderPlas, J. T. 2020, ApJ, 891, 155, doi: 10.3847/1538-4357/ab77ad
  * Foley & Driscoll (2016) Foley, B. J., & Driscoll, P. E. 2016, Geochemistry, Geophysics, Geosystems, 17, 1885, doi: 10.1002/2015GC006210
  * Foley & Smye (2018) Foley, B. J., & Smye, A. J. 2018, Astrobiology, 18, 873, doi: 10.1089/ast.2017.1695
  * France et al. (2020) France, K., Duvvuri, G., Egan, H., et al. 2020, AJ, 160, 237, doi: 10.3847/1538-3881/abb465
  * Fuentes et al. (2019) Fuentes, J. J., Crowley, J. W., Dasgupta, R., & Mitrovica, J. X. 2019, Earth and Planetary Science Letters, 511, 154, doi: 10.1016/j.epsl.2019.01.020
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109, doi: 10.3847/1538-3881/aa80eb
  * Gaillard & Scaillet (2014) Gaillard, F., & Scaillet, B. 2014, Earth and Planetary Science Letters, 403, 307, doi: 10.1016/j.epsl.2014.07.009
  * Garraffo et al. (2017) Garraffo, C., Drake, J. J., Cohen, O., Alvarado-Gómez, J. D., & Moschou, S. P. 2017, ApJ, 843, L33, doi: 10.3847/2041-8213/aa79ed
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456, doi: 10.1038/nature21360
  * Gordiets & Kulikov (1985) Gordiets, B. F., & Kulikov, Y. N. 1985, Advances in Space Research, 5, 113, doi: 10.1016/0273-1177(85)90278-9
  * Grasser et al. (2023) Grasser, N., Kislyakova, K., Scherf, M., Lammer, H., & Van Looveren, G. 2023, Earth and Planetary Science Letters, 623, 118442, doi: 10.1016/j.epsl.2023.118442
  * Greene et al. (2023) Greene, T. P., Bell, T. J., Ducrot, E., et al. 2023, Nature, 618, 39, doi: 10.1038/s41586-023-05951-7
  * Hamano et al. (2013) Hamano, K., Abe, Y., & Genda, H. 2013, Nature, 497, 607, doi: 10.1038/nature12163
  * Hier-Majumder & Hirschmann (2017) Hier-Majumder, S., & Hirschmann, M. M. 2017, Geochemistry, Geophysics, Geosystems, 18, 3078, doi: 10.1002/2017GC006937
  * Hirth & Kohlstedt (1996) Hirth, G., & Kohlstedt, D. L. 1996, Earth and Planetary Science Letters, 144, 93, doi: 10.1016/0012-821X(96)00154-9
  * Ih et al. (2023) Ih, J., Kempton, E. M. R., Whittaker, E. A., & Lessard, M. 2023, ApJ, 952, L4, doi: 10.3847/2041-8213/ace03b
  * Jakosky et al. (2018) Jakosky, B. M., Brain, D., Chaffin, M., et al. 2018, Icarus, 315, 146, doi: 10.1016/j.icarus.2018.05.030
  * Kaltenegger (2017) Kaltenegger, L. 2017, ARA&A, 55, 433, doi: 10.1146/annurev-astro-082214-122238
  * Kane et al. (2020) Kane, S. R., Roettenbacher, R. M., Unterborn, C. T., Foley, B. J., & Hill, M. L. 2020, , 1, 36, doi: 10.3847/PSJ/abaab5
  * Kasting (1988) Kasting, J. F. 1988, Icarus, 74, 472, doi: 10.1016/0019-1035(88)90116-9
  * Korenaga (2007) Korenaga, J. 2007, Journal of Geophysical Research (Solid Earth), 112, B05408, doi: 10.1029/2006JB004502
  * Kreidberg et al. (2019) Kreidberg, L., Koll, D. D. B., Morley, C., et al. 2019, Nature, 573, 87, doi: 10.1038/s41586-019-1497-4
  * Krissansen-Totton (2023) Krissansen-Totton, J. 2023, ApJ, 951, L39, doi: 10.3847/2041-8213/acdc26
  * Krissansen-Totton & Fortney (2022) Krissansen-Totton, J., & Fortney, J. J. 2022, ApJ, 933, 115, doi: 10.3847/1538-4357/ac69cb
  * Kubyshkina et al. (2018) Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018, ApJ, 866, L18, doi: 10.3847/2041-8213/aae586
  * Kulikov et al. (2006) Kulikov, Y. N., Lammer, H., Lichtenegger, H. I. M., et al. 2006, Planet. Space Sci., 54, 1425, doi: 10.1016/j.pss.2006.04.021
  * Lammer et al. (2008) Lammer, H., Kasting, J. F., Chassefière, E., et al. 2008, Space Sci. Rev., 139, 399, doi: 10.1007/s11214-008-9413-5
  * Lammer et al. (2007) Lammer, H., Lichtenegger, H. I. M., Kulikov, Y. N., et al. 2007, Astrobiology, 7, 185, doi: 10.1089/ast.2006.0128
  * Lammer et al. (2018) Lammer, H., Zerkle, A. L., Gebauer, S., et al. 2018, A&A Rev., 26, 2, doi: 10.1007/s00159-018-0108-y
  * Lincowski et al. (2023) Lincowski, A. P., Meadows, V. S., Zieba, S., et al. 2023, ApJ, 955, L7, doi: 10.3847/2041-8213/acee02
  * Luger et al. (2017) Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nature Astronomy, 1, 0129, doi: 10.1038/s41550-017-0129
  * Lustig-Yaeger et al. (2019) Lustig-Yaeger, J., Meadows, V. S., & Lincowski, A. P. 2019, AJ, 158, 27, doi: 10.3847/1538-3881/ab21e0
  * Lustig-Yaeger et al. (2023) Lustig-Yaeger, J., Fu, G., May, E. M., et al. 2023, Nature Astronomy, doi: 10.1038/s41550-023-02064-z
  * Ma et al. (2004) Ma, Y., Nagy, A. F., Sokolov, I. V., & Hansen, K. C. 2004, Journal of Geophysical Research (Space Physics), 109, A07211, doi: 10.1029/2003JA010367
  * Ma et al. (2013) Ma, Y. J., Nagy, A. F., Russell, C. T., et al. 2013, Journal of Geophysical Research (Space Physics), 118, 321, doi: 10.1029/2012JA018265
  * Magaudda et al. (2020) Magaudda, E., Stelzer, B., Covey, K. R., et al. 2020, A&A, 638, A20, doi: 10.1051/0004-6361/201937408
  * Moran et al. (2023) Moran, S. E., Stevenson, K. B., Sing, D. K., et al. 2023, ApJ, 948, L11, doi: 10.3847/2041-8213/accb9c
  * Morin et al. (2010) Morin, J., Donati, J. F., Petit, P., et al. 2010, MNRAS, 407, 2269, doi: 10.1111/j.1365-2966.2010.17101.x
  * Morley et al. (2017) Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, ApJ, 850, 121, doi: 10.3847/1538-4357/aa927b
  * Odert et al. (2018) Odert, P., Lammer, H., Erkaev, N. V., et al. 2018, Icarus, 307, 327, doi: 10.1016/j.icarus.2017.10.031
  * Raymond et al. (2004) Raymond, S. N., Quinn, T., & Lunine, J. I. 2004, Icarus, 168, 1, doi: 10.1016/j.icarus.2003.11.019
  * Raymond et al. (2022) Raymond, S. N., Izidoro, A., Bolmont, E., et al. 2022, Nature Astronomy, 6, 80, doi: 10.1038/s41550-021-01518-6
  * Regenauer-Lieb et al. (2001) Regenauer-Lieb, K., Yuen, D. A., & Branlund, J. 2001, Science, 294, 578, doi: 10.1126/science.1063891
  * Reiners & Basri (2010) Reiners, A., & Basri, G. 2010, ApJ, 710, 924, doi: 10.1088/0004-637X/710/2/924
  * Salz et al. (2015) Salz, M., Schneider, P. C., Czesla, S., & Schmitt, J. H. M. M. 2015, A&A, 576, A42, doi: 10.1051/0004-6361/201425243
  * Sleep & Zahnle (2001) Sleep, N. H., & Zahnle, K. 2001, J. Geophys. Res., 106, 1373, doi: 10.1029/2000JE001247
  * Sreejith et al. (2020) Sreejith, A. G., Fossati, L., Youngblood, A., France, K., & Ambily, S. 2020, A&A, 644, A67, doi: 10.1051/0004-6361/202039167
  * Teixeira & Foley (2023) Teixeira, K., & Foley, B. 2023, TRAPPIST-1c Code, 1.0.0, Zenodo, doi: 10.5281/zenodo.10052104.  https://doi.org/10.5281/zenodo.10052104
  * Tian (2009) Tian, F. 2009, ApJ, 703, 905, doi: 10.1088/0004-637X/703/1/905
  * Tikoo & Elkins-Tanton (2017) Tikoo, S. M., & Elkins-Tanton, L. T. 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20150394, doi: 10.1098/rsta.2015.0394
  * Unterborn et al. (2022) Unterborn, C. T., Foley, B. J., Desch, S. J., et al. 2022, ApJ, 930, L6, doi: 10.3847/2041-8213/ac6596
  * Unterborn et al. (2015) Unterborn, C. T., Johnson, J. A., & Panero, W. R. 2015, ApJ, 806, 139, doi: 10.1088/0004-637X/806/1/139
  * van der Holst et al. (2014) van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81, doi: 10.1088/0004-637X/782/2/81
  * West et al. (2006) West, A. A., Bochanski, J. J., Hawley, S. L., et al. 2006, AJ, 132, 2507, doi: 10.1086/508652
  * Wood et al. (2021) Wood, B. E., Müller, H.-R., Redfield, S., et al. 2021, ApJ, 915, 37, doi: 10.3847/1538-4357/abfda5
  * Young et al. (2023) Young, E. D., Shahar, A., & Schlichting, H. E. 2023, Nature, 616, 306, doi: 10.1038/s41586-023-05823-0
  * Zahnle & Catling (2017) Zahnle, K. J., & Catling, D. C. 2017, ApJ, 843, 122, doi: 10.3847/1538-4357/aa7846
  * Zieba et al. (2022) Zieba, S., Zilinskas, M., Kreidberg, L., et al. 2022, A&A, 664, A79, doi: 10.1051/0004-6361/202142912
  * Zieba et al. (2023) Zieba, S., Kreidberg, L., Ducrot, E., et al. 2023, Nature, 620, 746, doi: 10.1038/s41586-023-06232-z