# Water absorption in the transmission spectrum of the water-world candidate
GJ 9827 d

Pierre-Alexis Roy Department of Physics and Trottier Institute for Research on
Exoplanets, Université de Montréal, Montreal, QC, Canada  Björn Benneke
Department of Physics and Trottier Institute for Research on Exoplanets,
Université de Montréal, Montreal, QC, Canada  Caroline Piaulet Department of
Physics and Trottier Institute for Research on Exoplanets, Université de
Montréal, Montreal, QC, Canada  Michael A. Gully-Santiago Department of
Astronomy, University of Texas at Austin, Austin, TX 78712, USA  Ian J.M.
Crossfield  The University of Kansas, Department of Physics and Astronomy,
Malott Room 1082, 1251 Wescoe Hall Drive, Lawrence, KS, 66045, USA  Caroline
V. Morley Department of Astronomy, University of Texas at Austin, Austin, TX
78712, USA  Laura Kreidberg Max Planck Institute for Astronomy, Königstuhl 17,
D-69117 Heidelberg, Germany  Center for Astrophysics, Harvard & Smithsonian,
60 Garden Street, Cambridge, MA, 02138, USA  Thomas Mikal-Evans Max Planck
Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany  Jonathan
Brande Department of Physics and Astronomy, University of Kansas, 1082 Malott,
1251 Wescoe Hall Dr., Lawrence, KS 66045, USA  Simon Delisle Department of
Physics and Trottier Institute for Research on Exoplanets, Université de
Montréal, Montreal, QC, Canada  Thomas P. Greene Space Science and
Astrobiology Division, NASA Ames Research Center, MS 245-6, Moffett Field, CA
94035, USA  Kevin K. Hardegree-Ullman Steward Observatory, The University of
Arizona, Tucson, AZ 85721, USA  Travis Barman Lunar and Planetary Laboratory,
The University of Arizona, 1640 E. University Blvd, Tucson, AZ 85721, USA
Jessie L. Christiansen Caltech/IPAC, M/S 100-22, 1200 E. California Blvd.,
Pasadena, CA 91125, USA  Diana Dragomir Department of Physics and Astronomy,
University of New Mexico, 1919 Lomas Blvd. NE, Albuquerque, NM 87131, USA
Jonathan J. Fortney Department of Astronomy & Astrophysics, University of
California, Santa Cruz, CA 95064, USA  Andrew W. Howard Department of
Astronomy, California Institute of Technology, Pasadena, CA 91125, USA  Molly
R. Kosiarek Department of Astronomy and Astrophysics, University of
California, Santa Cruz, CA 95064, USA  Joshua D. Lothringer Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21210, USA

###### Abstract

Recent work on the characterization of small exoplanets has allowed us to
accumulate growing evidence that the sub-Neptunes with radii greater than
$\sim 2.5\,R_{\oplus}$ often host H2/He-dominated atmospheres both from
measurements of their low bulk densities and direct detections of their low
mean-molecular-mass atmospheres. However, the smaller sub-Neptunes in the
1.5-2.2 R⊕ size regime are much less understood, and often have bulk densities
that can be explained either by the H2/He-rich scenario, or by a volatile-
dominated composition known as the “water world” scenario. Here, we report the
detection of water vapor in the transmission spectrum of the $1.96\pm 0.08$ R⊕
sub-Neptune GJ 9827 d obtained with the Hubble Space Telescope. We observed 11
HST/WFC3 transits of GJ 9827 d and find an absorption feature at 1.4$\mu$m in
its transit spectrum, which is best explained (at 3.39$\sigma$) by the
presence of water in GJ 9827 d’s atmosphere. We further show that this feature
cannot be caused by unnoculted star spots during the transits by combining an
analysis of the K2 photometry and transit light-source effect retrievals. We
reveal that the water absorption feature can be similarly well explained by a
small amount of water vapor in a cloudy H2/He atmosphere, or by a water vapor
envelope on GJ 9827 d. Given that recent studies have inferred an important
mass-loss rate ($>0.5\,$M⊕/Gyr) for GJ 9827 d making it unlikely to retain a
H-dominated envelope, our findings highlight GJ 9827 d as a promising water
world candidate that could host a volatile-dominated atmosphere. This water
detection also makes GJ 9827 d the smallest exoplanet with an atmospheric
molecular detection to date.

Exoplanets (498); Exoplanet atmospheres (487); Planetary atmospheres (1244)

##  1 Introduction

While many questions remain regarding the nature of sub-Neptune exoplanets,
the last decade of transmission spectroscopy with the Hubble Space Telescope
(HST) has shown that the larger sub-Neptunes are often best described by
H2/He-dominated atmospheres (e.g., Benneke et al., 2019a, b, Mikal-Evans et
al., 2020, Kreidberg et al., 2022). However, this picture is much less clear
when considering sub-Neptunes that are in the smaller 1.5-2.2 R⊕ size-regime,
near the radius valley (Fulton et al., 2017, Fulton & Petigura, 2018, Van
Eylen et al., 2018, Hardegree-Ullman et al., 2020). These planets have bulk
densities than can be explained by either the H2/He-rich sub-Neptune scenario,
or by a volatile-dominated composition, where water (or another molecule of
similar mean-molecular-weight) supplants H2 and He as the most abundant
atmospheric species (Rogers & Seager, 2010, Luque & Pallé, 2022, Rogers et
al., 2023). This type of exoplanet has been long theorized and is referred to
as “water world” (Adams et al., 2008, Acuña et al., 2022).

These smaller sub-Neptunes, which are often inconsistent with extended
H-dominated atmospheres with large scale heights (Aguichine et al., 2021,
Piaulet et al., 2022), are also found in a smaller mass regime than the larger
sub-Neptunes, making these close-in planets much more exposed to mass-loss
processes (Owen, 2019), and thus more likely to have lost their H2 and He
envelope over their lifetime. A recent study found a first line of evidence
for the existence of such volatile-rich water worlds in the super-Earth
Kepler-138 d, by combining a thorough interior analysis of the planet with
mass-loss estimates, effectively showing that this super-Earth cannot be
purely rocky, but that it also cannot retain a hydrogen layer (Piaulet et al.,
2022). However, the direct spectroscopic confirmation of a volatile-rich high
mean-molecular-weight atmosphere on a water world candidate still eludes us,
and such a result would provide a new line of evidence for the water worlds.

The discovery of the transiting sub-Neptune GJ 9827 d (Niraula et al., 2017,
Rodriguez et al., 2018) represents a rich opportunity to characterize the
atmosphere of a warm sub-Neptune via transmission spectroscopy and to deepen
our understanding of this potential water world (Aguichine et al., 2021).
Rapidly orbiting (6.2 days) a low-mass K6V star with a size of $1.96\pm
0.08\,R_{\oplus}$, a mass of $3.4\pm 0.6\,M_{\oplus}$ (Kosiarek et al., 2021),
and a zero-albedo equilibrium temperature of $680\pm 25\,$K (Rodriguez et al.,
2018), GJ 9827 d allows us to obtain a high signal-to-noise ratio (S/N) in
transmission spectroscopy, and to add a precious new target to the sample of
sub-Neptunes with transit spectra. While JWST now allows to observe the
eclipses and phase curves of small exoplanets deeper in the infrared (e.g.,
Kempton et al., 2023), transit spectroscopy remains the best method to obtain
in-depth looks into the atmospheres of sub-Neptunes and potential water worlds
with HST, as they rarely are hot enough to provide high S/N in the near-
infrared (hot Neptune desert; Owen & Lai, 2018).

While the average density of GJ 9827 d has now been constrained ($>3\sigma$)
by numerous RV studies of the system (Prieto-Arranz et al., 2018, Rice et al.,
2019, Kosiarek et al., 2021), there still remains ambiguity regarding its bulk
composition, as its density could be explained by a range of compositions from
an extended H2/He layer to a water world composition with a $\sim$20 % water
mass fraction (Aguichine et al., 2021). However, the high irradiation of the
planet, the old age of the system (Kosiarek et al., 2021) and the non-
detections of HeI and H$\alpha$ absorption from the ground (Kasper et al.,
2020, Carleo et al., 2021, Krishnamurthy et al., 2023) make it unlikely that
GJ 9827 d would have retained a primordial H-dominated envelope to date.

In this work, we present the most precise look yet at GJ 9827 d via
transmission spectroscopy with the Hubble Space Telescope Wide Field Camera 3
(HST/WFC3), and reveal a water absorption feature in its transmission
spectrum. In Section 2, we present the observations obtained for this study
and we describe the data analysis in Section 3. Section 4 presents our
atmospheric analysis of the HST transit spectrum and the related results are
presented in Section 5. We end by discussing our findings and presenting our
conclusions in Section 6.

Figure 1: All 10 HST/WFC3 broadband light-curve fits of the transits of GJ
9827 d. Left: Systematics-corrected and normalized broadband light curves for
the 10 transits of GJ 9827 d (data points). Each visit is centered around the
fitted transit time for that visit. The best-fitting transit model is also
shown as the grey line. Right: Residuals of the broadband light-curve fits
shown on the left.

##  2 Observations and data reduction

GJ 9827 d was observed transiting its host star 11 times between December 2017
and December 2020 with the Wide Field Camera 3 (WFC3) onboard the Hubble Space
Telescope (HST) as part of the mini-Neptune atmosphere diversity survey (GO
15333: PI Crossfield). The G141 grism was used in order to obtain the
transmission spectrum of the sub-Neptune over the 1.1-1.7 $\mu$m range.

Figure 2: Observed mid-transit time (black points) for each of the 10 transits
compared to a fitted linear ephemeris to all transit timings presented (except
the two timings with large uncertainties; $T_{0,BJD}=2459185.987\pm 0.002$ ,
$P=6.20186\pm 0.00002$ d) . The ephemeris derived from the K2 campaign is
shown (blue; Niraula et al., 2017) along with the Spitzer transit for this
planet (red; Kosiarek et al., 2021). Our 10 HST/WFC3 transits suggest that
there are statistically significant TTVs in the orbit of GJ 9827 d of the
order of 5-10 minutes. Figure 3: Individual relative transit depths compared
to the relative final transmission spectrum. For each of the 10 transmission
spectra, we obtain the relative transit spectrum by subtracting the average
across wavelengths. We then subtract the relative final transmission spectrum
from the individual relative spectra, and display these transit depths for all
visits in each channel. The relative transit depths are consistent across the
10 visits with the final relative spectrum. Points that are significantly away
from zero (dotted lines) have larger error bars which makes them less
important in the weighted average and thus do not bias our transmission
spectrum. The uncertainty on the combined spectrum is also displayed as the
grey region. In the top left panel, the broadband transit depths are shown,
centered on the average across the 10 visits, highlighting the variability in
the observed broadband transit depths from the different visits.

Each of the 11 HST/WFC3 transit observations consisted of $\sim$3 hours of
observing time spanning three telescope orbits with $\sim$1 hour gaps between
them. Each transit observation is thus composed of one orbit before, during,
and after transit. The transit time series were obtained with the G141 grism
using the spatial scan mode. In order to optimize the duty cycle of our
observation, we used both the forward and backward detector scans. We discard
one of the transits from our analysis (November 1st 2019) since a pointing
maneuver cut orbit 1 short before the ramp was stabilized, effectively
carrying an unusually strong ramp to orbit 2 and polluting the in-transit
observations.

We reduce the observations following standard procedures for HST/WFC3
observations (details in Benneke et al., 2019a, b). In order to minimize the
background contribution, we subtract consecutive reads up-the-ramp and then
add together the background subtracted frames. We then construct flat-fielded
images from the flat-field data product provided by STScI. We use a normalized
row-added flux template in order to remove and replace outlier pixels in our
frames. We follow Benneke et al. (2019a) in order to correct for the slight
slanted shape of the trace on the detector, which is introduced by the spatial
scan mode, using a trapezoidal shape integration scheme for the wavelength
bins, which we choose to be 30 nm wide. Our flux integration does not perform
presmoothing and does account for partial pixels along the trapezoidal bin
boundaries. Finally, in order to account for the small drift of the star
across the detector during the observations, we account for a small position
shift which is measured in each frame.

##  3 Data Analysis

We perform the light-curve fitting of our 10 transits of GJ 9827 d
individually using the ExoTEP framework (Benneke et al., 2017, 2019a, 2019b).
We use ExoTEP to jointly fit the transit model with a systematics model and a
photometric noise parameter in a Markov Chain Monte Carlo scheme (Foreman-
Mackey et al., 2013). We decide to fit the transits individually since they
display variability in the transit timings (see Figures 1, 2, and Section
6.1).

Each visit in our data set consists of three HST orbits which do not cover the
full transit duration of GJ 9827 d (Figure 1). The in-transit observations
only occur during the second orbit of each visit, either observing the
ingress, the middle of transit, or the egress (Figure 1). For that reason, we
cannot obtain reliable constraints on the orbital parameters out of the
partial transit observations and decide to use a fixed orbital solution during
the fits ($b$=0.91, $a/R_{\star}=19.88$).

Since the visits do not include a burn-in orbit, we cannot follow the standard
procedure to discard the first orbit, which displays a stronger ramp in time
as the detector is still settling (e.g., Kreidberg et al., 2014). We rather
choose to keep the 2-4 last points of orbit 1 in each visit (Figure 1), as the
strong ramp has stabilized by then, and it provides essential baseline
information, especially for visits that only have mid-transit or egress data
in orbit 2 (Figure 1). For all orbits, we discard the first forward and
backward scans.

Because GJ 9827 is a close-in, near-resonant system, some of the visits in our
data set also include transits of GJ 9827 b. Given the partial coverage of our
visits, we simply remove the points where GJ 9827 b is expected to transit,
which affects visits 5 and 10. Visits 6 and 7 also include a transit of planet
b, but it happens in the first orbit which is already mostly discarded.
Finally, we remove the last 5 points of visit 3 since they are clear outliers.

###  3.1 White-light-curve fit

We fit for systematics trends in the normalized transit light curves
simultaneously with the transit model using an analytical model that allows
for a linear slope throughout the visit duration and an exponential ramp in
each orbit. Following previous work (e.g., Kreidberg et al., 2014, Benneke et
al., 2019b) we use the following parameterization to account for these
systematics:

$S_{\mathrm{model}}(t)=(c\,S(t)+v\,t_{v})\times(1-e^{(-a\,t_{\mathrm{orb}}-b-d)}).$
(1)

Here, $c$ is the normalization constant, $v$ is the linear slope throughout
the visit, $a$ and $b$ are the rate and amplitude of the exponential ramp in
each orbit, and $d$ is an offset only for first orbit reads. $S(t)$ is a
function that is equal to 1 for forward scans and to $s$ for backward scans,
allowing to correct for the offset between forward and backward scans. Finally
$t_{v}$ is the time since the start of the visit, while $t_{orb}$ is the time
since the start of the orbit.

We produce our transit light-curve models using the Batman package (Kreidberg,
2015). Since we are not trying to fit the orbital solution, the only two
astrophysical parameters that we are fitting in the transit light curve are
the transit depth, and the mid-transit time (Figure 2). For the limb
darkening, we use the Exotic-LD package (Grant & Wakeford, 2022) to compute
the coefficients using 1D stellar models (Kurucz, 1993) and the quadratic limb
darkening law. The impact parameter and semi-major axis are set to the values
in Niraula et al. (2017). We obtain the posterior distribution on our
parameters by running a MCMC analysis individually for each of the 10 visits.
We use four walkers per parameter and all priors are uniform. The 10 white-
light-curve fits are shown in Figure 1.

Instrument | Wavelength | Depth |  $\pm$1$\sigma$
---|---|---|---
| [$\mu$m] | [ppm] | [ppm]
HST/WFC3 | 1.10 – 1.13 | 941.3 | 59.7
| 1.13 – 1.16 | 1003.3 | 28.3
| 1.16 – 1.19 | 939.6 | 24.5
| 1.19 – 1.22 | 945.5 | 23.4
| 1.22 – 1.25 | 987.7 | 22.4
| 1.25 – 1.28 | 965.7 | 18.8
| 1.28 – 1.31 | 921.0 | 22.3
| 1.31 – 1.34 | 998.7 | 20.3
| 1.34 – 1.37 | 1018.9 | 22.1
| 1.37 – 1.40 | 986.8 | 21.9
| 1.40 – 1.43 | 1015.4 | 21.2
| 1.43 – 1.46 | 960.4 | 22.1
| 1.46 – 1.49 | 982.5 | 23.0
| 1.49 – 1.52 | 992.2 | 22.5
| 1.52 – 1.55 | 953.8 | 21.3
| 1.55 – 1.58 | 935.3 | 22.8
| 1.58 – 1.61 | 919.3 | 23.5
| 1.61 – 1.64 | 944.0 | 19.1
| 1.64 – 1.67 | 941.2 | 21.9
Table 1:  HST/WFC3 near-infrared combined spectrum of GJ 9827 d

###  3.2 Spectroscopic light-curve fit

We use the white-light-curve fits results in terms of the systematics model to
pre-correct the light curves in each spectral bin (divide-white method;
Stevenson et al., 2014). Thus, we divide the spectroscopic light curves by the
white-light-curve best-fitting systematics model before starting the fitting.
We produce our spectroscopic transit light-curve models similarly as in the
white-light-curve case, but we now keep the mid-transit time fixed to the
best-fitting value found by the white-light-curve fit. The limb darkening is
again modelled with Exotic-LD, and this time, our systematics model is a
three-parameter linear slope with trace position (measured during the data
reduction, see Section 2). We thus obtain 10 transmission spectra, one from
each visit, by running an MCMC analysis on each. We again use four walkers per
parameter and uniform priors on all parameters.

Figure 4: Water detection in the transmission spectrum of GJ 9827 d. Left:
Transmission spectrum of GJ 9827 d (black points) shown with our model
transmission spectra constraints from the nested sampling atmosphere retrieval
(blue) and from the photometry-informed “transit light-source effect”
retrieval (orange). The dark blue and light blue shaded regions show the
1$\sigma$ and 2$\sigma$ Bayesian credible intervals from the atmosphere
retrieval respectively. The atmospheric median transmission model is shown in
blue and the best-fitting model is shown in red. The best-fitting TLS model is
shown in orange along with the 1$\sigma$ and 2$\sigma$ Bayesian credible
intervals in light orange. Right: Joint constraints on the cloud-top pressure
versus the water mixing ratio derived from our Scarlet well-mixed retrieval.
The colored shading describes the normalized probability density as a function
of the water mixing ratio (assuming uniform vertical profiles) of the
atmosphere, and of the cloud-top pressure. The black contours highlight the 1,
2 and 3$\sigma$ Bayesian credible regions. The water abundance relative to a
solar composition atmosphere is shown on the top axis. The posterior
probability distribution allows for multiple atmospheric scenarios ranging
from H2/He envelopes with small amounts of water to water-dominated envelopes.
The blue points identify two representative samples of these two scenarios
which are displayed in Figure 6.

###  3.3 Combining all the visits together

We compute a weighted average of our 10 individual transmission spectra to
obtain our final transmission spectrum of GJ 9827 d. In order to verify the
robustness of our spectrum, and to ensure that it is not affected by that
variability in the observations, we compare the relative transit depth in each
channel, for each visit (Figure 3). From each spectrum, we subtract the
weighted average (across wavelengths) of said spectrum, essentially making it
a relative transit spectrum centered around zero. We then subtract the
relative combined spectrum of all visits (also centered around zero) from each
individual spectrum to effectively center each spectroscopic transit depth
around zero. We then inspect this relative transit depth for each
spectroscopic bin and for each visit, in order to ensure that the points in
our final transmission spectrum are not affected by outliers (Figure 3). We
find that for each spectroscopic channel, all visits mostly agree with the
weighted average within error bars, and points that are inconsistent have
larger error bars, which makes them much less important in the weighted
average, since the weight of each points is inversely proportional to the
uncertainty squared (Figure 3). The final average transmission spectrum is
presented in Table 1 and in Figure 4. We decide to discard the last
spectroscopic channel (1.67-1.70 $\mu$m) since it is systematically lower than
the rest of the spectrum, and is near the edge of the trace on the detector
where the data is less reliable.

##  4 Atmospheric Modeling

We perform atmosphere retrievals on our GJ 9827 d transmission spectrum using
the SCARLET framework (Benneke & Seager, 2012, 2013, Knutson et al., 2014,
Kreidberg et al., 2014, Benneke, 2015, Benneke et al., 2019a, b, Pelletier et
al., 2021, Roy et al., 2022). SCARLET parameterizes the molecular abundances,
the cloud deck pressure and the temperature to fit our spectrum. SCARLET uses
a Bayesian nested sampling analysis (with single ellipsoïd sampling; Skilling,
2004, 2006) to obtain the posterior probability distribution of our parameter
space and the Bayesian evidence of our models.

Figure 5: Posterior probability distributions of the free parameters used in
the SCARLET free chemistry retrieval. The diagonal panels show the
marginalized probability distributions of all individual parameters, whereas
the off-diagonal panels show the marginalized probability distributions for
each pair of parameters as colored shading. The 1, 2, and 3$\sigma$ contours
are shown in the 2D posteriors. Water is the only molecule detected in our
retrieval analysis.

For each set of parameters, SCARLET produces a forward atmosphere model in
hydrostatic equilibrium (Benneke, 2015), computes the opacity associated to
each molecule throughout the 40 vertical (pressure) layers of the model,
computes the transmission spectrum for that model and finally performs the
likelihood evaluation. The model transmission spectra produced at each step
have a resolution of 16000 and are then convolved to the wavelength bins of
the observed spectrum. The molecules considered in our retrieval are H2O, CH4,
CO2, CO and N2, as well as H2 and He which are not parameterized and rather
fill up the atmosphere (Benneke & Seager, 2013).

Retrieval model | Evidence | Bayes Factor | Nσ
---|---|---|---
| ln(Zi) | B = Zref/Zi |
All molecules | -90.68 | ref. | ref.
\+ clouds |  |  |
H2O removed | -94.96 | 72.52 | 3.39
CH4 removed | -90.24 | 0.65 | 0.90
CO2 removed | -90.45 | 0.79 | 0.90
CO removed | -90.60 | 0.92 | 0.90
N2 removed | -90.73 | 1.06 | 1.14
Clouds removed | -90.78 | 1.10 | 1.23
H2O, CH4 removed | -94.88 | 66.92 | 3.36
Flat spectrum | -97.11 | 620.76 | 4.01
Table 2:  Bayesian model comparison results from our SCARLET atmosphere
retrievals in the free chemistry setting

We assume well-mixed vertical chemical profiles, where the abundances of
molecules do not vary throughout the atmosphere. We choose a log-uniform prior
on the abundance of each molecule ranging from 10-10 to 1 in volume mixing
ratio. We assume a constant temperature structure throughout the atmosphere
and use a Gaussian prior centered on the planet’s equilibrium temperature (680
$\pm$ 100 K; Rodriguez et al., 2018) on that parameter. The parameterization
also includes a cloud deck top pressure, which blocks all light rays going
through that pressure level. We again use a log-uniform prior on that
parameter from 10-4 mbar to 105 mbar. Thus, our atmosphere retrieval includes
seven free parameters in total (or less when molecules are removed, see Table
2).

Figure 6: HST/WFC3 transmission spectrum of GJ 9827 d (data points) along with
SCARLET forward atmosphere models (colored lines). Top: Two samples from our
well-mixed retrieval analysis (Figure 4) are shown, representing the mini-
Neptune scenario with a cloudy H2/He atmosphere composed of $\sim$1% water
(pale blue) and a water world scenario with a water-rich atmosphere (dark
blue). Middle: A secondary atmosphere model for a pure methane envelope is
also shown (red) in order to highlight the methane absorption features.
Chemically consistent models (still assuming a uniform temperature profile)
are shown for a cloud-free, solar composition case (solar metallicity, solar
C/O; orange) and for a cloudy case with C/O=0.1 and a 100 $\times$ solar
metallicity (green). The observed spectrum is inconsistent with cloud-free
low-metallicity scenarios and prefers water absorption features to methane
absorption features, mainly around 1.2 and 1.65 $\mu$m. The strength of the
features in the spectrum is also displayed in units of H/He scale heights
(right axis). Bottom: The best-fit model from the retrieval analysis is shown
(pale blue), along with the transmission spectrum of the same model once the
water opacity is turned off (grey). The contribution of water opacity to the
spectral signatures is highlighted in blue. We also present a binned version
of the transmission spectrum where points are binned together by pair with the
exception of the blue-most channel.

##  5 Results

The observed transit spectrum of GJ 9827 d displays a water absorption feature
at 1.4 $\mu$m (Figure 4). Qualitatively, the transit depths in the spectrum
are deeper in the 1.4 $\mu$m water band followed by a downward slope that
follows the wing of the water absorption feature. Quantitatively, a Bayesian
model comparison analysis of our well-mixed retrievals (Benneke & Seager,
2013) prefers models that include the presence and absorption of water with a
significance of 3.39$\sigma$ (Bayes Factor = 72.52; Table 2) compared to
models that do not include water.

###  5.1 Metallicity-clouds degeneracy

Our free chemistry retrievals show that the data can both be explained by a
water-rich atmosphere, where water is the most abundant molecule, as well as
with a H2/He-dominated atmosphere that still contains a small amount of water
(Figure 4). At 1$\sigma$, models with a water mixing ratio between 0.02% and
80% are preferred by the spectrum. When compared to the amount of water in a
solar metallicity envelope, we see that this interval in abundance ranges from
1 $\times$ solar metallicity models, which are dominated by H2 and He gas, to
1000 $\times$ solar metallicity models, where water is now the principal
species (Figure 4).

In the cases where water is present in small amounts, a cloud deck is needed
to explain the observed transit spectrum (Figures 4, 5). This is due to the
fact that the observed spectrum does not display the large amplitude expected
from cloud-free models (Figure 6). This need for clouds in low mean-molecular-
weight atmospheres is also seen in the marginalized probability distributions
of the other molecules (Figure 5), since their spectral features are
inconsistent with the observed spectrum, and thus must be muted by clouds in
order to yield high-likelihood models. In cases where the water is more
abundant and becomes the principal molecule, the spectra naturally display
muted features because of the high density of the atmospheres and lower
atmospheric scale height, thus removing the need for high clouds. In this
water-rich scenario, the constraint on the cloud-top pressure disappears,
explaining the observed posterior distribution (all deep cloud decks become
equally consistent; Figure 4).

###  5.2 Upper limits on other molecules

While methane is known to display a similar 1.4 $\mu$m absorption feature as
water (Bézard et al., 2022), it remains disfavored in our free chemistry
retrieval analysis. The spectral signatures of methane not only include a
feature around 1.4 $\mu$m, but also at 1.2 $\mu$m and at 1.7 $\mu$m, as shown
by SCARLET forward atmosphere models for a pure methane envelope and for a
solar composition atmosphere (Figure 6). However, the observed transmission
spectrum of GJ 9827 d does not display these absorption features at 1.2 and
1.7 $\mu$m (Figure 6) and is in better agreement with the water models that
display a smaller feature at 1.2 $\mu$m, no absorption at 1.7 $\mu$m, and a
broader feature at 1.4 $\mu$m. In order to obtain chemically consistent models
that agree with the observed transit spectrum, the carbon-to-oxygen (C/O)
ratio must be decreased in order to favor O-bearing molecules (here, water)
and a cloud deck must be included to mute the spectral amplitude of the water
absorption features. Such models (e.g. C/O = 0.1, Metallicity = 100 $\times$
solar and pCloud = 1 mbar) yield qualitatively and quantitatively similar
transmission spectra to those favored by our free chemistry retrieval (Figures
4, 6).

No other molecule besides water is statistically detected by our retrievals
(Table 2, Figure 5). However, we can derive upper limits in their abundances
based on our results from the free chemistry retrievals, either from the non-
detection of specific spectral features (e.g., CH4) or because too much of
anyone species increases the mean-molecular-weight of the atmosphere,
eventually yielding a flat spectrum (e.g., N2). Thus, our retrievals allow us
to constrain the upper limits on the mixing ratios of CH4, CO2, CO and N2 to
3.04, 19.4, 52.5 and 60.4 % at 3$\sigma$ significance.

###  5.3 Significance of a featureless spectrum

In order to evaluate how our spectrum deviates from a featureless spectrum, we
compute the deviation from the best-fitting straight line using $\chi^{2}$
statistics. Using the binned version (bottom of Figure 6), we obtain that the
transmission spectrum of GJ 9827 d deviates from a straight line at
3.24$\sigma$. In order to assess how our water detection compares to a
featureless flat spectrum within the Bayesian paradigm, we compute the
Bayesian evidence of a one-parameter flat line model
$\mathcal{Z}_{\textrm{flat}}$. Given the simplicity of this one-parameter
model, we do not need to use SCARLET nested sampling to obtain that value, and
rather numerically estimate it via the following analytical solution:

$\displaystyle\mathcal{Z}_{\textrm{flat}}$
$\displaystyle=\left(\frac{1}{\sqrt{2\pi}}\right)^{N}\times\left(\prod_{i=1}^{N}\frac{1}{\sigma_{i}}\right)\times\left(\frac{1}{\theta_{2}-\theta_{1}}\right)$
(2)
$\displaystyle\qquad\qquad\times\int_{\theta_{1}}^{\theta_{2}}\exp\left(-\sum_{i=1}^{N}\frac{(D_{i}-\theta)^{2}}{2\sigma_{i}^{2}}\right)d\theta,$

where $N$ is the number of points in the spectrum, $\sigma_{i}$ is the
uncertainty of the $i^{\mathrm{th}}$ point of the spectrum denoted $D_{i}$,
and where $\theta_{1}$ and $\theta_{2}$ are the limits of our uniform prior on
the transit depth parameter $\theta$. This allows us to show that our
atmosphere model is preferred to the flat spectrum model at 4.01$\sigma$
(Bayes Factor = 620.76; Table 2).

###  5.4 Ruling out stellar contamination

The Transit Light Source Effect (TLS, Rackham et al., 2018) can mimic water
features at the 20 ppm level or more for modestly spotted stars under certain
observational configurations, and could thus create the feature observed in
our transmission spectrum. The best constraint available for the starspot
coverage and contrast for GJ 9827 comes from the K2 Campaign 12 (C12)
lightcurve, with an SFF-derived (Vanderburg & Johnson, 2014) peak-to-valley
amplitude of 0.45%, slightly lower than the typical 0.7% for K6 spectral types
(Rackham et al., 2019). Coarse scaling laws can relate the observed K2
amplitude to the surface coverage of starspots, under assumptions of size,
number, and location of spots on a stellar surface (Rackham et al., 2018, Guo
et al., 2018). For a K6 spectral type, 0.45% amplitude variations relate
(conservatively) to spot-covering fractions of 1-4% (Rackham et al., 2019).
Thus, we adopt a spot contrast typical of a K6 star and a filling factor of
1-4% (Rackham et al., 2019). Under these assumptions, a planet with a 1%
transit depth could expect an H2O contamination of $<$15 ppm from unocculted
spots. GJ 9827 d’s much smaller transit depth of $<0.1\%$ would therefore
yield a negligible TLS water contamination of $<1.5$ ppm.

We perform a retrieval on our combined transit spectrum in which we fit for
TLS parameters rather than for atmospheric parameters. We use the 1-4% spots
filling factor as a prior in this TLS retrieval, and the other parameters are
the photospheric temperature, the spots’ temperature difference and a scaling
factor. Our TLS models simulate unocculted star spots by averaging the stellar
spectrum of the star’s photoshere with the spectrum of the cooler spots
(weighted by the spots coverage) based on Phoenix stellar models (Husser et
al., 2013). It then simulates the transit of an airless planet to obtain the
stellar-contaminated transit spectrum. When using the photometry-derived
prior, we find that the TLS fitting is restricted to flat models which
indicate that stellar contamination is not expected for that system (Figure
4).

Repeating the same TLS retrieval with nonrestrictive uniform priors on all
parameters further demonstrates that the signal cannot be caused by the star.
We run the same TLS retrieval as described above, but without using the 1-4%
spots filling factor prior, to see under what stellar conditions the signal
can be explained by the star. We find that, in order for the TLS to reproduce
the signal in our spectrum, not only does the spots parameters need to be
unrealistically large (73% spots coverage and $<-794\,$K spot temperature
difference), but the model needs to adopt a strong positive ramp towards short
wavelengths (as in Moran et al., 2023), which becomes inconsistent with the K2
transit depth measurement. We thus conclude that stellar contamination cannot
explain the feature in the transmission spectrum of GJ 9827 d, and that the
water detection comes from the planet’s atmosphere.

##  6 Discussion and Conclusions

The HST/WFC3 transmission spectrum of GJ 9827 d presented in this work
provides a precious target in the population of sub-Neptune exoplanets for
which we have precise transmission spectra; and highlights GJ 9827 d as the
smallest exoplanet with an unambiguous atmospheric molecular detection to
date. Compared to the other characterized sub-Neptunes, our detection of a
$\sim$1 H/He scale height water feature (Figure 6) makes it stronger than for
similarly hot sub-Neptunes, although it remains broadly consistent with the
previously observed trend where hotter sub-Neptunes display stronger H2O
amplitudes (Crossfield & Kreidberg, 2017). Our analysis of the 10 observations
of GJ 9827 d’s transits also revealed some variability in this planet’s orbit,
which is to be considered for further monitoring of the system with state-of-
the-art facilities such as JWST. Finally, our detection of a water feature in
GJ 9827 d’s transit spectrum provides the first detection of water in the
atmosphere of a potential water world, which, when combined with GJ 9827 d’s
large mass-loss rate, provides a first line of evidence for this sub-Neptune
hosting a water-steam dominated atmosphere.

###  6.1 Variability in the transits of GJ 9827 d

The analysis of the 10 transits of GJ 9827 d with HST revealed a significant
variability in the transit timings observed from one visit to the other
(Figure 2). While this variation is not surprising for a near-resonant system
and did not impact the features observed in the transit spectrum (as shown by
our analysis of the relative spectra; Figure 3), it still is in contrast with
the previously observed TTVs for this system which were of the order of
$\sim$3 minutes (Niraula et al., 2017). However, the 5-10 minutes TTVs
observed for planet d in this work are consistent with an independent study of
the TTVs of the GJ 9827 system combining all photometric and radial velocity
data (Livingston et al., in prep.).

The limited number of in-transit data points in the time series in this
program could also explain the range of transit depths observed in our
results. As described earlier, each HST orbit displays an exponential ramp in
time that is fitted by our systematics model. This ramp has a much stronger
effect in the first few integrations than in the last few integrations of each
orbit, when it has settled. Thus, HST/WFC3 observations inherently provide
better quality observations towards the end of each orbit. When considering
the individual visits in our program, it seems that egress visits yield deeper
transit depths than mid-transit or ingress visits (Figures 1, 3). This could
then be explained by the fact that the different visits in our data set have
varying in-transit data quality depending on whether the late-orbit
integrations are in the transit (ingress and mid-transit visits) or are in the
baseline (egress visits). For instance, visit 6, which is an egress
observation, displays a deeper transit depth than in the other visits.
However, the relative shape of the transmission spectrum, which is the
quantity of interest for this work, is consistent with the other visits
(Figure 3).

Another potential source of the TTV and transit depth variability observed in
our program is star-spot crossing. GJ 9827 has been shown to display quasi-
periodic flux ($\sim$0.45%) variations with a period of $\sim$30 days
(Rodriguez et al., 2018, Teske et al., 2018, Prieto-Arranz et al., 2018, Rice
et al., 2019, Kosiarek et al., 2021). If these stellar variations were caused
by stellar spots, then using a fixed orbital solution and a transit model that
does not include the effect of spots in our light-curve fits could lead to
biases in our retrieved parameters.

Because of the variability discussed above, we decided to fix the orbital
solution and limb darkening coefficients for the light-curve fits in our
program (Section 3). In order to ensure that the limb darkening coefficients
and stellar parameters chosen do not affect our atmospheric inference, we
repeat our light-curve fits for multiple assumptions on the limb darkening.
Using quadratic limb darkening coefficients, we reproduce the same fitting but
using a 3D stellar model (Magic et al., 2015) when computing the coefficients.
We further try the light-curve fits by varying the effective stellar
temperature to the +1$\sigma$ and $-1\sigma$ values of that parameter
(Kosiarek et al., 2021). In all cases, we find that the limb darkening and
choice of stellar parameters only affect the retrieved spectrum with a
constant offset throughout the wavelength range, and that the relative spectra
all show the water absorption feature and are all consistent within 1$\sigma$.
This thus shows that our choice for the stellar and limb darkening parameters
does not affect the shape of the transmission spectrum, and subsequently, our
atmospheric analysis.

Similarly, given the difficult observational setting, we test the robustness
of the spectrum to the systematics models used by trying two alternative
light-curve fitting methods. First, we start from the divide-white corrected
spectroscopic light curves and jointly fit a relative transmission spectrum.
To do so, we jointly fit (across visits) the relative transit depth in each
channel, where the broadband average of the individually fitted spectra is
subtracted for each visit (since there sometimes are discrepancies between the
white-light-curve transit depths, and average spectral depths in HST/WFC3
data). Secondly, we repeat the method presented in Section 3, but using the
RECTE systematics model and orbit 1 (and not using the divide-white method;
Zhou et al., 2017) for the seven visits which are not affected by transits of
planet b in orbits 1 or 2. This gives us 7 transmission spectra that we
combine with a weighted average. Both of these methods produce relative
transmission spectra that are consistent within uncertainties with the one
presented in Table 1. The spectrum obtained from the RECTE models has
increased uncertainties, both from the smaller number of visits and from
increased fitted scatter in some visits, but is still in agreement with the
other two spectra. The spectrum presented in this work is thus robust to the
choice of systematics model.

###  6.2 Water in the envelope of a potential water world

The water detection in the transit spectrum of GJ 9827 d makes it the first
water world candidate with an atmospheric water detection consistent with a
water-rich envelope. It thus positions itself in the sample of potential water
worlds, with other small sub-Neptunes and super-Earths such as Kepler-138 d
(Piaulet et al., 2022), L 98-59 d (Kostov et al., 2019), TOI-1685 b (Bluhm et
al., 2021), GJ 3090 d (Almenara et al., 2022), TOI-270 d (Günther et al.,
2019, Mikal-Evans et al., 2023). In contrast, a water feature was also
detected in the transit spectrum of TOI-270 d (Mikal-Evans et al., 2023), but
the analysis revealed that the H-rich atmosphere scenario was favored for this
sub-Neptune, showing that the line can be fine between a mini-Neptune and a
water world.

With its small mass of 3.42 M⊕ (Kosiarek et al., 2021) and its proximity to
its host star (6.2 d orbit), the estimated mass-loss rate of GJ 9827 d is
$>$0.5 M⊕/Gyr (Krishnamurthy et al., 2023). With an estimated age around 6 Gyr
(Kosiarek et al., 2021), GJ 9827 d is unlikely to retain an extended H2/He
envelope today. Furthermore, monitoring of GJ 9827 d’s spectrum in the search
of H$\alpha$ and HeI signatures with Keck/NIRSPEC (Kasper et al., 2020)
CARMENES (Carleo et al., 2021) and IRD (Krishnamurthy et al., 2023) resulted
in no evidence of an extended H2/He atmosphere around the planet. Hence, the
H-rich scenario with a smaller water abundance would provide a somewhat
contradictory statement to the previous studies that observed GJ 9827 d from
the ground. However, the water-rich scenario can both explain the observed HST
transit spectrum, as well as the non-detection of H$\alpha$/HeI lines from
ground-based studies. The water-dominated envelope is thus the compositional
scenario that explains all of the data at hand on this system in the most
natural way.

In this water-rich scenario, GJ 9827 d would thus represent a larger, hotter,
close-in version of the icy moons of the giant planets in the solar system.
Indeed, water is believed to be the dominant volatile of the icy moons of the
solar system (Schubert et al., 2004). GJ 9827 d could then have formed outside
of the water ice line, where water ice is available in large amounts as a
planetary building block (Mousis et al., 2019, Venturini et al., 2020). It
could then have migrated towards its current stable near-resonant orbit,
during which the increasingly important stellar irradiation would have driven
an important H2/He loss, and it would be observed today with a high mean-
molecular-weight water vapor atmosphere due to its warm temperature (Adams et
al., 2008, Pierrehumbert, 2023) and its H2/He depletion.

While our transmission spectrum cannot unambiguously distinguish between the
H-rich and H-depleted scenarios, we have provided the first water detection in
the envelope of a water world candidate, making it a key target for further
monitoring with JWST. Transmission spectroscopy of GJ 9827 d with NIRISS/SOSS
and NIRSpec/G395H would provide the high-precision continuous viewing of the
full transit of the planet that is needed to explain the variability observed
with HST, as well as provide a more precise transmission spectrum that could
not only probe the water absorption bands, but also probe for the presence of
carbon bearing species like CO and CO2 above 4 $\mu$m. A JWST transmission
spectrum of GJ 9827 d would thus lift the degeneracy observed in our study
(Figure 4) and potentially confirm the water world nature of this sub-Neptune,
simultaneously yielding the first direct detection of a water vapor dominated
envelope.

All of the data presented in this paper were obtained from the Mikulski
Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.
The specific observations analyzed can be accessed via
https://doi.org/10.17909/dvqh-2r48 (catalog 10.17909/dvqh-2r48). We wish to
thank the reviewer for the insightful comments which enhanced the quality of
our manuscript. P.-A. R. further thanks L. Bazinet, L.-P. Coulombe and S.
Pelletier for their help, comments and ideas throughout the multiple
iterations of this work. This work is based on observations with the NASA/ESA
HST, obtained at the Space Telescope Science Institute (STScI) operated by
AURA, Inc. P.-A.R. and B.B. acknowledge financial support from the Natural
Sciences and Engineering Research Council (NSERC) of Canada. P.-A.R. further
acknowledges support from the University of Montreal, and from the Trottier
institute for exoplanets (iREx). B.B. also acknowledges financial support from
the Fond de Recherche Québécois-Nature et Technologie (FRQNT; Québec).

## References

  * Acuña et al. (2022) Acuña, L., Lopez, T. A., Morel, T., et al. 2022, Astronomy & Astrophysics, 660, A102, doi: 10.1051/0004-6361/202142374
  * Adams et al. (2008) Adams, E. R., Seager, S., & Elkins‐Tanton, L. 2008, The Astrophysical Journal, 673, 1160, doi: 10.1086/524925
  * Aguichine et al. (2021) Aguichine, A., Mousis, O., Deleuil, M., & Marcq, E. 2021, The Astrophysical Journal, 914, 84, doi: 10.3847/1538-4357/abfa99
  * Almenara et al. (2022) Almenara, J. M., Bonfils, X., Otegi, J. F., et al. 2022, Astronomy & Astrophysics, 665, A91, doi: 10.1051/0004-6361/202243975
  * Benneke (2015) Benneke, B. 2015, arXiv:1504.07655 [astro-ph].  http://arxiv.org/abs/1504.07655
  * Benneke & Seager (2012) Benneke, B., & Seager, S. 2012, The Astrophysical Journal, 753, 100, doi: 10.1088/0004-637X/753/2/100
  * Benneke & Seager (2013) —. 2013, The Astrophysical Journal, 778, 153, doi: 10.1088/0004-637X/778/2/153
  * Benneke et al. (2017) Benneke, B., Werner, M., Petigura, E., et al. 2017, The Astrophysical Journal, 834, 187, doi: 10.3847/1538-4357/834/2/187
  * Benneke et al. (2019a) Benneke, B., Knutson, H. A., Lothringer, J., et al. 2019a, Nature Astronomy, 3, 813, doi: 10.1038/s41550-019-0800-5
  * Benneke et al. (2019b) Benneke, B., Wong, I., Piaulet, C., et al. 2019b, The Astrophysical Journal Letters, 887, L14, doi: 10.3847/2041-8213/ab59dc
  * Bluhm et al. (2021) Bluhm, P., Pallé, E., Molaverdikhani, K., et al. 2021, Astronomy & Astrophysics, 650, A78, doi: 10.1051/0004-6361/202140688
  * Bézard et al. (2022) Bézard, B., Charnay, B., & Blain, D. 2022, Nature Astronomy, 6, 537, doi: 10.1038/s41550-022-01678-z
  * Carleo et al. (2021) Carleo, I., Youngblood, A., Redfield, S., et al. 2021, The Astronomical Journal, 161, 136, doi: 10.3847/1538-3881/abdb2f
  * Crossfield & Kreidberg (2017) Crossfield, I. J. M., & Kreidberg, L. 2017, The Astronomical Journal, 154, 261, doi: 10.3847/1538-3881/aa9279
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306, doi: 10.1086/670067
  * Fulton & Petigura (2018) Fulton, B. J., & Petigura, E. A. 2018, The Astronomical Journal, 156, 264, doi: 10.3847/1538-3881/aae828
  * Fulton et al. (2017) Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, The Astronomical Journal, 154, 109, doi: 10.3847/1538-3881/aa80eb
  * Grant & Wakeford (2022) Grant, D., & Wakeford, H. R. 2022, Exo-TiC/ExoTiC-LD: ExoTiC-LD v3.0.0, v3.0.0, Zenodo, doi: 10.5281/zenodo.7437681
  * Guo et al. (2018) Guo, Z., Gully-Santiago, M., & Herczeg, G. J. 2018, The Astrophysical Journal, 868, 143, doi: 10.3847/1538-4357/aaeb9b
  * Günther et al. (2019) Günther, M. N., Pozuelos, F. J., Dittmann, J. A., et al. 2019, Nature Astronomy, 3, 1099, doi: 10.1038/s41550-019-0845-5
  * Hardegree-Ullman et al. (2020) Hardegree-Ullman, K. K., Zink, J. K., Christiansen, J. L., et al. 2020, The Astrophysical Journal Supplement Series, 247, 28, doi: 10.3847/1538-4365/ab7230
  * Husser et al. (2013) Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, Astronomy & Astrophysics, 553, A6, doi: 10.1051/0004-6361/201219058
  * Kasper et al. (2020) Kasper, D., Bean, J. L., Oklopčić, A., et al. 2020, The Astronomical Journal, 160, 258, doi: 10.3847/1538-3881/abbee6
  * Kempton et al. (2023) Kempton, E. M.-R., Zhang, M., Bean, J. L., et al. 2023, Nature, 1, doi: 10.1038/s41586-023-06159-5
  * Knutson et al. (2014) Knutson, H. A., Benneke, B., Deming, D., & Homeier, D. 2014, Nature, 505, 66, doi: 10.1038/nature12887
  * Kosiarek et al. (2021) Kosiarek, M. R., Berardo, D. A., Crossfield, I. J. M., et al. 2021, The Astronomical Journal, 161, 47, doi: 10.3847/1538-3881/abca39
  * Kostov et al. (2019) Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, The Astronomical Journal, 158, 32, doi: 10.3847/1538-3881/ab2459
  * Kreidberg (2015) Kreidberg, L. 2015, Publications of the Astronomical Society of the Pacific, 127, 1161, doi: 10.1086/683602
  * Kreidberg et al. (2014) Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Nature, 505, 69, doi: 10.1038/nature12888
  * Kreidberg et al. (2022) Kreidberg, L., Mollière, P., Crossfield, I. J. M., et al. 2022, The Astronomical Journal, 164, 124, doi: 10.3847/1538-3881/ac85be
  * Krishnamurthy et al. (2023) Krishnamurthy, V., Hirano, T., Gaidos, E., et al. 2023, Monthly Notices of the Royal Astronomical Society, 521, 1210, doi: 10.1093/mnras/stad404
  * Kurucz (1993) Kurucz, R.-L. 1993, Kurucz CD-Rom, 13
  * Luque & Pallé (2022) Luque, R., & Pallé, E. 2022, Science, 377, 1211, doi: 10.1126/science.abl7164
  * Magic et al. (2015) Magic, Z., Chiavassa, A., Collet, R., & Asplund, M. 2015, Astronomy & Astrophysics, 573, A90
  * Mikal-Evans et al. (2020) Mikal-Evans, T., Crossfield, I. J. M., Benneke, B., et al. 2020, The Astronomical Journal, 161, 18, doi: 10.3847/1538-3881/abc874
  * Mikal-Evans et al. (2023) Mikal-Evans, T., Madhusudhan, N., Dittmann, J., et al. 2023, The Astronomical Journal, 165, 84, doi: 10.3847/1538-3881/aca90b
  * Moran et al. (2023) Moran, S. E., Stevenson, K. B., Sing, D. K., et al. 2023, The Astrophysical Journal Letters, 948, L11, doi: 10.3847/2041-8213/accb9c
  * Mousis et al. (2019) Mousis, O., Ronnet, T., & Lunine, J. I. 2019, The Astrophysical Journal, 875, 9, doi: 10.3847/1538-4357/ab0a72
  * Niraula et al. (2017) Niraula, P., Redfield, S., Dai, F., et al. 2017, The Astronomical Journal, 154, 266, doi: 10.3847/1538-3881/aa957c
  * Owen (2019) Owen, J. E. 2019, Annual Review of Earth and Planetary Sciences, 47, 67, doi: 10.1146/annurev-earth-053018-060246
  * Owen & Lai (2018) Owen, J. E., & Lai, D. 2018, Monthly Notices of the Royal Astronomical Society, 479, 5012, doi: 10.1093/mnras/sty1760
  * Pelletier et al. (2021) Pelletier, S., Benneke, B., Darveau-Bernier, A., et al. 2021, The Astronomical Journal, 162, 73, doi: 10.3847/1538-3881/ac0428
  * Piaulet et al. (2022) Piaulet, C., Benneke, B., Almenara, J. M., et al. 2022, Nature Astronomy, 1, doi: 10.1038/s41550-022-01835-4
  * Pierrehumbert (2023) Pierrehumbert, R. T. 2023, The Astrophysical Journal, 944, 20, doi: 10.3847/1538-4357/acafdf
  * Prieto-Arranz et al. (2018) Prieto-Arranz, J., Palle, E., Gandolfi, D., et al. 2018, Astronomy & Astrophysics, 618, A116, doi: 10.1051/0004-6361/201832872
  * Rackham et al. (2018) Rackham, B. V., Apai, D., & Giampapa, M. S. 2018, The Astrophysical Journal, 853, 122, doi: 10.3847/1538-4357/aaa08c
  * Rackham et al. (2019) —. 2019, The Astronomical Journal, 157, 96, doi: 10.3847/1538-3881/aaf892
  * Rice et al. (2019) Rice, K., Malavolta, L., Mayo, A., et al. 2019, Monthly Notices of the Royal Astronomical Society, 484, 3731, doi: 10.1093/mnras/stz130
  * Rodriguez et al. (2018) Rodriguez, J. E., Vanderburg, A., Eastman, J. D., et al. 2018, The Astronomical Journal, 155, 72, doi: 10.3847/1538-3881/aaa292
  * Rogers et al. (2023) Rogers, J. G., Schlichting, H. E., & Owen, J. E. 2023, The Astrophysical Journal Letters, 947, L19, doi: 10.3847/2041-8213/acc86f
  * Rogers & Seager (2010) Rogers, L. A., & Seager, S. 2010, The Astrophysical Journal, 716, 1208, doi: 10.1088/0004-637X/716/2/1208
  * Roy et al. (2022) Roy, P.-A., Benneke, B., Piaulet, C., et al. 2022, The Astrophysical Journal, 941, 89, doi: 10.3847/1538-4357/ac9f18
  * Schubert et al. (2004) Schubert, G., Spohn, T., & McKinnon, W. B. 2004
  * Skilling (2004) Skilling, J. 2004, in AIP Conference Proceedings, Vol. 735 (Garching (Germany): AIP), 395–405, doi: 10.1063/1.1835238
  * Skilling (2006) Skilling, J. 2006, Bayesian Analysis, 1, 833, doi: 10.1214/06-BA127
  * Stevenson et al. (2014) Stevenson, K. B., Bean, J. L., Seifahrt, A., et al. 2014, The Astronomical Journal, 147, 161, doi: 10.1088/0004-6256/147/6/161
  * Teske et al. (2018) Teske, J. K., Wang, S., Wolfgang, A., et al. 2018, The Astronomical Journal, 155, 148, doi: 10.3847/1538-3881/aaab56
  * Vanderburg & Johnson (2014) Vanderburg, A., & Johnson, J. A. 2014, Publications of the Astronomical Society of the Pacific, 126, 948, doi: 10.1086/678764
  * Van Eylen et al. (2018) Van Eylen, V., Agentoft, C., Lundkvist, M. S., et al. 2018, Monthly Notices of the Royal Astronomical Society, 479, 4786, doi: 10.1093/mnras/sty1783
  * Venturini et al. (2020) Venturini, J., Guilera, O. M., Haldemann, J., Ronco, M. P., & Mordasini, C. 2020, Astronomy & Astrophysics, 643, L1, doi: 10.1051/0004-6361/202039141
  * Zhou et al. (2017) Zhou, Y., Apai, D., Lew, B. W. P., & Schneider, G. 2017, The Astronomical Journal, 153, 243, doi: 10.3847/1538-3881/aa6481