# Atmospheric Escape From Three Terrestrial Planets in the L 98-59 System

Emeline F. Fromont Department of Astronomy, University of Maryland, College
Park, MD 20742, USA  NASA Goddard Space Flight Center, Greenbelt, MD 20771,
USA  John P. Ahlers NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
University of Maryland, Baltimore County, Baltimore, MD 21250, USA  Laura N.
R. do Amaral Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de
México, Cto. Exterior S/N, C.U., Coyoacán, 04510 Ciudad de Mexico, CDMX  Rory
Barnes Department of Astronomy, University of Washington, Seattle, WA 98105,
USA  Emily A. Gilbert Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA  Elisa V. Quintana
NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA  Sarah Peacock
University of Maryland, Baltimore County, Baltimore, MD 21250, USA  NASA
Goddard Space Flight Center, Greenbelt, MD 20771, USA  Thomas Barclay NASA
Goddard Space Flight Center, Greenbelt, MD 20771, USA  University of Maryland,
Baltimore County, Baltimore, MD 21250, USA  Allison Youngblood NASA Goddard
Space Flight Center, Greenbelt, MD 20771, USA

###### Abstract

A critically important process affecting the climate evolution and potential
habitability of an exoplanet is atmospheric escape, in which high-energy
radiation from a star drives the escape of hydrogen atoms and other light
elements from a planet’s atmosphere. L 98-59 is a benchmark system for
studying such atmospheric processes, with three transiting terrestrial-size
planets receiving Venus-like instellations (4–25 S⊕) from their M3 host star.
We use the VPLanet model (Barnes et al., 2020) to simulate the evolution of
the L 98-59 system and the atmospheric escape of its inner three small
planets, given different assumed initial water quantities. We find that,
regardless of their initial water content, all three planets accumulate
significant quantities of oxygen due to efficient water photolysis and
hydrogen loss. All three planets also receive enough XUV flux to drive rapid
water loss, which considerably affects their developing climates and
atmospheres. Even in scenarios of low initial water content, our results
suggest that the James Webb Space Telescope (JWST) will retained oxygen on the
L 98-59 planets in its future scheduled observations, with planets b and c
being the most likely targets to possess an extended atmosphere. Our results
constrain the atmospheric evolution of these small rocky planets, and they
provide context for current and future observations of the L 98-59 system to
generalize our understanding of multi-terrestrial planet systems.

##  1 Introduction

To date, both ground-based radial velocity (RV) surveys and space-based
transit surveys have found small (R$<$1.6 R⊕; Cloutier & Menou, 2020) planets
at higher occurrence rates around M dwarfs than for hotter stars (Bonfils et
al., 2013; Dressing & Charbonneau, 2013, 2015; Hardegree-Ullman et al., 2019;
Ment & Charbonneau, 2023). As of June 2023, NASA’s Exoplanet Archive lists 181
confirmed and 116 candidate small planets around M dwarf stars. Additionally,
many of these planets orbit nearby, bright stars, making them ideal candidates
for atmospheric characterization by JWST (Gardner et al., 2006).

While hundreds of multi-planet systems have been discovered, only a limited
number are amenable to follow-up observations. Most multi-planet systems have
host stars that are too faint or have predicted signal amplitudes that are too
small to be useful for follow-up mass and atmospheric measurements.

Situated at the border of JWST’s continuous viewing zone (Demangeon et al.,
2021), the L 98-59 system is an excellent benchmark target for JWST
observations. As shown in Figure 1, the L 98-59 planets occupy an enticing
region of parameter space for follow-up study. The system is composed of four
confirmed planets, L 98-59 b, L 98-59 c, L 98-59 d, and L 98-59 e, which orbit
around their bright (K = 7.1 mag) and nearby (10.6 pc; Kostov et al., 2019) M3
dwarf host star. The inner three small exoplanets, L 98-59 b, L 98-59 c, and L
98-59 d reside interior to the system’s habitable zone (HZ). Our work focuses
only on these three planets, as L 98-59 e likely does not transit (Demangeon
et al., 2021).

Figure 1:

There is considerable value in modeling the evolution of the L 98-59 system,
particularly because of a unique combination of characteristics including the
presence of multiple terrestrial planets, their likely exposure to flares, and
the system’s overall ease of observation. These characteristics will allow us
to better understand the interactions of planets around M dwarfs, climate
development on Earth-analogs, and potential habitability conditions that arise
from such an environment.

The L 98-59 system presents an excellent opportunity to study the evolution
and atmospheric characterization of multiple, small planets that formed in the
same stellar environment (Greene et al., 2016; Morley et al., 2017; Demangeon
et al., 2021). L 98-59 is the brightest and nearest M dwarf star system that
has at least two measured planetary masses and radii, as seen in Figure 2.
While L 98-59 appeared to be a quiet star when it was first observed by TESS
(Ricker et al., 2015; Kostov et al., 2019), subsequent observations revealed
stellar activity in the form of white-light flares. M dwarf stars produce
frequent flares across the electromagnetic spectrum (e.g. Muirhead et al.,
2018; do Amaral et al., 2022) and remain active for long portions of their
lifetimes (West et al., 2008, 2015). This activity likely affects planet
atmospheres and must be accounted for in our atmospheric evolution models.
Studying the atmospheres of planets orbiting active stars such as L 98-59 will
provide us with a greater understanding of the effects of high-energy, or
short-wavelength, radiation on the atmospheric retention and evolution of
terrestrial planets over time.

Figure 2: The distance from the Sun and V-band magnitude of Figure 1. Most of
these systems are very faint and/or very distant, making them difficult to
observe as atmospheric characterization targets. L 98-59 is both nearby and
bright, making it a prime target for JWST observations.

Atmospheric escape is an important process to consider as part of the
evolution and potential habitability of exoplanets. Strong X-ray and extreme
ultraviolet (XUV; 1-1000 Å, Ribas et al. (2005)) radiation from the host star
drives heating and ionization of the upper atmospheres of highly irradiated
planets, leading to the escape of gases (Luger et al., 2015; Koskinen et al.,
2022). XUV-driven escape is believed to have strongly affected the atmospheric
evolution of solar system planets, such as Venus, Earth, and Mars (Lammer et
al., 2008), and is considered to be a likely culprit for sculpting the
observed exoplanet population (Owen et al., 2020). Indeed, planets orbiting
close to their host stars, like the L 98-59 planets, are particularly
vulnerable to atmospheric escape. Transit observations by the Hubble Space
Telescope (HST) in the far-UV (FUV) spectrum confirm that many planets close
to their host stars lose mass to hydrodynamic escape (e.g., HD 209458b (Vidal-
Madjar et al., 2003), HD 189733b (Lecavelier des Etangs et al., 2012), GJ 436b
(Ehrenreich et al., 2015), GJ 3470b (Bourrier et al., 2018)), which can have
important consequences not only for planets’ atmospheric evolution, but for
their structure and composition as well (Luger et al., 2015).

If the inner three planets of the system contained any volatiles in their
initial composition, they should experience atmospheric escape because of
their proximity to their star. To study the impact of XUV-driven escape on the
L 98-59 planets, we use the modeling tool VPLanet (Barnes et al., 2020) to
evolve the system over billions of years. By studying the atmospheric loss of
this system, we hope to better understand the effects of an M dwarf’s stellar
environment on the development of terrestrial-size planets, which may impact
their promise as potential targets for biosignature detection.

In the following sections, we describe our findings on the atmospheric
evolution of the three inner planets in the L 98-59 system. We start by
discussing the observational data and system properties of the L 98-59 system
in Section 2. We then describe our planetary simulation model in Section 3. We
interpret our model’s results in Section 4 and discuss their implications for
future studies of the system in Section 5. Finally, we conclude with an
overview of our results and their importance in the context of a new era of
JWST observations in Section 6.

##  2 System Observations and Properties

###  2.1 Observations

TESS observed L 98-59 (TIC 307210830, TOI 175) . Kostov et al. (2019)
initially reported the system’s discovery, in which they confirmed the
presence of the inner three terrestrial planets (b, c, d). Since then,
Cloutier et al. (2019) and Demangeon et al. (2021) have provided planetary
masses and eccentricities (Section 2.3). Although Demangeon et al. (2021) has
recently confirmed the fourth planet, L 98-59 e, (likely a rocky planet or a
water world), we do not consider this planet in our model and analysis because
it does not transit. We set our stellar and planetary parameters using the
values displayed in Table 1.

owing to L 98-59’s relative brightness for a host of small planets, it is an
excellent transmission spectroscopy target with HST and JWST . No atmospheric
signal was conclusively detected by HST/WFC3 in five transits of planet b and
one transit each for planets c and d. The data from planet b did not
demonstrate any evidence for an atmosphere (Damiano et al., 2022), and the
data collected for planet c was inconclusive, showing modulation in the
transmission spectra that could be interpreted as evidence of an atmosphere,
but at low significance (Barclay et al., 2023). Moreover, there were
indications that the signal observed could be a result of contaminating signal
from inhomogeneities in the stellar photosphere (Barclay et al., 2021, 2023).
Two upcoming HST/WFC3 transits of planet c may resolve the ambiguity in the
results for L 98-59 c. Through Guaranteed Time Observations, JWST will observe
L 98-59 c with NIRISS and L 98-59 d with NIRSpec.

Table 1: Stellar and planetary system parameters Parameter | Value adopted for this paper | Source
---|---|---
L 98-59 |  |
Mass (M☉) | $0.273\pm 0.030$ | Demangeon et al. (2021)
Radius (R☉) | $0.303\pm 0.0245$ | Demangeon et al. (2021)
Temperature (K) | $3415\pm 135$ | Demangeon et al. (2021)
XUV saturation fraction | $0.002\pm 0.001$ | Peacock et al. (2020)
XUV beta exponent | $0.882\pm 0.139$ | Peacock et al. (2020)
XUV saturation time (Myr) | 650 | Peacock et al. (2020)
Age (Gyr) | $>1$ | Demangeon et al. (2021)
Min flare energy (erg) | $1.73\pm 0.39\times 10^{31}$ | Barclay et al. (2023)
Max flare energy (erg) | $7.14\pm 0.88\times 10^{31}$ | Barclay et al. (2023)
L 98-59 b |  |
Mass (M⊕) | $0.4\pm 0.155$ | Demangeon et al. (2021)
Radius (R⊕) | $0.850\pm 0.054$ | Demangeon et al. (2021)
Density (g/cm3) | $3.6$ | Demangeon et al. (2021)
Flux (S⊕) | $24.7$ | Demangeon et al. (2021)
Period (days) | $2.2531136$ | Demangeon et al. (2021)
Semi-major axis (AU) | $0.0218$ | Demangeon et al. (2021)
Eccentricity | $0.103$ | Demangeon et al. (2021)
Inclination | 0.0 | Demangeon et al. (2021)
Argument of periastron (deg) | $192$ | Demangeon et al. (2021)
Transit time (BJDTDB \- 2457000) | $1366.17067$ | Demangeon et al. (2021)
L 98-59 c |  |
Mass (M⊕) | $2.22\pm 0.255$ | Demangeon et al. (2021)
Radius (R⊕) | $1.385\pm 0.085$ | Demangeon et al. (2021)
Density (g/cm3) | $4.57$ | Demangeon et al. (2021)
Flux (S⊕) | $12.8$ | Demangeon et al. (2021)
Period (days) | $3.6906777$ | Demangeon et al. (2021)
Semi-major axis (AU) | $0.0303$ | Demangeon et al. (2021)
Eccentricity | $0.103$ | Demangeon et al. (2021)
Inclination | 0.0 | Demangeon et al. (2021)
Argument of periastron (deg) | $261$ | Demangeon et al. (2021)
Transit time (BJDTDB \- 2457000) | $1367.27375$ | Demangeon et al. (2021)
L 98-59 d |  |
Mass (M⊕) | $1.94\pm 0.28$ | Demangeon et al. (2021)
Radius (R⊕) | $1.521\pm 0.1085$ | Demangeon et al. (2021)
Density (g/cm3) | $2.95$ | Demangeon et al. (2021)
Flux (S⊕) | $5.01$ | Demangeon et al. (2021)
Period (days) | $7.4507245$ | Demangeon et al. (2021)
Semi-major axis (AU) | $0.0484$ | Demangeon et al. (2021)
Eccentricity | $0.0740$ | Demangeon et al. (2021)
Inclination | 0.0 | Demangeon et al. (2021)
Argument of periastron (deg) | $180$ | Demangeon et al. (2021)
Transit time (BJDTDB \- 2457000) | $1362.73974$ | Demangeon et al. (2021)
All planets |  |
XUV absorption efficiency | $0.15-0.3$ | Luger & Barnes (2015)
XUV absorption radius/planet radius | 1.0 | Luger & Barnes (2015)

###  2.2 Host Star Properties

L 98-59 is a small, main-sequence M3 dwarf star with a luminosity about 100
times fainter than our sun (0.012 L☉; Cloutier et al., 2019; Kostov et al.,
2019; Demangeon et al., 2021). Due to its slow rotation, Kostov et al. (2019)
estimated L 98-59’s age to be $>$1 Gyr. Follow-up calculations by Demangeon et
al. (2021) agree with this result, stating the age of L 98-59 to be over 800
Myr. This star was also thought to be relatively quiet with little stellar
activity, based on lack of H$\alpha$ emission in optical spectra
(Pidhorodetska et al., 2021; Kostov et al., 2019). However, recent studies
have confirmed the detection of flares on L 98-59 (Stelzer et al., 2022;
Howard, 2022; Barclay et al., 2023).

###  2.3 Planetary System Properties

There are four planets in the L 98-59 system, two of which are . Initially,
Cloutier et al. (2019) were unable to constrain the mass of L 98-59 b, instead
defining an upper limit of 1.01 Earth masses (M⊕). Demangeon et al. (2021)
were later able to refine the planets’ masses to be $0.4^{+0.16}_{-0.15}$,
$2.22^{+0.26}_{-0.25}$, and $1.94\pm 0.28$ M⊕ for planets b, c, and d
respectively. As part of their RV analysis, Cloutier et al. (2019) also
derived densities for L 98-59 c and d, determining planet c to have a rock-
dominated composition and planet d to likely have a significant amount of
water or gaseous atmosphere due to its lower density. Following up on their
results, Demangeon et al. (2021) determined that the three inner planets have
small iron cores, with densities $3.6^{+1.4}_{-1.5}$ g/cm3 (b),
$4.57^{+0.77}_{-0.85}$ g/cm3 (c), $2.95^{+0.79}_{-0.51}$ g/cm3 (d). Their
analysis showed that planets b and c have very similar compositions, likely
rocky and with a small mass fraction of water and a low gas mass, if any. On
the other hand, their model favored a richer gas and water content for planet
d.

Quick et al. (2020) reached similar conclusions in their volcanic and
cryovolcanic analysis. Using estimates of the planets’ total internal heat,
which includes both tidal heating and radiogenic heating, they estimated that
L 98-59 b and c are likely to exhibit extreme volcanic activity at their
surfaces. They also suggested that L 98-59 d may be an ocean world with some
volcanic activity, based on the planet’s density and effective temperature.

The L 98-59 system resembles other multi-planet systems around M dwarfs, such
as TRAPPIST-1 and Kepler-186 (Gillon et al., 2017; Quintana et al., 2014), in
that M dwarfs frequently host compact, multi-planet systems (Kostov et al.,
2019; Muirhead et al., 2015). At 0.0218 AU (b), 0.0303 AU (c), and 0.0484 AU
(d), the three inner planets reside interior to both the system’s conservative
and optimistic HZ, as calculated from Kopparapu et al. (2013, 2014) using the
values from Table 1. These HZ limit values are listed in Table 2. The system
also follows the “peas in a pod” configuration initially described by Weiss et
al. (2018), where within a multi-planet system, planets tend to have similar
sizes and in systems with more than three planetary bodies, planets tend to be
evenly spaced in their orbits.

Table 2: Optimistic and Conservative Habitable Zones of L 98-59 HZ Limit | Flux | Distance
---|---|---
Conservative |  |
Inner HZ | 0.9298 S⊕ | 0.109 AU
Outer HZ | 0.244 S⊕ | 0.212 AU
Optimistic |  |
Inner HZ | 1.492 S⊕ | 0.086 AU
Outer HZ | 0.22 S⊕ | 0.224 AU

##  3 Numerical Methods and Models

###  3.1 System Evolution Model

To study the atmospheric development of the L 98-59 planets, we simulate the
evolution of the system over several billion years assuming an initial water
content between 1 and 100 terrestrial oceans (TO) and accounting for the host
star’s radiation and flares. We use the modeling software package VPLanet
(Virtual Planet Simulator; Barnes et al., 2020) to model the atmospheric
escape of the three innermost, terrestrial planets of the L 98-59 system.
VPLanet is an open-source software for simulating the evolution of a planetary
system, focusing specifically on habitability. By coupling various models, it
allows users to incorporate specific modules and features into their planetary
system, such as tidal evolution or a circumbinary system. We will discuss each
module we use in our model in Sections 3.1.1 (STELLAR), 3.1.2 (FLARE), and
3.1.3 (AtmEsc). We evolve the system from before the star’s main sequence (MS)
phase (5 Myr) until the upper limit for the star’s age at the present day (13
Gyrs), in order to visualize the full timescale of this system’s evolution

We make three important assumptions in the initial setup of our model. First,
we assume that the planets are in their currently-observed orbits at the start
of each simulation. We thus begin the simulations with the planets in their
present orbits around L 98-59, and assume that any planet migration in the
early stages of the system’s evolution has already taken place. This first
assumption ties into our second one, in which we begin the simulations before
the star’s main sequence phase, at a stellar age of 5 Myr. Pre-MS M dwarfs
emit higher amounts of flux in the XUV as compared to MS stars (Ramirez &
Kaltenegger, 2014; Luger & Barnes, 2015; Tian & Ida, 2015; do Amaral et al.,
2022). Excess XUV flux can cause more water loss and oxygen accumulation early
in a planet’s evolution, so it is important to include that early evolution
time in our simulations; by not doing so, we might severely underestimate the
desiccation of the L 98-59 planets. Third, we assume that all three planets
start with the same water mass at the beginning of each simulation. This
allows us to directly compare the effects of L 98-59’s stellar activity and
evolution on each planet in the system, and in turn, how water content changes
on each of them.

We investigate the evolution and outcome of three different atmospheric escape
scenarios by varying the initial water content of the L 98-59 planets. do
Amaral et al. (2022) found that the amount of escaping water is independent
from the starting quantity of water on a planet, and that the star mass is
inversely correlated with the amount of water lost. In other words, planets
tend to lose the same amount of water regardless of their initial budget, and
smaller stars will contribute most to this water loss due to their longer pre-
main-sequence (PMS) phase. Given we do not know the underlying initial water
budget for the L 98-59 planets, we vary the starting water budgets, in order
to gain a broader understanding of the water retention on each planet over
time.

We vary the water budgets, running sets of simulations with initial water
masses of 1, 10, and 100 TO following the procedures of Morbidelli et al.
(2000); Raymond et al. (2006); Chassefière et al. (2012); Luger & Barnes
(2015); do Amaral et al. (2022) in order to account for the possibility of the
planets forming in situ, or migrating in from further out in the disk where
more water may have been available at planet formation (Mordasini et al.,
2012; Tian & Ida, 2015; Unterborn et al., 2018; Pidhorodetska et al., 2021).
Our parameter variation method is further discussed in Section 3.1.4.

We describe the specific modules used in our model in the following
subsections. Using the VPLanet modules STELLAR and FLARE, we account for both
XUV and flare activity from the L 98-59 star in our model. We also use the
AtmEsc module to simulate the planets’ atmospheric and oceanic evolution over
time.

####  3.1.1 STELLAR module

As a star evolves and ages, the changes in its stellar activity affect the
formation, evolution, and potential habitability of its surrounding planets.
Early M stars, such as L 98-59, emit near-constant elevated levels of XUV flux
for the first few hundred million years before decreasing with time ($\sim
t^{-1}$) by two orders of magnitude towards field age (Peacock et al., 2020).
This decrease in emission is linked to the spin down of the star reducing the
dynamo production of its magnetic field (West et al., 2015).

The STELLAR module allows users to include the star’s changing parameters into
the system evolution model, including its XUV parameters and stellar evolution
model (Baraffe et al., 2015). Following the results from Peacock et al.
(2020), we adopt an XUV saturation timescale of 650 Myr for L 98-59. To
determine the XUV saturation fraction and $\beta$ exponent for L 98-59, we use
evolutionary models of 0.35 M☉ stars from Peacock et al.
(2020)111https://stdatu.stsci.edu/hlsp/hazmat in combination with The
Röntgensatellit (ROSAT) X-ray measurements of proxy stars from Shkolnik &
Barman (2014) to create full panchromatic spectra representative of L 98-59 at
ages between 10 Myr and 5 Gyr2225 Gyr is the oldest spectrum available from
Peacock et al. (2020). At X-ray wavelengths (1-100 Å), we adopt a single value
flux consistent with the median ROSAT measurement for each given age bin.
Since the XUV flux declines in a consistent manner beyond 650 Myr, calculating
the XUV saturation fraction and $\beta$ exponent with models between 650 Myr
and 5 Gyr (rather than 13 Gyr) yields a sufficiently long baseline to use.. We
use the Peacock et al. (2020) models for wavelengths $>$100 Å and the ROSAT
measurements for wavelengths $<$100 Å. We compute the XUV (1 – 1000 Å) and
bolometric luminosities from these spectra and fit a power law to the
decreasing LXUV/Lbol to determine the $\beta$ exponent (the coefficient in
this power-law fit). We do this with the median models for each age and the
inner quartile models to provide an uncertainty: $\beta$ = $0.882\pm 0.139$.
The XUV saturation fraction is LXUV/Lbol at the saturation time of 650 Myr,
LXUV/Lbol = $0.002\pm 0.001$.

####  3.1.2 FLARE module

The presence of flares in a planetary system can have a strong impact on the
habitability of a planet, particularly on the retention and detection of its
atmosphere. Flares can alter the chemical composition of a planet’s
atmosphere, and repeated exposure to flares can cause a planet to quickly lose
any atmosphere it might have accumulated (Chen et al., 2021; Tilley et al.,
2019; Louca et al., 2023). Until now, little work has been done to study the
impacts of flare exposure on a planet’s ocean retention, even though the most
commonly observed stars for small, rocky exoplanet characterization – main
sequence M stars, or M dwarfs – are known to emit such radiation (Billings,
2011; Shields et al., 2016; Fujii et al., 2018; do Amaral et al., 2022). Our
evolution model thus includes the effects of flares on the atmospheric escape
and ocean retention of the L 98-59 planets.

The FLARE module in VPLanet enables users to specify basic information on the
star’s flare distribution. We identify 8 flares in Sector 11 of TESS
observations using a modified version of bayesflare (Pitkin et al., 2014) as
described in Gilbert et al. (2022). Using the identified flare parameters
(peak time, amplitude, FWHM) as identified in the flare detection, we model
these flares using xoflares (Barclay & Gilbert, 2020). We integrate under
these models to determine the equivalent duration of each flare, then scale
these to absolute energies using the convolution of an M3 spectral model of L
98-59 and the TESS bandpass. From here, we determine the mean energy of the
star and then use the model of the individual flares from Barclay et al.
(2023) to determine the minimum and maximum flare energies, which are
$1.73^{+0.38}_{-0.41}\times 10^{31}$ erg and $7.14^{+0.88}_{-0.88}\times
10^{31}$ erg, respectively. These uncertainties are averaged such that we use
a minimum flare energy of $1.73\pm 0.39\times 10^{31}$ erg and a maximum flare
energy of $7.14\pm 0.88\times 10^{31}$ erg. We use the “TESSUV” option for the
“sFlareBandPass” parameter, which follows Günther et al. (2020) to calculate
which fraction of the bolometric energy for observed flares fall into the
U-band (E${}_{U}\approx$ 7.6%Ebol).

####  3.1.3 AtmEsc module

Atmospheric escape is the process undergone by a planet when gases from its
atmosphere are lost to space. There are two broad categories of atmospheric
escape defined by the gas loss mechanism: thermal escape, and non-thermal
escape. Thermal escape includes Jeans escape, in which the temperature of a
gas accelerates its particles to above the escape velocity, and hydrodynamic
escape, in which escaping gas molecules drag along other molecules, creating a
fluid-like escape behavior. Non-thermal escape includes processes such as
photochemical loss, ion loss, and ionospheric outflow (among others), and
generally involves more complex interactions in ions and plasma (see Gronoff
et al. (2020) for more detail on atmospheric escape processes).

Thermal escape is considered the dominant mechanism for highly irradiated
planetary atmospheres, so the VPLanet model includes energy-limited and
diffusion-limited escape for H/He and water vapor atmospheres. Energy-limited
escape occurs as a hydrodynamic wind “blows away” hydrogen in the atmosphere,
and is driven by a fixed fraction of the incoming XUV energy. However, if not
all of the oxygen escapes and some remains in the atmosphere, the hydrogen
will have to escape by diffusing through a background of oxygen in the
atmosphere, thus causing diffusion-limited escape (Luger & Barnes, 2015).

The AtmEsc module in VPLanet describes the evolution of a planet’s atmosphere
and ocean retention. We specify each planet’s starting water budget as 1, 10,
or 100 TO, and vary this quantity in a way that is further discussed in
Section 3.1.4. We follow the water loss model described in Luger & Barnes
(2015) to simulate the evolution of water on each planet and to set some of
our parameters: we set the XUV absorption radius to be equal to the planet
radius, since XUV radiation is absorbed in the uppermost layers of a planet’s
atmosphere, and we set the XUV absorption efficiency parameter to vary between
0.15 - 0.3. We select the model from Bolmont et al. (2017) in the VPLanet code
to determine the XUV absorption efficiency parameter for water vapor. In these
simulations, we assume minimal O2 absorption at the surface (Lincowski et al.,
2018), so we do not include O2 loss processes other than the ones built into
the VPLanet integrator. We also assume no significant H/He envelopes. We keep
the Jeans Time, in which the ballistic escape of individual atoms drives
atmospheric mass loss in the low temperature limit (Luger et al., 2015), to
the default value of 1 Gyr.

####  3.1.4 Parameter Variation

In order to consider the full range of scenarios for the system’s evolution,
we vary the following parameters within their uncertainties: star mass (Ms),
planet mass for all three planets (Mp), planet radius for all three planets
(Rp), XUV saturation fraction, XUV beta exponent, minimum flare energy, and
maximum flare energy. Using Monte Carlo sampling, we draw 1000 values from a
Gaussian distribution for each of these parameters. We also include each
planet’s XUV absorption efficiency in our sampling as a flat distribution from
0.15 to 0.3, which is considered the range appropriate for planets with
hydrogen-rich atmospheres around M dwarfs (Luger & Barnes, 2015). We then use
these varying parameters to generate 1000 unique simulations, spanning the
parameter space of the L 98-59 system. We will refer to these 1000 simulations
as a set of simulations. We fix all other parameters in the model, excluding
the semi-major axis which varies indirectly due to Kepler’s third law.

We perform a set of simulations for three initial water masses: 1, 10, and 100
TO on each planet. For each initial water budget, we simulate the instellation
and each planet’s climatic response using input parameters as described above
to build a distribution of water loss and atmospheric escape scenarios. In our
model, we consider “desiccated” to mean that a planet’s water content has
dropped Thus, we model 1000 total scenarios per starting water budget, with
which we determine the range of most likely scenarios within 3$\sigma$ for
water retention and oxygen accumulation on the L 98-59 planets. That is, given
the distribution of our simulation results, we find the water evolution
scenario at +3$\sigma$ to represent the “most” or “longest” water retention
and the water evolution scenario at -3$\sigma$ to represent the “least” or
“shortest” water retention, giving us a range of most likely water and oxygen
evolution behaviors for each planet.

Figure 3: 1 TO These plots show the range of outcomes for oxygen accumulation
(top row, green) and water retention (bottom row, blue) for the set of
simulations with 1 TO initial water content on planets b (left column), c
(middle column), and d (right column). The dark and light blue lines represent
the simulations corresponding to the +3$\sigma$ and -3$\sigma$ (respectively)
of the water desiccation times distribution for each planet at 1 TO. These
distributions can be found in Figure 6, top row. The dark and light green
lines correspond to the same simulations as the ones chosen in the water
behavior plots. Given 1 TO, water isn’t likely to survive past 1 Gyr for any
planet, which makes water observations more unlikely. However, the quantities
of oxygen produced are still significant, up to bars. Notice that for planet
b, simulations that retain water longer are the ones that accumulate the least
oxygen, contrary to planets c and d. These b simulations also show a quick
initial drop in water content, and then a slower loss over time than the
simulations that retain water for shorter times. This is due to a high XUV
flux for b, which increases radiative cooling and decreases escape efficiency.
Further detail is provided in Section 4. Figure 4: 10 TO These plots show the
range of outcomes for oxygen accumulation (top row, green) and water retention
(bottom row, blue) for the set of simulations with 10 TO initial water content
on each of the planets. The distributions for the water behavior can be found
in Figure 6, middle row. Given 10 TO, all three planets can accumulate
hundreds of bars of oxygen, and planet b is likely to retain water for longer
than planets c and d (past 1 Gyr) making water detections more likely for this
planet. However, depending on the system’s age, all three planets may still
have some water left. Details on planet b’s behavior, similar to Figure 3, are
given in Section 4. Figure 5: 100 TO These plots show the range of outcomes
for oxygen accumulation (top row, green) and water retention (bottom row,
blue) for the set of simulations with 100 TO initial water content on each of
the planets. The distributions for the water behavior can be found in Figure
6, bottom row. Given 100 TO, all three planets are likely to retain water past
1 Gyr, making observations of water relatively likely. In addition, thousands
of bars of oxygen may be produced for each planet in this scenario. Details on
planet b’s behavior, similar to Figure 3, are given in Section 4.

##  4 Results

The main results from our simulations, showing the water retention and oxygen
accumulation over time for the L 98-59 planets at different starting water
budgets, are summarized in Figures 3, 4, and 5. We find that all three planets
are experiencing strong atmospheric escape due to their proximity and exposure
to their star’s XUV and flare activity. All three planets are able to
accumulate significant quantities of oxygen in their atmospheres over time.
While all three planets display significant water loss, planet b shows the
potential for water retention when given a medium to high initial water
content, but planets c and d are unlikely to retain significant amounts of
water unless given a high initial water content.

As seen in the histograms from Figure 6, when given 100 TO of starting water,
of the simulations predict that planet b will not desiccate and will instead
retain some fraction of water, as opposed to planets c and d which still tend
to desiccate in this situation. At higher XUV fluxes, such as ones that planet
b experiences, the escape efficiency decreases due to increased radiative
cooling. This means that less hydrogen is escaping, but that which is escaping
leaves the atmosphere at a high velocity and drags along large amounts of
oxygen with it (Tian et al., 2008; Bolmont et al., 2017).

Planets c and d are both less likely than planet b to retain significant
amounts of water over long periods of time, unless they begin evolving with
high quantities of water. Indeed, most of our VPLanet simulations predict that
in cases of 1 TO and 10 TO of initial water, planets c and d completely
desiccate before 1 Gyr, as seen in Figures 3, 4, unless they start their
evolution with a significant amount of water. Given 100 TO, planet c is able
to retain water for close to Gyr (Figure 6, bottom middle), and planet d shows
a few cases of retaining water for more than 13 Gyr (Figure 6, bottom right).

The atmospheric escape resulting from the star’s intense XUV radiation causes
the planets to accumulate significant quantities of oxygen, in both low and
high initial water content cases. Given 1 TO, the planets can accumulate from
around bars up to more than bars of O2. In addition, given 100 TO the planets
can accumulate thousands of bars of O2. Unless the planets began their orbits
completely desiccated and stripped of any volatiles, there is a strong chance
of observing secondary atmospheres dominated by oxygen on all three of them.
We note that these simulations do not take into account potential volcanic
activity on the planets, which as discussed in Section 2.3, is a likely
scenario for planets b and c. This implies that the planets may also
accumulate other volatiles produced by volcanic outgassing, besides O2. Our
simulations also do not account for potential sinks, like magma oceans – this
is further discussed in Section 5.4.

Examining the water evolution displayed in Figures 3, 4, and 5, we find that
the water loss and oxygen accumulation over time consistently do not behave
the same way for planet b as they do for planets c and d. For planet b, the
simulations that seem to retain water for shorter periods of time have a
rather fast and steady water loss rate right from the beginning, which only
slows down when the planet reaches insignificant levels of water. However, the
simulations that seem to retain water for longer begin their evolution with an
instantaneous drop in water quantity, then continue with a much slower water
loss rate. These simulations also accumulate much less oxygen than simulations
with shorter water retention times. By contrast, planets c and d display a
more constant water loss rate in all simulations throughout the sets, and
their simulations that retain water the longest also accumulate the most
oxygen.

This opposing behavior in water loss and oxygen retention is related to how
the escape efficiency depends on the XUV flux, for the same reason why planet
b is more likely than the other planets to retain some water, given a high
initial water content. Planet b experiences a higher XUV flux than the other
planets, such that it has a lower escape efficiency and could then retain
water for longer. But with a high XUV flux, the hydrolyzed oxygen that does
accumulate in the atmosphere will be dragged away by quickly escaping
hydrogen. Planets c and d, with a higher escape efficiency, will experience a
steeper and more constant water loss. Exposure to a less intense XUV flux also
means that they can accumulate more oxygen, since it does not get carried away
by hydrogen as dramatically as for planet b.

Figure 6: These plots show the distribution of the desiccation times for the
set of simulations run at 1 TO (top row), 10 TO (middle row), and 100 TO
(bottom row) for planets b (left column, green), c (middle column, light
blue), and d (right column, dark blue). The vertical green line and the
vertical black lines in each plot show the log mean and $\pm 3\sigma$ of the
results distribution, respectively. In all three initial water cases, planet b
on average retains water the longest. In the bottom left histogram, the
fullest bin on the far right of the plot shows all cases of water survival for
planet b at 100 TO. Note that the cases of water survival only , as the
simulations all stopped at time 1.3 Gyr. A distribution of only the desiccated
for planet b can be found in Figure 7. Figure 7: This plot shows the
distribution of the desiccated cases for the set of simulations run at 100 TO
for planet b, as seen on the bottom left plot of Figure 6. The vertical green
line and the vertical black lines show the log mean and $\pm 3\sigma$ of the
results distribution, respectively. is interesting to consider that, depending
on the L 98-59 system’s age, planet b may still have some water present.

##  5 Discussion

###  5.1 L 98-59 planets are in a runaway greenhouse

The L 98-59 planets have mean semi-major axes values of 0.0218 AU (b), 0.0303
AU (c), and 0.0484 AU (d), and thus have instellations ranging from about 4–25
times the instellation that the Earth receives from the Sun (Figure 1). These
planets therefore fall firmly in the Venus zone (Kane et al., 2014; Ostberg et
al., 2023), the region around a star within which a planet like Earth with
water on its surface would likely have been forced into a runaway greenhouse.
Most Venus analogs have been found around faint stars for which it is
difficult to obtain follow-up measurements. However, the L 98-59 planets’
location along with the brightness of their host star create a convenient
combination which would allow us to learn more about the evolution and
atmospheric development of Venus-like planets (Pidhorodetska et al., 2021). L
98-59 c and d, if they have atmospheres, are prime JWST targets to observe and
study potential Venus analogs and further study the origin and backstory of
our sister planet .

As seen in Figures 3, 4, and 5, the L 98-59 planets are likely experiencing a
runaway greenhouse, which makes them an excellent target for JWST observations
of evolving atmospheres and unstable climates. In a runaway greenhouse state,
the planets’ surface temperatures exceed the critical point of water (647 K)
and their oceans start evaporating (Kasting, 1988; Kasting et al., 1993; Abe,
1993). Our simulation results show that all three planets are constantly
losing water and accumulating high quantities of oxygen from water photolysis,
because they are so small and so close to their star. Schindler & Kasting
(2000) modeled synthetic spectra of hypothetical Earth-like atmospheres, which
showed that “Venus-like” planets could conceivably accumulate large amounts of
oxygen in their atmospheres. Luger & Barnes (2015) confirmed this behavior in
their planetary evolution and water loss model, demonstrating that planets
around the habitable zone of M dwarf stars may develop atmospheres with
hundreds to thousands of bars of O2. Thus, if the three planets had started
their evolution with any non-negligible water budget, it is likely that they
would display signs of water loss and volatile accumulation in their
atmospheres. Any detection of volatiles on the L 98-59 planets would show an
impermanent state in the lifetime of the system, but would also offer real-
time snapshots of a runaway greenhouse state and its effects on the climate
stability and evolving atmospheric composition of a planet.

###  5.2 Potential future volatile detections

Depending on the system’s age, there is a real possibility of observing
potential water signatures on L 98-59 b, as it seems prone to longer water
retention than L 98-59 c and d. Additionally, our simulations predict a
significant accumulation of oxygen in planet b’s atmosphere, especially in the
case of 100 TO. Such quantities of oxygen and water vapor, resulting from
water loss, may create a haze in the planet’s atmosphere, thus obscuring
potential water signal detections. Indeed, a recent transmission spectrum of L
98-59 b, found to be relatively flat, predicts that planet b either has no
atmosphere or an atmosphere with high-altitude clouds or haze (Damiano et al.,
2022). The authors’ data rule out a cloud-free and H2-dominated atmosphere,
which is consistent with our initial assumptions about the planet. However,
this transmission spectrum does not favor a water vapor-dominated atmosphere
either, although the possibility remains if clouds are included in the model.

JWST observations of planet b in parallel with planets c and d would provide
essential information on the evolution of different planets around the same M
dwarf. Not only would it allow us to compare the development and loss of
atmospheres on three small planets in the same system, but the detection of
H2O on L 98-59 b would allow us to place constraints on evolutionary
scenarios, by looking at how oceans and atmospheres behave when subjected to
such intense irradiation. For example, we may determine whether these planets
formed with a large initial amount of water, giving us a way to probe into the
mechanism of volatile accumulation during the initial formation stage of
planets.

###  5.3 Assumptions regarding flares and water delivery

As an M dwarf, L 98-59 is more prone to stellar variability and flares, which
presents a challenge for the planets’ ability to retain any significant
atmosphere or ocean. XUV radiation emitted by M dwarfs heats the exosphere of
a planet’s atmosphere, causing the exobase to expand, and therefore
facilitating escape (Murray-Clay et al., 2009; France et al., 2020). Thus,
continuous exposure to repeated flaring can significantly deplete an
atmosphere (Tilley et al., 2019; France et al., 2020). It is likely that any
water on the L 98-59 planets would be significantly affected, if not
desiccated, by the star’s flares and intense XUV radiation, as water-rich
atmospheres may be severely depleted during a star’s PMS phase (Tian & Ida,
2015; Luger & Barnes, 2015; do Amaral et al., 2022).

In addition, flares can generate up to hundreds of bars of additional O2 due
to water photolysis (do Amaral et al., 2022), and planets in an extended
runaway greenhouse in the HZ of M dwarfs can also accumulate hundreds to
thousands of bars of O2 (Luger & Barnes, 2015). So while these atmospheric
escape processes might create detectable amounts of oxygen on the L 98-59
planets, we note that this would be a false biosignature and an unreliable
measure of the habitability of these planets (Wordsworth & Pierrehumbert,
2014; Tian et al., 2014; Domagal-Goldman et al., 2014; Luger & Barnes, 2015).
Although M dwarfs like L 98-59 may seem unfavorable hosts to habitable
planets, with the combination of many evolutionary factors, there is still a
chance for them to host volatile-rich worlds, and they remain important
targets to study terrestrial planets subject to frequent high levels of XUV
radiation.

The presence of water on the L 98-59 planets is highly dependent on their
composition, which is in turn dependent on their formation process and
location. There are three possible formation scenarios: 1) formation at their
present location, 2) formation further out in the disk and inward migration to
their current locations, and 3) causing the planets to shift their orbits to
their current locations . The L 98-59 planets orbit extremely close to their
star, so it is most likely that they have migrated inward during the early
stages of the system’s formation (Ogihara & Ida, 2009; Pu & Lai, 2019).
Simulations by Raymond et al. (2007) also show that it is difficult to form
planets larger than 0.3 M⊕ so close to a low-mass star. Given the predictions
of studies by Ogihara & Ida (2009) and Ciesla et al. (2015), we consider there
to be a significant likelihood for the L 98-59 planets to be composed of a
non-negligible proportion of volatiles, including water, by the time they have
fully formed and reached their current orbits.

###  5.4 Limitations of this work

While our simulations show a range of outcomes that fall relatively in line
with predictions from previous works, we note that they may not fully capture
the specific details of the L 98-59 system’s evolution. As established by
Barnes et al. (2020), the AtmEsc module of the VPLanet model is only an
approximate description of atmospheric escape, and does not include some
processes present in small, terrestrial planets, including the wavelength
dependence of upper atmosphere heating and its variation with composition and
atmospheric temperature structure, line cooling mechanisms, and other non-
thermal escape processes. It also does not include CO2, which is a very common
volatile and likely to be found in the atmospheres of planets. But while the
VPLanet model does not account for some known effects in the evolution of
terrestrial planets, we believe that the results obtained give us a good
estimate of the general state of the L 98-59 planets. More importantly, it
provides us with a strong baseline for comparison with future observations of
this system, or others like it. This general knowledge can 1) be a benchmark
for what to expect when observing similar systems, 2) help us to interpret
these future observations, and 3) refine our models and understandings of the
evolution of small, multi-planet systems around M dwarfs.

Our simulations assume that no magma ocean is present on the planets. However
like many planetary-sized bodies, the L 98-59 planets were likely molten right
after their formation, a state known as a magma ocean (Solomon, 1979;
Wetherill, 1990; Lammer et al., 2018). This would affect their initial
volatile budget, how much oxygen may accumulate, and how fast water is lost.
During a magma ocean state, oxygen not only accumulates in the atmosphere but
may also enter the melted surface, the latter which could outgas H2O as a
result (see Barth et al. (2021)). These potential composition differences
would affect the quantity of oxygen and water in the atmosphere, the quantity
of oxygen stored in the planet, and the quantity of water produced by the
magma ocean. Perhaps this would decrease the amount of oxygen in the
atmosphere, and allow any water present to be “protected” by the magma ocean
and retained longer. We note that future works could include this magma ocean
state for a more complete picture of the L 98-59 planets’ evolution.

Finally, in their analysis of a recent planet b transmission spectrum, Damiano
et al. (2022) describe other potential models (besides cloud-free,
H2-dominated, and water-dominated) that include volatiles such as HCN, CO2,
CO, N2, and CH4, which our own simulations do not take into account. So while
our results point to the possibility of a water vapor and oxygen-dominated
atmosphere on L 98-59 b, our predictions may be limited because our model
assumes a high initial water content on the planet to achieve this result and
does not account for a varied set of volatiles.

##  6 Conclusion

In this paper, we model the evolution of the L 98-59 system to study the
effects of atmospheric escape on the potential volatiles present on the inner
three planets. Consisting of three small, transiting, and likely rocky planets
(and a fourth non-transiting planet) orbiting their M-dwarf star, this system
is an ideal candidate to better characterize the atmospheric evolution of
multiple terrestrial planets within the same stellar environment.

We run simulations with a parameter sweep of the system’s characteristics,
giving the planets an initial water budget of either 1, 10, or 100 TO each and
evolving them over several billion years. We gather the following results:

  1. 1.

strong atmospheric escape due to their proximity to L 98-59 and the exposure
to their star’s XUV and flare activity.

  2. 2.

The smallest and closest planet, b, may counter-intuitively retain water
longer than c and d, because its exposure to a higher XUV flux leads to
increased radiative cooling, and thus a lower escape efficiency.

  3. 3.

Given enough starting water (100 TO) and depending on the system’s age, all
three planets may still have water leftover today, with the highest chances of
detection being for planet b.

  4. 4.

All three planets can accumulate significant quantities of oxygen created by
water photolysis, from bars up to thousands of bars.

Multi-planet systems can provide us with important insight on planet formation
and evolution, orbital dynamics, and planetary architectures. Moreover,
planets around bright, nearby stars such as L 98-59 are ideal targets for
atmospheric characterization with emission and transmission spectroscopy
(Kostov et al., 2019). This system serves as an excellent case to study the
development of multiple terrestrial planets and their atmospheres.

Further, the presence of an atmosphere on any of the L 98-59 planets is highly
dependent on the way the planets evolved in the presence of their star’s
stellar activity, and specifically, the initial composition with which they
began evolving. Therefore, our work studying the effects of flares and strong
XUV radiation on volatile loss or accumulation for these planets may help us
to better constrain certain characteristics of the system, such as stellar
age, initial planetary composition, and planetary formation processes, in
future observations with JWST.

## Acknowledgements

The material is based upon work supported by NASA under award number
80GSFC21M0002. This research was carried out in part at the Jet Propulsion
Laboratory, California Institute of Technology, under a contract with the
National Aeronautics and Space Administration (80NM0018D0004). RB acknowledges
support from NASA grant 80NSSC23K0261. L.N.R.A. acknowledge the support of
UNAM DGAPA PAPIIT project IN110420. L.N.R.A. thanks CONACYT’s graduate
scholarship program for its support. E.F.F thanks Thaddeus D. Komacek for
insightful discussions regarding simulation results.

Software: VPLanet (Barnes et al., 2020), mercury (Chambers & Murison, 2000),
jupyter (Kluyver et al., 2016), NumPy (Harris et al., 2020), SciPy (Virtanen
et al., 2020), Matplotlib (Hunter, 2007), Pandas (pandas development team,
2020), SymPy (Meurer et al., 2017), seaborn (Waskom, 2021), fitter (Cokelaer
et al., 2023).

## References

  * Abe (1993) Abe, Y. 1993, Lithos, 30, 223, doi: 10.1016/0024-4937(93)90037-D
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42, doi: 10.1051/0004-6361/201425481
  * Barclay & Gilbert (2020) Barclay, T., & Gilbert, E. 2020, mrtommyb/xoflares, 0.2.1, Zenodo, doi: 10.5281/zenodo.4156285
  * Barclay et al. (2021) Barclay, T., Kostov, V. B., Colón, K. D., et al. 2021, AJ, 162, 300, doi: 10.3847/1538-3881/ac2824
  * Barclay et al. (2023) Barclay, T., Sheppard, K. B., Latouf, N., et al. 2023, arXiv e-prints, arXiv:2301.10866, doi: 10.48550/arXiv.2301.10866
  * Barnes et al. (2020) Barnes, R., Luger, R., Deitrick, R., et al. 2020, PASP, 132, 024502, doi: 10.1088/1538-3873/ab3ce8
  * Barth et al. (2021) Barth, P., Carone, L., Barnes, R., et al. 2021, Astrobiology, 21, 1325, doi: 10.1089/ast.2020.2277
  * Billings (2011) Billings, L. 2011, Nature, 470, 27, doi: 10.1038/470027a
  * Bolmont et al. (2017) Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728, doi: 10.1093/mnras/stw2578
  * Bonfils et al. (2013) Bonfils, X., Delfosse, X., Udry, S., et al. 2013, A&A, 549, A109, doi: 10.1051/0004-6361/201014704
  * Bourrier et al. (2018) Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., et al. 2018, A&A, 620, A147, doi: 10.1051/0004-6361/201833675
  * Burke et al. (2020) Burke, C. J., Levine, A., Fausnaugh, M., et al. 2020, TESS-Point: High precision TESS pointing tool, Astrophysics Source Code Library.  http://ascl.net/2003.001
  * Chambers & Murison (2000) Chambers, J. E., & Murison, M. A. 2000, AJ, 119, 425, doi: 10.1086/301161
  * Chassefière et al. (2012) Chassefière, E., Wieler, R., Marty, B., & Leblanc, F. 2012, Planet. Space Sci., 63, 15, doi: 10.1016/j.pss.2011.04.007
  * Chen et al. (2021) Chen, H., Zhan, Z., Youngblood, A., et al. 2021, Nature Astronomy, 5, 298, doi: 10.1038/s41550-020-01264-1
  * Chiang & Laughlin (2013) Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444, doi: 10.1093/mnras/stt424
  * Ciesla et al. (2015) Ciesla, F. J., Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 804, 9, doi: 10.1088/0004-637X/804/1/9
  * Cloutier & Menou (2020) Cloutier, R., & Menou, K. 2020, AJ, 159, 211, doi: 10.3847/1538-3881/ab8237
  * Cloutier et al. (2019) Cloutier, R., Astudillo-Defru, N., Bonfils, X., et al. 2019, A&A, 629, A111, doi: 10.1051/0004-6361/201935957
  * Cokelaer et al. (2023) Cokelaer, T., Kravchenko, A., lahdjirayhan, et al. 2023, cokelaer/fitter: v1.5.2, v1.5.2, Zenodo, doi: 10.5281/zenodo.7497983
  * Damiano et al. (2022) Damiano, M., Hu, R., Barclay, T., et al. 2022, AJ, 164, 225, doi: 10.3847/1538-3881/ac9472
  * Demangeon et al. (2021) Demangeon, O. D. S., Zapatero Osorio, M. R., Alibert, Y., et al. 2021, A&A, 653, A41, doi: 10.1051/0004-6361/202140728
  * 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
  * Domagal-Goldman et al. (2014) Domagal-Goldman, S. D., Segura, A., Claire, M. W., Robinson, T. D., & Meadows, V. S. 2014, ApJ, 792, 90, doi: 10.1088/0004-637X/792/2/90
  * Dressing & Charbonneau (2013) Dressing, C. D., & Charbonneau, D. 2013, ApJ, 767, 95, doi: 10.1088/0004-637X/767/1/95
  * Dressing & Charbonneau (2015) —. 2015, ApJ, 807, 45, doi: 10.1088/0004-637X/807/1/45
  * Ehrenreich et al. (2015) Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459, doi: 10.1038/nature14501
  * Engle & Guinan (2023) Engle, S. G., & Guinan, E. F. 2023, ApJ, 954, L50, doi: 10.3847/2041-8213/acf472
  * France et al. (2020) France, K., Duvvuri, G., Egan, H., et al. 2020, AJ, 160, 237, doi: 10.3847/1538-3881/abb465
  * Fujii et al. (2018) Fujii, Y., Angerhausen, D., Deitrick, R., et al. 2018, Astrobiology, 18, 739, doi: 10.1089/ast.2017.1733
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485, doi: 10.1007/s11214-006-8315-7
  * Gilbert et al. (2022) Gilbert, E. A., Barclay, T., Quintana, E. V., et al. 2022, AJ, 163, 147, doi: 10.3847/1538-3881/ac23ca
  * Gillon et al. (2017) Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456, doi: 10.1038/nature21360
  * Greene et al. (2016) Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17, doi: 10.3847/0004-637X/817/1/17
  * Gronoff et al. (2020) Gronoff, G., Arras, P., Baraka, S., et al. 2020, Journal of Geophysical Research (Space Physics), 125, e27639, doi: 10.1029/2019JA027639
  * Günther et al. (2020) Günther, M. N., Zhan, Z., Seager, S., et al. 2020, AJ, 159, 60, doi: 10.3847/1538-3881/ab5d3a
  * Hardegree-Ullman et al. (2019) Hardegree-Ullman, K. K., Cushing, M. C., Muirhead, P. S., & Christiansen, J. L. 2019, AJ, 158, 75, doi: 10.3847/1538-3881/ab21d2
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hord et al. (2023) Hord, B. J., Kempton, E. M. R., Mikal-Evans, T., et al. 2023, arXiv e-prints, arXiv:2308.09617, doi: 10.48550/arXiv.2308.09617
  * Howard (2022) Howard, W. S. 2022, MNRAS, 512, L60, doi: 10.1093/mnrasl/slac024
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Kane et al. (2014) Kane, S. R., Kopparapu, R. K., & Domagal-Goldman, S. D. 2014, ApJ, 794, L5, doi: 10.1088/2041-8205/794/1/L5
  * Kane et al. (2019) Kane, S. R., Arney, G., Crisp, D., et al. 2019, Journal of Geophysical Research (Planets), 124, 2015, doi: 10.1029/2019JE005939
  * Kasting (1988) Kasting, J. F. 1988, Icarus, 74, 472, doi: 10.1016/0019-1035(88)90116-9
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108, doi: 10.1006/icar.1993.1010
  * Kempton et al. (2018) Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401, doi: 10.1088/1538-3873/aadf6f
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Kopparapu et al. (2014) Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J., et al. 2014, ApJ, 787, L29, doi: 10.1088/2041-8205/787/2/L29
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131, doi: 10.1088/0004-637X/765/2/131
  * Koskinen et al. (2022) Koskinen, T. T., Lavvas, P., Huang, C., et al. 2022, ApJ, 929, 52, doi: 10.3847/1538-4357/ac4f45
  * Kostov et al. (2019) Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32, doi: 10.3847/1538-3881/ab2459
  * 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. (2018) Lammer, H., Zerkle, A. L., Gebauer, S., et al. 2018, A&A Rev., 26, 2, doi: 10.1007/s00159-018-0108-y
  * Lecavelier des Etangs et al. (2012) Lecavelier des Etangs, A., Bourrier, V., Wheatley, P. J., et al. 2012, A&A, 543, L4, doi: 10.1051/0004-6361/201219363
  * Lincowski et al. (2018) Lincowski, A. P., Meadows, V. S., Crisp, D., et al. 2018, ApJ, 867, 76, doi: 10.3847/1538-4357/aae36a
  * Louca et al. (2023) Louca, A. J., Miguel, Y., Tsai, S.-M., et al. 2023, MNRAS, 521, 3333, doi: 10.1093/mnras/stac1220
  * Luger & Barnes (2015) Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119, doi: 10.1089/ast.2014.1231
  * Luger et al. (2015) Luger, R., Barnes, R., Lopez, E., et al. 2015, Astrobiology, 15, 57, doi: 10.1089/ast.2014.1215
  * Ment & Charbonneau (2023) Ment, K., & Charbonneau, D. 2023, arXiv e-prints, arXiv:2302.04242, doi: 10.48550/arXiv.2302.04242
  * Meurer et al. (2017) Meurer, A., Smith, C. P., Paprocki, M., et al. 2017, PeerJ Computer Science, 3, e103, doi: 10.7717/peerj-cs.103
  * Morbidelli et al. (2000) Morbidelli, A., Chambers, J., Lunine, J. I., et al. 2000, , 35, 1309, doi: 10.1111/j.1945-5100.2000.tb01518.x
  * Mordasini et al. (2012) Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A&A, 547, A112, doi: 10.1051/0004-6361/201118464
  * 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
  * Muirhead et al. (2018) Muirhead, P. S., Dressing, C. D., Mann, A. W., et al. 2018, AJ, 155, 180, doi: 10.3847/1538-3881/aab710
  * Muirhead et al. (2015) Muirhead, P. S., Mann, A. W., Vanderburg, A., et al. 2015, ApJ, 801, 18, doi: 10.1088/0004-637X/801/1/18
  * Murray-Clay et al. (2009) Murray-Clay, R. A., Chiang, E. I., & Murray, N. 2009, ApJ, 693, 23, doi: 10.1088/0004-637X/693/1/23
  * Ogihara & Ida (2009) Ogihara, M., & Ida, S. 2009, ApJ, 699, 824, doi: 10.1088/0004-637X/699/1/824
  * Ostberg et al. (2023) Ostberg, C., Kane, S. R., Li, Z., et al. 2023, arXiv e-prints, arXiv:2302.03055, doi: 10.48550/arXiv.2302.03055
  * Owen et al. (2020) Owen, J. E., Shaikhislamov, I. F., Lammer, H., Fossati, L., & Khodachenko, M. L. 2020, Space Sci. Rev., 216, 129, doi: 10.1007/s11214-020-00756-w
  * pandas development team (2020) pandas development team, T. 2020, pandas-dev/pandas: Pandas, latest, Zenodo, doi: 10.5281/zenodo.3509134
  * Peacock et al. (2020) Peacock, S., Barman, T., Shkolnik, E. L., et al. 2020, ApJ, 895, 5, doi: 10.3847/1538-4357/ab893a
  * Pidhorodetska et al. (2021) Pidhorodetska, D., Moran, S. E., Schwieterman, E. W., et al. 2021, AJ, 162, 169, doi: 10.3847/1538-3881/ac1171
  * Pitkin et al. (2014) Pitkin, M., Williams, D., Fletcher, L., & Grant, S. D. T. 2014, MNRAS, 445, 2268, doi: 10.1093/mnras/stu1889
  * Pu & Lai (2019) Pu, B., & Lai, D. 2019, MNRAS, 488, 3568, doi: 10.1093/mnras/stz1817
  * Quick et al. (2020) Quick, L. C., Roberge, A., Mlinar, A. B., & Hedman, M. M. 2020, PASP, 132, 084402, doi: 10.1088/1538-3873/ab9504
  * Quintana et al. (2014) Quintana, E. V., Barclay, T., Raymond, S. N., et al. 2014, Science, 344, 277, doi: 10.1126/science.1249403
  * Ramirez & Kaltenegger (2014) Ramirez, R. M., & Kaltenegger, L. 2014, ApJ, 797, L25, doi: 10.1088/2041-8205/797/2/L25
  * Raymond et al. (2006) Raymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265, doi: 10.1016/j.icarus.2006.03.011
  * Raymond et al. (2007) Raymond, S. N., Scalo, J., & Meadows, V. S. 2007, ApJ, 669, 606, doi: 10.1086/521587
  * Ribas et al. (2005) Ribas, I., Guinan, E. F., Güdel, M., & Audard, M. 2005, ApJ, 622, 680, doi: 10.1086/427977
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003, doi: 10.1117/1.JATIS.1.1.014003
  * Schindler & Kasting (2000) Schindler, T. L., & Kasting, J. F. 2000, Icarus, 145, 262, doi: 10.1006/icar.2000.6340
  * Shields et al. (2016) Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Phys. Rep., 663, 1, doi: 10.1016/j.physrep.2016.10.003
  * Shkolnik & Barman (2014) Shkolnik, E. L., & Barman, T. S. 2014, AJ, 148, 64, doi: 10.1088/0004-6256/148/4/64
  * Solomon (1979) Solomon, S. C. 1979, Physics of the Earth and Planetary Interiors, 19, 168, doi: 10.1016/0031-9201(79)90081-5
  * Stelzer et al. (2022) Stelzer, B., Bogner, M., Magaudda, E., & Raetz, S. 2022, A&A, 665, A30, doi: 10.1051/0004-6361/202142088
  * Tian et al. (2014) Tian, F., France, K., Linsky, J. L., Mauas, P. J. D., & Vieytes, M. C. 2014, Earth and Planetary Science Letters, 385, 22, doi: 10.1016/j.epsl.2013.10.024
  * Tian & Ida (2015) Tian, F., & Ida, S. 2015, Nature Geoscience, 8, 177, doi: 10.1038/ngeo2372
  * Tian et al. (2008) Tian, F., Kasting, J. F., Liu, H.-L., & Roble, R. G. 2008, Journal of Geophysical Research (Planets), 113, E05008, doi: 10.1029/2007JE002946
  * Tilley et al. (2019) Tilley, M. A., Segura, A., Meadows, V., Hawley, S., & Davenport, J. 2019, Astrobiology, 19, 64, doi: 10.1089/ast.2017.1794
  * Unterborn et al. (2018) Unterborn, C. T., Desch, S. J., Hinkel, N. R., & Lorenzo, A. 2018, Nature Astronomy, 2, 297, doi: 10.1038/s41550-018-0411-6
  * Vidal-Madjar et al. (2003) Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J. M., et al. 2003, Nature, 422, 143, doi: 10.1038/nature01448
  * Villanueva et al. (2018) Villanueva, G. L., Smith, M. D., Protopapa, S., Faggi, S., & Mandell, A. M. 2018, J. Quant. Spec. Radiat. Transf., 217, 86, doi: 10.1016/j.jqsrt.2018.05.023
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Waskom (2021) Waskom, M. L. 2021, Journal of Open Source Software, 6, 3021, doi: 10.21105/joss.03021
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48, doi: 10.3847/1538-3881/aa9ff6
  * West et al. (2008) West, A. A., Hawley, S. L., Bochanski, J. J., et al. 2008, AJ, 135, 785, doi: 10.1088/0004-6256/135/3/785
  * West et al. (2015) West, A. A., Weisenburger, K. L., Irwin, J., et al. 2015, ApJ, 812, 3, doi: 10.1088/0004-637X/812/1/3
  * Wetherill (1990) Wetherill, G. W. 1990, Annual Review of Earth and Planetary Sciences, 18, 205, doi: 10.1146/annurev.ea.18.050190.001225
  * Wordsworth & Pierrehumbert (2014) Wordsworth, R., & Pierrehumbert, R. 2014, ApJ, 785, L20, doi: 10.1088/2041-8205/785/2/L20