# JWST/NIRCam Transmission Spectroscopy of the Nearby Sub-Earth GJ 341b

James Kirk Department of Physics, Imperial College London, Prince Consort
Road, London, SW7 2AZ, UK James Kirk<EMAIL_ADDRESS>Kevin B.
Stevenson Johns Hopkins APL, Laurel, MD, 20723, USA  Guangwei Fu Department of
Physics & Astronomy, Johns Hopkins University, Baltimore, MD, USA  Jacob
Lustig-Yaeger Johns Hopkins APL, Laurel, MD, 20723, USA  Sarah E. Moran
Department of Planetary Sciences and Lunar and Planetary Laboratory,
University of Arizona, Tuscon, AZ, USA  Sarah Peacock University of Maryland,
Baltimore County, MD 21250, USA  NASA Goddard Space Flight Center, Greenbelt,
MD 20771, USA  Munazza K. Alam Carnegie Earth & Planets Laboratory,
Washington, DC, 20015, USA  Space Telescope Science Institute, Baltimore, MD
21218, USA  Natasha E. Batalha NASA Ames Research Center, Moffett Field, CA,
USA  Katherine A. Bennett Department of Earth & Planetary Sciences, Johns
Hopkins University, Baltimore, MD, USA  Junellie Gonzalez-Quiles Department of
Earth & Planetary Sciences, Johns Hopkins University, Baltimore, MD, USA
Mercedes López-Morales Center for Astrophysics ${\rm\mid}$ Harvard &
Smithsonian, 60 Garden St, Cambridge, MA 02138, USA  Joshua D. Lothringer
Department of Physics, Utah Valley University, Orem, UT, 84058 USA  Ryan J.
MacDonald Department of Astronomy, University of Michigan, 1085 S. University
Ave., Ann Arbor, MI 48109, USA  NHFP Sagan Fellow  E. M. May Johns Hopkins
APL, Laurel, MD, 20723, USA  L. C. Mayorga Johns Hopkins APL, Laurel, MD,
20723, USA  Zafar Rustamkulov Department of Earth & Planetary Sciences, Johns
Hopkins University, Baltimore, MD, USA  David K. Sing Department of Earth &
Planetary Sciences, Johns Hopkins University, Baltimore, MD, USA  Department
of Physics & Astronomy, Johns Hopkins University, Baltimore, MD, USA  Kristin
S. Sotzen Johns Hopkins APL, Laurel, MD, 20723, USA  Jeff A. Valenti Space
Telescope Science Institute, Baltimore, MD 21218, USA  Hannah R. Wakeford
School of Physics, HH Wills Physics Laboratory, University of Bristol,
Bristol, UK

###### Abstract

We present a JWST/NIRCam transmission spectrum from 3.9–5.0 m of the recently-
validated sub-Earth GJ 341b ($\mathrm{R_{P}}=0.92$ R⊕, $\mathrm{T_{eq}}=540$
K) orbiting a nearby bright M1 star ($\mathrm{d}=10.4$ pc,
$\mathrm{K_{mag}}=5.6$). We use three independent pipelines to reduce the data
from the three JWST visits and perform several tests to check for the
significance of an atmosphere. Overall, our analysis does not uncover evidence
of an atmosphere. Our null hypothesis tests find that none of our pipelines’
transmission spectra can rule out a flat line, although there is weak evidence
for a Gaussian feature in two spectra from different pipelines (at 2.3 and 2.9
$\sigma$). However, the candidate features are seen at different wavelengths
(4.3 m vs 4.7 m), and our retrieval analysis finds that different gas species
can explain these features in the two reductions (CO2 at $3.1\sigma$ compared
to O3 at $2.9\sigma$), suggesting that they are not real astrophysical
signals. Our forward model analysis rules out a low mean molecular weight
atmosphere ($<350\times$ solar metallicity) to at least $3\sigma$, and
disfavors CH4-dominated atmospheres at $1-3\sigma$, depending on the
reduction. Instead, the forward models find our transmission spectra are
consistent with no atmosphere, a hazy atmosphere, or an atmosphere containing
a species that does not have prominent molecular bands across the NIRCam/F444W
bandpass, such as a water-dominated atmosphere. Our results demonstrate the
unequivocal need for two or more transit observations analyzed with multiple
reduction pipelines, alongside rigorous statistical tests, to determine the
robustness of molecular detections for small exoplanet atmospheres.

JWST, Terrestrial Exoplanet Atmospheres, Transmission Spectroscopy

††journal: AJ††facilities: JWST(NIRCam)††software:  Astropy (Astropy
Collaboration et al., 2013, 2018), batman (Kreidberg, 2015), CHIMERA (Line &
Yung, 2013; Line et al., 2014), Dynesty (Speagle, 2020), emcee (Foreman-Mackey
et al., 2013), Eureka! (Bell et al., 2022), ExoCTK (Bourque et al., 2021),
ExoTiC-LD (Grant & Wakeford, 2022) FIREFLy (Rustamkulov et al., 2022),
Forecaster (Chen & Kipping, 2017), IPython (Pérez & Granger, 2007), jwst
(Bushouse et al., 2022), Matplotlib (Hunter, 2007), NumPy (van der Walt et
al., 2011; Harris et al., 2020), PHOENIX (Allard et al., 2012) PICASO (Batalha
et al., 2019), PyMC3 (Salvatier et al., 2016), pysynphot(STScI Development
Team, 2013), SciPy (Virtanen et al., 2020), smarter (Lustig-Yaeger et al.,
2022, 2023a), Tiberius (Kirk et al., 2019, 2021)

##  1 Introduction

The M-dwarf opportunity is one of the most promising routes to measure the
atmospheres of terrestrial, and potentially habitable, exoplanets (e.g.
Charbonneau & Deming, 2007). This opportunity arises from the favorable
planet-to-star radius ratios, close-in habitable zones, and frequent transits
of rocky planets transiting M dwarfs. However, the M-dwarf opportunity must
contend with the high and prolonged X-ray and EUV emission from M-dwarf host
stars (e.g., Pizzolato et al., 2003; Preibisch & Feigelson, 2005; Shkolnik &
Barman, 2014; Peacock et al., 2019), which could cause significant exoplanet
atmospheric loss (e.g., Lammer et al., 2003; Owen & Jackson, 2012; Owen &
Mohanty, 2016).

Indeed, atmosphere observations of terrestrial planets (R${}_{\mathrm{P}}\leq
1.6$ R⊕) around M dwarfs have so far revealed no conclusive detections of
atmospheres (e.g., de Wit et al., 2016, 2018; Wakeford et al., 2019; Diamond-
Lowe et al., 2018, 2022; Kreidberg et al., 2019), with recent JWST studies of
TRAPPIST-1b and -c being consistent with no, or minimal, atmospheres (Greene
et al., 2023; Lim et al., 2023; Zieba et al., 2023).

In order to assess the M-dwarf opportunity, it is necessary to determine
which, if any, terrestrial planets around M dwarfs are capable of retaining an
atmosphere. One way to address this question is to determine whether M dwarf
exoplanets are separated by a ‘cosmic shoreline’, a division in
insolation–escape velocity space (I–vesc) separating bodies with atmospheres
from those without (Zahnle & Catling, 2017). In the Solar System, Zahnle &
Catling (2017) demonstrated that this dividing line follows an
$\mathrm{I}\propto\mathrm{v_{esc}}^{4}$ relation.

Following the successful launch and commissioning of JWST (Gardner et al.,
2023; McElwain et al., 2023; Rigby et al., 2023), we are now able to expand
the I–vesc parameter space over which we can address the existence of a cosmic
shoreline to terrestrial planets around M-dwarf stars. This is the focus of
JWST program 1981 (PIs: K. Stevenson and J. Lustig-Yaeger), which is obtaining
JWST transmission spectra of five exoplanets around M dwarfs to determine
whether a cosmic shoreline exists for M-dwarf systems.

Our first study (Lustig-Yaeger & Fu et al. 2023b) ruled out Earth-like,
hydrogen/helium, water- and methane-dominated clear atmospheres for the Earth-
sized exoplanet LHS 475b. This is consistent with predictions from the cosmic
shoreline. This study also confirmed LHS 475b as a bona fide exoplanet,
promoting it from a planet candidate in parallel with the efforts of Ment et
al. (2023).

Our second study (Moran & Stevenson et al. 2023) was of the super-Earth GJ
486b which sits close to the shoreline predicted by Zahnle & Catling (2017)
and thus provides an important test of its existence. We tentatively detected
water in GJ 486b’s transmission spectrum, however, our NIRSpec/G395H
observations could not determine the origin of the water, be it in the
planet’s atmosphere or the stellar atmosphere.

Most recently, in May & MacDonald et al. 2023, we observed two transits of the
super-Earth GJ 1132b, another planet near the shoreline. Our study revealed
the importance of multi-visit observations for small transiting exoplanets,
since one visit was suggestive of spectral features while the other produced a
flat transmission spectrum. Given the inconsistent transmission spectra
between our two visits, we were unable to make a conclusive determination of
the planet’s likely atmospheric composition, although neither visit was
consistent with a previous report of HCN (Swain et al., 2021).

Here, we extend our program to GJ 341b (TOI 741.01), a planet that was
recently validated by DiTomasso et al. (2023) and which orbits a nearby bright
M1 star ($K=5.6$, 10.4 pc, Gaia Collaboration et al. 2016, 2023). The TESS-
derived planetary radius is $0.92\pm 0.05$ R⊕, making this a sub-Earth, with a
$3\sigma$ upper limit on the planet’s mass of 4.5 M⊕ (DiTomasso et al., 2023).
The mass–radius relation of Chen & Kipping (2017) gives a planet mass of
$0.72\pm 0.14$ M⊕. Assuming a Bond albedo of zero and uniform heat
redistribution, the planet’s equilibrium temperature is 540 K, while its
irradiation temperature is 760 K. Given the planet’s escape velocity of 10 km
s-1, the cosmic shoreline relation predicts that this planet has not retained
a substantial atmosphere. Our study seeks to test this prediction.

The paper is laid out as follows: in Section 2 we describe the observational
setup. In Section 3 we explain our data reduction and light curve fitting,
using three independent pipelines. We present our transmission spectra, flat
line rejection tests, atmospheric forward models and atmospheric retrievals in
Section 4. We discuss our results and present our conclusions in Section 5.

##  2 Observations

We observed three transits of GJ 341b across three visits with JWST’s Near
Infrared Camera (NIRCam, Rieke et al. 2023) as part of GO program 1981 (PIs:
K. Stevenson and J. Lustig-Yaeger). Unlike the other GO 1981 targets that were
observed with NIRSpec, NIRCam was selected for GJ 341 because this bright star
($\mathrm{K_{mag}}=5.6$) saturates NIRSpec in one group. The three visits
occurred on March 2, March 10, and April 17, 2023 UT, with each visit lasting
5.25 hours. We used the F444W filter, covering a wavelength range from 3.8–5.1
m at a native spectral resolution of $R\approx 1650$, and the SUBGRISM128
subarray, giving us 2048 pixel columns and 128 pixel rows. Given the
brightness of our target, we used 3 groups per integration and acquired 4652
integrations during each of our three visits. The effective integration time
was 3.38 s.

##  3 Data reduction

As with the previous studies in our program (Lustig-Yaeger & Fu et al. 2023b;
Moran & Stevenson et al. 2023; May & MacDonald et al. 2023), we used three
different pipelines to reduce our data independently: Eureka!, Tiberius, and
Tswift. We describe the approach taken with each pipeline in more detail below
and include a table in the appendix that summarizes each pipeline’s approach
(Table 6).

###  3.1 Eureka!

The Eureka! data reduction pipeline (Bell et al., 2022) is well established,
with published JWST analyses demonstrating reliable results from numerous
observations (The JWST Transiting Exoplanet Community Early Release Science
Team et al. 2022; Ahrer et al. 2023; Alderson et al. 2023; Rustamkulov et al.
2023; Lustig-Yaeger & Fu et al. 2023b; Moran & Stevenson et al. 2023). As
described below, we adopt similar methods for the analysis of GJ 341b.

Starting from the uncal.fits files, we masked a region within 30 pixels of the
measured trace from pixel column 640 to 2048. We then performed a double-
iteration, group-level background subtraction by first fitting and removing
the mean flux on a column-by-column basis and, second, fitting and removing
the mean flux on a row-by-row basis. Unlike NIRSpec and NIRISS, NIRCam
detectors exhibit 1/f noise along pixel rows, thus necessitating this
additional row-by-row step. This orientation of the 1/f noise is because
NIRCam’s rapid readout direction aligns with the wavelength dispersion
direction (Schlawin et al., 2020). In Stage 1, we skipped the jump detection
step because our observations only had 3 groups per integration, which makes
the jump step unreliable. We processed the rateints.fits files through the
regular jwst Stage 2 pipeline.

In Stage 3, Eureka! performed a full-frame outlier search with a 4$\sigma$
rejection threshold. We then computed a clean median frame and corrected the
spectrum’s curvature by rolling pixel columns to bring the trace to the center
of the detector. We skipped any additional background subtraction and applied
optimal spectral extraction (Horne, 1986) using an aperture that was within 3
pixels of the measured trace. We tested a range of aperture sizes and found
that a 3-pixel half-width minimized the white light curve residuals. To
calculate our wavelength solution, we used a polynomial from Flight Program
1076 (E. Schlawin, private communication).

We extracted the flux from 3.9–5.0 m, splitting the light into 55
spectroscopic channels, each 0.02 m in width. For each light curve, we applied
a box-car filter of length 50 integrations and flagged any 3.5$\sigma$
outliers. This procedure removed $<2$% of all integrations for a given light
curve.

We fitted all three white light curves simultaneously. We used batman
(Kreidberg, 2015) to fit the system parameters and fixed the quadratic limb-
darkening coefficients to those provided by ExoTiC-LD (Grant & Wakeford, 2022)
from 3D stellar models (Magic et al., 2015). For our stellar parameters, we
used $\mathrm{T_{eff}}=3745$ K and $\log g=4.71$ g cm-2 and assumed a solar
metallicity ([Fe/H] = 0) (ExoMAST). Across the three visits, we fitted for a
single planet-star radius ratio ($R_{P}/R_{*}$), orbital inclination ($i$),
semi-major axis ($a/R_{*}$), and transit time ($T_{0}$). The period (P) was
held fixed to 7.5768707 days (ExoMAST). For each visit, we fitted an
independent exponential times linear ramp model to account for light curve
systematics. We also renormalized the error bars from each visit in order to
achieve a reduced chi-squared of unity. In total, there were 19 free
parameters. Table 1 lists our best-fit system parameters.

When fitting the spectroscopic light curves, we fixed the orbital inclination,
semi-major axis, and transit time to the white light curve best-fit values.
Again, we fixed the quadratic limb-darkening parameters to the values provided
by ExoTiC-LD for each spectroscopic channel. The planet-to-star radius ratio
was the only physical free parameter. For each spectroscopic channel, we
fitted an independent exponential times linear ramp model to account for light
curve systematics. Including the error bar renormalization term, each
spectroscopic fit had six free parameters. For both the white and
spectroscopic light curves, we clipped the first 300 integrations from each
visit to avoid fitting the exponential decrease in flux. Figure 1 shows our
fits to the Eureka! white and spectroscopic light curves, with the Allan
variance plot of the residuals shown in Figure 2.

Reduction | $T_{0}(\mathrm{BJD_{TDB}})$ | $i(^{\circ})$ | $a/R_{*}$ | $R_{P}/R_{*}$
---|---|---|---|---
Eureka! | $2460006.420147^{+0.000152}_{-0.000141}$ | $88.17^{+0.64}_{-0.51}$ | $20.12^{+2.93}_{-2.19}$ | $0.016079^{+0.000271}_{-0.000254}$
Tiberius | $2460006.420152^{+0.000193}_{-0.000196}$ | $88.77^{+0.88}_{-0.97}$ | $22.88^{+2.97}_{-4.37}$ | $0.015315^{+0.000316}_{-0.000297}$
Tswift | $2460006.420143^{+0.000241}_{-0.000245}$ | $88.72^{+0.41}_{-0.46}$ | $23.48^{+1.44}_{-1.95}$ | $0.016130^{+0.000295}_{-0.000337}$
DiTomasso et al. | $2459301.771\pm 0.002$ | $89.22^{+0.54}_{-0.96}$ | $24.50^{+2.07}_{-4.01}$ | $0.0161^{+0.0007}_{-0.0008}$
Table 1: GJ 341b’s system parameters from the three independent reductions’
white light curve fits. We also include the parameters from the independent
TESS and radial velocity analysis (DiTomasso et al., 2023). Figure 1: The
light curves and best fit models resulting from the Eureka! pipeline. The
columns show, from left to right: transit 1, transit 2, transit 3. The rows
show, from top to bottom: 1) the white light curve and residuals, 2) the
spectroscopic light curves, 3) the best-fit models to the spectroscopic light
curves, 4) the residuals from the spectroscopic light curve fitting. Figure 2:
The Allan variance plot from the white light curve fits to each pipeline’s
reduction and for all three visits. The amplitude of the residuals from
Tiberius (blue), Eureka! (orange) and Tswift (green) are shown as a function
of bin size in minutes. In the absence of red noise, the residuals would
follow the dashed black line.

###  3.2 Tiberius

The Tiberius reduction began with the uncalibrated fits files (uncal.fits). We
ran these through the stage 1 detector processing stages of the jwst pipeline.
However, we did not perform the jump step as this has been found to increase
the noise in time series spectrophotometry (e.g., Rustamkulov et al., 2023).
Furthermore, unlike in the application of Tiberius to NIRSpec data (e.g.,
Rustamkulov et al. 2023; Lustig-Yaeger & Fu et al. 2023b; Moran & Stevenson et
al. 2023) and in our Eureka! reductions of the data presented here, we did not
apply a 1/f correction to the data. This is because the 1/f noise is parallel
to the dispersion dimension for NIRCam data, complicating its removal.

Following the stage 1 detector processing, we worked with the integration
level images (gainscalestep.fits files). While these were processed through
the gain_scale step, the jwst pipeline does not apply a gain correction if the
default gain setting is used. For this reason, we worked in units of DN/s
throughout the rest of our analysis. Since we work with normalized light
curves, and rescale our photometric uncertainties, the gain correction is not
necessary.

Before extracting the stellar flux from each image, each pixel was compared to
the median value for that pixel across the time series and $>5\sigma$ outliers
were replaced by the median for that pixel. We then oversampled each pixel by
a factor of 10 in the spatial dimension using a flux-conserving interpolation,
which allows us to more precisely place the aperture boundaries. We found this
step to reduce white noise by 10 % in previous tests on JWST data (Rustamkulov
et al., 2023). Next, we traced the stellar spectrum using a fourth order
polynomial between columns 700 and 2042 (inclusive, zero-indexed) and
performed standard aperture extraction. We subtracted the background at every
column using the median value across the full column, after masking a total of
45 pixels centered on the stellar trace, and summed up the stellar flux in an
aperture with a full width of 5 pixels. We experimented with wider apertures
and found these led to consistent transmission spectra and light curves but
with greater noise.

We fitted our white light curves from the three visits simultaneously after
subtracting one orbital period (7.576863 days, DiTomasso et al. 2023) from
transit 2’s time array and six orbital periods from transit 3’s time array so
that all time arrays could be fitted with a single $T_{0}$, which corresponded
to transit 1. The 16 free parameters in our white light curve fits were
$T_{0}$, $i$, $a/R_{*}$, $R_{P}/R_{*}$ and the 12 parameters defining our
systematics model. We defined this model as a single component exponential
ramp multiplied by a linear-in-time polynomial for each visit. The period was
held fixed to 7.576863 days (DiTomasso et al., 2023). We used batman
(Kreidberg, 2015) to generate our analytic transit light curves.

For our limb darkening calculations, we adopted $\mathrm{T_{eff}}=3770\pm 40$
K and $\log g=4.72\pm 0.02$ g cm-2 (DiTomasso et al., 2023) and assumed a
solar metallicity ([Fe/H] = 0). We used a quadratic limb darkening law with
coefficients held fixed to values calculated from 3D stellar atmosphere models
(Magic et al., 2015) using ExoTiC-LD (Grant & Wakeford, 2022) and our own
custom instrument throughput. This was constructed by dividing our observed
stellar spectrum by the appropriate Phoenix model spectrum (Husser et al.,
2013) and smoothing the result with a running median, following the approach
of Alderson et al. (2023).

We explored the parameter space with a Markov chain Monte Carlo (MCMC) using
the emcee Python package (Foreman-Mackey et al., 2013). We ran three sets of
chains each with 50,000 steps and 200 walkers. After the first run, we
discarded the first 25,000 steps and calculated the median values from the
second 25,000 steps of our chains which were used as starting points for the
second run. After the second run, we again discarded the first 25,000 steps
and used the best-fitting model to rescale our photometric uncertainties to
give $\chi^{2}_{\nu}=1$. Our third and final run used the rescaled photometric
uncertainties. The medians, 16th and 84th percentiles of our posteriors from
this third run are reported in Table 1. These posteriors were constructed
after discarding the first 100 steps and thinning by a factor of 10.

Our spectroscopic light curves were created using 57 20-pixel-wide wavelength
bins between 3.906 and 5.004 m and we adopted the same wavelength solution as
used in the Eureka! reduction. For our spectroscopic light curve fits, we
fixed the system parameters ($T_{0}$, $i$ and $a/R_{*}$) to the Tiberius
values in Table 1 and fitted each visit independently, not jointly. Prior to
fitting, we divided each spectroscopic light curve by a visit-dependent common
noise model. This noise model was defined as WLCv/tmv where WLCv is the visit-
dependent white light curve and tmv is the visit-dependent best-fitting
transit (batman) model. For our spectroscopic light curves, the exponential
coefficients did not converge for every wavelength bin across each visit,
despite experimenting with different priors. For this reason, at the
spectroscopic light curve stage, we clipped the first 600 integrations (40
minutes) from each visit, after removing the common noise, and fitted the
light curves with a linear polynomial and no exponential. This gave comparable
residual red noise to the fits where the exponential coefficients did
converge, indicating that the ramp was adequately removed by clipping the
first 600 integrations. Figure 11 in the Appendix shows the fits to our white
and spectroscopic light curves, with the Allan variance plot of the residual
red noise from our white light curve fits shown in Figure 2.

###  3.3 Tswift

Transit swift (Tswift) is a JWST data reduction routine developed to swiftly
analyze time-series observation of transiting exoplanet atmospheres. Starting
from the uncal.fits images, we used the default jwst pipeline with steps.jump
turned off to produce stage one rateints.fits images. Next, we performed row-
by-row and column-by-column background subtraction (Schlawin et al., 2020) for
each integration by subtracting out the median of the non-illuminated regions
from each row. Next, we cross-correlated the spectral line spread function
with each column to determine the vertical position of the spectral trace for
that column. We then used a third-order polynomial to fit the vertical
position versus the column number, determining the spectral trace location for
each integration. The spectrum was extracted using a full-width 6-pixel wide
aperture centered on this polynomial fit as it minimizes the scatter in the
white light curve. Given the brief exposure time for each integration and the
high number of integrations per visit, cosmic rays or hot pixels have a
minimal impact on the transit light curve. Furthermore, we identified and
discarded any outliers in time exceeding 5$\sigma$ based on the transit light
curve baseline scatter for each column. Finally, we combined all wavelength
channels to produce the white light curve of each transit.

The fitting of each individual white light curve employs a combination of
batman (Kreidberg, 2015) and emcee (Foreman-Mackey et al., 2013). This
approach incorporates eight parameters: $T_{0}$, $R_{P}/R_{*}$, $a/R_{*}$,
$i$, a linear time slope, a constant scaling factor, and two coefficients for
the exponential ramp. The averaged best-fit parameters from the three white
light curves are listed in Table 1. To determine the quadratic limb darkening
parameters, we used the 3D stellar models (Magic et al., 2015) from the
ExoTiC-LD (Grant & Wakeford, 2022) package and subsequently fixed both
quadratic terms, using the same stellar parameters as in the Eureka! and
Tiberius analyses. We tested fitting for the second quadratic term while
fixing the first one, and it had no effect on the final transmission spectrum.

The NIRCam detector has four amplifiers and the spectrum is dispersed across
the second to fourth amplifier in the F444W observation mode (Schlawin et al.,
2020). As already discussed, there is 1/f noise along the dispersion direction
in our NIRCam data. In an ideal scenario, this noise could be subtracted if
every row in each amplifier had unilluminated pixels. However, due to the
continuous dispersion of spectra across various amplifiers, it is not feasible
to correct the entire spectrum. As a result, we observed lingering 1/f noise
residuals that showed correlation across wavelengths in the light curve fits.
We identified that these residual patterns are predominantly consistent across
all amplifiers, with a few minor patterns unique to individual amplifiers. To
better remove this correlated noise, we performed a white light fit within
each amplifier, and the residuals were then used as a common-mode correction
during the spectroscopic light curve fits.

For our spectroscopic light curve fits, we used batman and
scipy.optimize.curvefit. The model consists of five parameters: $R_{P}/R_{*}$,
a constant scaling factor, two coefficients for the exponential ramp, and a
scaling factor for the common-mode correction. We held $T_{0}$, $a/R_{*}$, and
$i$ fixed to the corresponding averaged best-fit values from the white light
curves (Table 1). Figure 12 shows our fits to the Tswift white and
spectroscopic light curves, with the Allan variance plot of the residuals
shown in Figure 2.

###  3.4 Engineering parameters

As can be seen in the Allan variance plots from our reductions (Figure 2)
there is residual red noise which becomes apparent at time scales longer than
1 minute. In order to remove this red noise, we explored detrending our light
curves with engineering telemetry data timeseries from the JWST Engineering
Database. The telemetry data are identified by mnemonics, which follow the
naming convention X_Y_Z, where X corresponds to the system (I for Integrated
Science Instrument Module, S for spacecraft), Y corresponds to the subsystem,
instrument, controller, or component (1–4 letters), and Z corresponds to the
telemetry name (2–20 characters). Throughout the rest of this section, we
refer to these timeseries telemetry data sets as ‘mnemonics’.

We began by whittling down the 6176 NIRCam-related mnemonics to the 596
mnemonics of time-varying floats (as opposed to booleans and non-varying
constants). We then resampled each mnemonic onto the same time-sampling as
each of our three visits and proceeded to sequentially incorporate each of the
596 mnemnonics into our white light curve systematics model, for the Tiberius
reduction, to test whether the addition of a mnemnonics model decreased the
residual red noise.

For these fits we used a Levenberg-Marquadt alogorithm, with the systematics
model being the summation of a linear polynomial describing the mnemonic and a
linear-in-time polynomial, after clipping the first 600 integrations. We then
compared the residual red noise in each case to our fiducial model: a linear-
in-time polynomial. We defined the residual red noise as being the RMS of the
residuals in bins of 200 integrations divided by the white noise expectation
(equivalent to dividing the color lines in our Allan variance plots by the
dashed black lines at $x=14.5$ mins, Figure 2). Of the 596 possible mnemonics,
we found that 112 mnemonics led to a decrease in the residual red noise
compared to the fiducial model for all three visits. The ten mnemonics that
led to the largest decrease in the RMS are given in Table 2. As shown in Table
2, these mnemonics mostly correspond to temperature.

Top 10 mnemonics (all visits) | Short description | Residual red noise
---|---|---
INRC_FA_PSC_LVCBRD_TEMP | LVC_BOARD_TEMP | 4.79
INRC_FA_PSC_SCRBRD_TEMP | SCR_BOARD_TEMP | 4.90
INRC_FB_PSC_SCRBRD_TEMP | SCR_BOARD_TEMP | 5.16
INRC_FA_ACE2_SC_TEMP_HI | SC_ASIC temperature high gain | 5.20
IGDP_NRC_FA_ACE2_SCTEMP | ASIC Temperature read by ACE card | 5.20
INRC_FB_PSC_LVCBRD_TEMP | LVC_BOARD_TEMP | 5.23
INRC_FA_ACE1_RR_RD | Voltage of ACE reference reference | 5.37
IGDP_NRC_FB_ACE1_SCTEMP | ASIC Temperature read by ACE card | 5.39
INRC_FB_ACE1_SC_TEMP_HI | SC_ASIC temperature high gain | 5.39
INRC_FB_ACE1_D2P5_V | Voltage of ACE Digital 2.5V voltage | 5.39
Fiducial | Linear-in-time polynomial | 5.50
Table 2: The ten engineering mnemonics that led to the largest improvement in
the cumulative residual red noise from all three white light curve fits. These
are ordered from largest improvement (top) to least improvement (bottom). With
regards to the mnemonics, ‘INRC’ is NIRCam telemetry reported by ISIM, the
Integrated Science Instrument Module, ‘IGDP’ are ‘ground data points’ computed
from downlinked engineering telemetry, and ‘TEMP’ corresponds to temperature.
We understand that ‘A’ and ‘B’ correspond to module A and module B, ‘ACE’
corresponds to Application Specific Integrated Circuit (ASIC) control
electronics, and ‘LV’ corresponds to low voltage. We define the residual red
noise in the text.

Following this test, we ran a second test where we increased the number of
mnemonics used to fit the white light curves, going from one to ten mnemonics
in the order presented in Table 2. As can be seen in Figure 3, including
additional mnemonics does not necessarily lead to an improvement in the
residual red noise. This is likely due to some mnemonics representing similar
telemetry data (e.g., INRC_FA_PSC_LVCBRD_TEMP and INRC_FA_PSC_SCRBRD_TEMP).
This figure also demonstrates that the use of a mnemonic significantly
improves the residual red noise in visit 2 while the changes in visits 1 and 3
are more subtle. Additionally, a combination of 10 mnemonics does not
eradicate the residual red noise. As a further test, we fitted all 596
mnemonics jointly with a linear-in-time polynomial. In this case, the residual
red noise was below the white noise expectations. We interpreted this as
evidence that including a large number of mnemonics can lead to over-fitting.

While independently fitting each visit’s light curve with mnemonics
demonstrated the residual red noise could be improved upon, we ended up not
using the mnemonic data in our final white light curve fits. This was because
we could not obtain acceptable fits when we tried to fit the three visits
simultaneously with mnemonics, likely due to the fact that the amplitude of
the systematics in our data is comparable to the transit depth. We advocate
for a further exploration of mnemonics for exoplanets with deeper transits
where the transit and systematics models are less degenerate.

Figure 3: The improvement in residual red noise in the white light curves by
adding an additional mnemonic / telemetry timeseries to the systematics model.
The blue, orange and green lines indicate the residuals from visit 1, 2 and
3’s white light curves. The fiducial model is a linear-in-time polynomial.
Each x-tick label indicates an additional mnemonic is added. The dashed gray
line indicates the white noise level that the residual noise would achieve if
no red noise is present.

###  3.5 Short wavelength photometry

Since our observations were conducted with NIRCam, we additionally have access
to the short wavelength (SW) channel between 0.6–2.3 m. This channel acquired
photometry simultaneously with the long wavelength spectroscopy described
above. We extracted the short wavelength photometry using Tiberius to
determine whether we could use these additional white light curves to improve
the precision of our system parameters.

As with the long wavelength data, we began by processing the uncal.fits files
through stage 1 of the jwst pipeline. We then extracted the time-series
photometry using Tiberius with an aperture 100 pixels wide and a center fixed
at a pixel column of 1795. We also subtracted the background at every pixel
row, which was calculated as the median of two 50-pixel-wide regions separated
by 100 pixels from each side of the target aperture. The resulting SW light
curves for each visit are shown in Figure 4. Given the significant noise in
these light curves, we did not use these in our parameter estimation.

Figure 4: The short wavelength photometric light curves for each visit,
revealing significant pre-transit noise (visit 2) and in-transit noise (visit
3). These have been smoothed with a running median calculated over a sliding
box with a width of 51 integrations. The dashed vertical lines indicate the
first and fourth contact points of the transit.

###  3.6 Measuring flux-calibrated spectra

Following our study of GJ 486b (Moran & Stevenson et al. 2023), we extracted
the flux calibrated stellar spectrum of GJ 341. We did this in an attempt to
constrain the covering fractions of star spots and faculae and hence determine
the impact of stellar contamination on our transmission spectra. To do this,
we processed our stellar spectra through stage 2 of the jwst pipeline.
Specifically, taking the gainscalestep.fits files, we processed these through
the assign_wcs, flat_field, extract_2d, srctype, pathloss and photom steps. We
did not perform the background step, and instead used an alternate version of
our stage 1 reduction where we did perform a group-level background
subtraction.

We then used Tiberius to perform standard aperture photometry on our
photomstep.fits files using an aperture with a full width of 20 pixels.
Despite running the srctype command with the \--source_type POINT argument
defined, the resulting units of our spectra were MJy/sr. We therefore had to
multiply our spectra by the solid angle of a pixel, listed as PIXAR_SR in the
photomstep.fits header, to convert to MJy. We defined the flux uncertainties
in the flux calibrated spectra as being the standard deviation of each pixel’s
flux across the time series. We found no significant variability in the flux-
calibrated spectra from visit-to-visit ($\leq 0.1$ %) and within each visit
($0.24$ %).

We then computed multi-component forward models using the Allard et al. (2012)
PHOENIX stellar models to quantify star spot and faculae covering fractions of
the stellar surface for each of the three visits. We employed a weighted
linear combination of three PHOENIX models to represent the background
photosphere, spots ($T_{\rm eff}$ $\leq$ $T_{\rm eff,photosphere}$ \- 100 K),
and faculae ($T_{\rm eff}$ $\geq$ $T_{\rm eff,photosphere}$ \+ 100 K). We
assumed that all spots have the same $T_{\rm eff}$, log($g$), and metallicity,
as do the faculae, and each feature is required to not exceed 45% of the total
surface. The grid of PHOENIX models used for our analysis covered $T_{\rm
eff}$ = 2800 – 4000 K, log($g$) = 4 – 5.5 cm s-2, [Fe/H] = -0.5 – 0.5
providing extensive coverage of possible spot and faculae temperatures for GJ
341 assuming photospheric values near the literature quoted $T_{\rm eff}$=
3770 $\pm$ 40 K and log($g$) = 4.72 $\pm$ 0.02 cm s-2 (DiTomasso et al.,
2023). The models were scaled by R${}_{*}^{2}$/dist2 using literature values
for GJ 341: R∗=0.506 R⊙ (DiTomasso et al., 2023) and $d$ = 10.45 pc (Gaia
Collaboration et al., 2021). We also smoothed and interpolated the models to
be the same resolution as the observations before calculating a
reduced-$\chi^{2}$. In our reduced-$\chi^{2}$ calculations, we considered 868
wavelength points for each visit and eight fitted parameters (the $T_{\rm
eff}$, log($g$), and [Fe/H] of the photosphere, the $T_{\rm eff}$ and coverage
fraction of both spots and faculae, and a scaling factor). The scaling factor
was multiplied by the R${}_{*}^{2}$/$d$2 term to account for uncertainty in
either measured quantity and varied from 0.9 to 1.3.

Unlike in our application of this procedure to the stellar spectra of GJ 486
(Lustig-Yaeger & Fu et al. 2023b) and GJ 1132 (May & MacDonald et al. 2023),
which used NIRSpec/G395H, we could not obtain an acceptable fit to the flux-
calibrated NIRCam/F444W spectra ($\chi^{2}_{\nu}\sim 60$), as shown in Figure
5. We do not know the reason for these poor fits, however, the shape of our
extracted spectrum changes significantly when we apply the photometric
calibration (photom) step in the jwst pipeline, suggesting that the NIRCam
PHOTOM reference file may be responsible. Given the poor fits in Figure 5, and
lack of contamination evident in our transmission spectra (see section 4), we
do interpret the results from our analysis of the flux calibrated spectrum as
evidence for stellar contamination in our data.

Figure 5: A comparison between our measured flux calibrated spectra (black)
and superpositions of PHOENIX stellar atmosphere models: photosphere with no
active region contribution (purple), photosphere with contribution from spots
(orange), and photosphere with contributions from spots and faculae (green).
Given the poor $\chi^{2}_{\nu}$, and lack of contamination evident in our
transmission spectra, we do interpret these model results as evidence for
stellar contamination in our data.

##  4 GJ 341b’s transmission spectrum

The transmission spectra we measure for GJ 341b from the three different
reductions are shown in Figure 6. This figure demonstrates that the three
independent reductions produce consistent transmission spectra. For this
figure, we offset the Eureka! spectrum by -15 ppm so that its median transit
depth matches that of Tiberius and Tswift between 4.0 and 4.8 m. We include
the transmission spectra measured for each visit by each reduction in the
Appendix (Figure 13) which demonstrates that the spectra are also consistent
between visits for each reduction. In the following subsections, we further
investigate the transmission spectrum of GJ 341b with flat line rejection
tests, atmosphere forward models, and retrievals.

Figure 6: GJ 341b’s transmission spectrum as measured by Tiberius (blue),
Eureka! (orange) and Tswift (green). These spectra are the the weighted means
of the three visits’ spectra. The Eureka! spectrum has been offset by -15 ppm
to make the median transit depths equal to Tiberius and Tswift between 4.0 and
4.8 m.

###  4.1 Gaussian flat line rejection tests

We begin our assessment of the GJ 341b transmission spectrum by quantifying
the statistical significance of any potential features in the spectrum
relative to the null hypothesis. In this case, the null hypothesis is that the
planet possesses a flat, featureless transmission spectrum that is devoid of
atmospheric information, as would be the case if the planet does not possess
an atmosphere or the atmospheric signal falls below our detection sensitivity
(Kreidberg et al., 2014; Lustig-Yaeger & Fu et al. 2023b).

Following the approaches of Moran & Stevenson et al. (2023) and May &
MacDonald et al. (2023), we perform two statistical tests on the transmission
spectra from each data reduction pipeline. The first test quantifies the
(in)consistency of our measurements with the null hypothesis given the
observational errors, while the second test quantifies the evidence favoring
an agnostic Gaussian absorption feature in the spectrum over the flat null
hypothesis. As both tests use slightly different approaches to quantify the
significance of one or more features in the spectrum, they are both useful
diagnostics. However, these tests are limited to statistical statements about
the spectrum as opposed to physical explanations, which require the more
detailed modeling presented in the forthcoming subsections.

To understand whether the observed spectrum is consistent with the null
hypothesis, we generated 100,000 random spectra on the same wavelength grid
subject to the assumption that the underlying true planet spectrum is the null
hypothesis and the observational errors are unchanged from those observed. We
calculated the reduced chi-squared ($\chi^{2}_{\nu}$) for each sample and
build up the expected distribution of $\chi^{2}_{\nu}$ under the null
hypothesis. That is, if the null hypothesis were true then we would expect
that our measured $\chi^{2}_{\nu}$ would be a single sample from that
distribution. Thus, we can quantify the probability that the measured
$\chi^{2}_{\nu}$ value would be at least as large as the one we observed given
the assumption of the null hypothesis (i.e., a p-value).

To understand whether there is evidence of a spectral feature in the observed
spectrum, we first fit a flat featureless spectrum to the data and then fit an
agnostic Gaussian model to the spectrum. We then compare the Bayesian evidence
of the two fits and calculate Bayes factors to estimate which model is
preferred, as in Moran & Stevenson et al. (2023) and May & MacDonald et al.
(2023). The Gaussian model contains four fitting parameters, which are all
constrained by uniform priors. These include the wavelength independent
transit depth upon which the Gaussian feature is added
${\sim}\mathcal{U}(0.0,1.0)$, the wavelength of the Gaussian center
${\sim}\mathcal{U}(3.5,5.5)$ m, the width of the feature
${\sim}\mathcal{U}(\Delta x_{0},\Delta x_{1})$ m, where $\Delta x_{0}$ is the
minimum spectroscopic bin width and $\Delta x_{1}$ is the width of the entire
wavelength range of the data, and the amplitude of the Gaussian feature
${\sim}\mathcal{U}(-10\sigma_{y},10\sigma_{y})$, where $\sigma_{y}$ is the
standard deviation of the transmission spectrum measurements (i.e., features
should not greatly exceed the scatter in the data).

The results of these tests are shown in Table 3. This demonstrates that the
Eureka! spectrum is the flattest of the three spectra, ruling out our null
hypothesis at only $0.2\sigma$, while the Gaussian feature test favors a flat
line over a feature at $2.2\sigma$. The Tiberius and Tswift spectra are less
flat. However, in the case of the Tiberius spectrum we believe this is due to
the smaller uncertainties in its transmission spectrum (see Section 4.2). For
these reductions, the null hypothesis is in weak tension, at $1.3\sigma$ for
Tiberius and $1.4\sigma$ for Tswift. The Gaussian feature test results in
$2.3\sigma$ evidence for a feature in the Tiberius spectrum (broad and
centered near 4.3 m) and $2.9\sigma$ in the Tswift spectrum (width of
${\sim}0.2$ m centered near 4.7 m). While these look promising at face value,
as we see from these tests and demonstrate with our subsequent retrieval
analysis (section 4.4), the two candidate features occur at different
wavelengths between the Tiberius and Tswift spectra and are therefore less
likely to be astrophysical.

Table 3: Is it Flat?

Reduction | $\chi^{2}_{\nu}$ | p-value ($\sigma$) | Feature?
---|---|---|---
Tiberius | 1.15 | 20% (1.3$\sigma$) | 2.3$\sigma$
Eureka! | 0.82 | 83% (0.2$\sigma$) | -2.2$\sigma$
Tswift | 1.18 | 17% (1.4$\sigma$) | 2.9$\sigma$

Note. — $\chi^{2}_{\nu}$ is the reduced chi-squared resulting from the best
fitting featureless fit (null hypothesis) to the observed spectrum, “p-value”
refers to the probability that $\chi^{2}_{\nu}$ would be at least as extreme
as the observed value under the assumption that the null hypothesis is correct
and is displayed along with the corresponding “sigma” value, and the
“Feature?” column shows the level of confidence in the detection of an
agnostic Gaussian absorption feature in the spectrum over the null hypothesis.

###  4.2 Sensitivity to Spectrum Uncertainties

The results of the aforementioned flat line rejection tests depend not just on
potential features in the spectrum, but also on the size of the error bars,
which are not equivalent for each reduction. The 1$\sigma$ spectroscopic
transit depth uncertainties are shown in Figure 7 (lower panel). Since all
reductions used approximately the same spectral bin widths, we can readily
compare their 1$\sigma$ uncertainties and apply the errors from one reduction
to the others to estimate the role of uncertainties in driving the statistical
results of our null hypothesis tests. Table 4 shows $\chi^{2}_{\nu}$ values
for flat line fits to each reduction’s transmission spectrum data points
(rows) with errors interpolated from the other reductions (columns).

Table 4: Reduced-$\chi^{2}$ Values for a flat line fit when exchanging errors between reductions  |  Tiberius Errors |  Eureka! Errors |  Tswift Errors
---|---|---|---
Tiberius | 1.15 | 0.68 | 0.79
Eureka! | 1.28 | 0.82 | 0.75
Tswift | 1.81 | 1.10 | 1.18
Figure 7: Light curve residual RMS (upper panel) and transit depth
uncertainties (lower panel) for each reduction (color lines). The expected
noise (photon + readnoise + dark current) calculated from the data is shown in
the upper panel (black line) and the median RMS relative to the expected noise
is listed in the legend for each reduction. The predicted noise from PandExo
(Batalha et al., 2017) is shown in the lower panel (black line).

Figure 7 (lower panel) and Table 4 show that the differences in uncertainties
between reductions leads to statistically significant differences in
$\chi^{2}_{\nu}$. The Tiberius reduction has the smallest uncertainties and
when those are applied to the Eureka! and Tswift spectra, those spectra
deviate more from a flat line. On the other hand, the Tswift spectrum has
$\chi^{2}_{\nu}>1$ regardless of which errors are applied, consistent with the
larger deviations seen in the Tswift spectrum. The Eureka! spectrum has
relatively large errors, but relatively small spectroscopic deviations, which
together explain our statistical result and subsequent physical
interpretations.

The upper panel of Figure 7 shows the light curve residual RMS and the
expected noise for comparison, which includes the photon noise, read noise and
dark current. The Tiberius, Tswift, and Eureka! reductions achieve a median
precision of 1.03$\times$, 1.05$\times$, and 1.14$\times$ the noise limit,
respectively. This ordering is different to the white light curve residuals,
where Eureka! obtained the lowest residual RMS (Figure 2). The lower RMS of
Tiberius and Tswift at the spectroscopic light curve stage comes from those
pipelines’ use of a common mode correction. Although the Tiberius transit
depth uncertainties are lower than the predicted PandExo uncertainties, this
is likely due to improved performance relative to the model-based prediction,
since the Tiberius reduction still sits above the noise limit (Figure 7, upper
panel).

###  4.3 Atmosphere forward models

As with previous targets in program JWST-GO-1981, we compare our transmission
spectra of GJ 341b to a set of atmospheric forward models to infer whether the
data suggests that the planet possesses an appreciable atmosphere, and if so,
its composition (Lustig-Yaeger & Fu et al. 2023b; Moran & Stevenson et al.
2023; May & MacDonald et al. 2023). Since GJ 341b has no published mass
constraints beyond an upper limit ($\leq$ 3 M⊕ at 1$\sigma$; DiTomasso et al.,
2023), we generate forward models using a planetary mass of 0.73 M⊕. This mass
assumes that the planetary radius of 0.92 R⊕ is consistent with a rocky bulk
composition, as Luque & Pallé (2022) showed that very small ($\leq$ 1 R⊕)
planets around M-dwarfs, like GJ 341b, appear to have very tightly constrained
densities similar to Earth-like iron-silicate bulk compositions.

Our atmospheric forward models include two cases: 1) CHIMERA thermochemical
equilibrium models (Line et al., 2013, 2014) from 100-1000$\times$ solar
metallicity, using a parameterized temperature-pressure profile (Guillot,
2010) with an equilibrium temperature of 540 K and opacity contributions from
H2O, CH4, CO, CO2, NH3, HCN, H2S, H2, and He; and 2) PICASO (Batalha et al.,
2019) models of single-gas atmospheres with an isothermal temperature-pressure
profile at 540 K. The PICASO models include atmospheres of pure 1 bar H2O,
CO2, or CH4, as well as a flat line model representative of no atmosphere or a
uniform aerosol layer. For all model configurations, we compute synthetic
transmission spectra with PICASO’s radiative transfer routine using molecular
opacities resampled from R = 60,000 (Batalha et al., 2020), rebinned to the
resolution of each data reduction. We compute best fits by calculating the
$\chi^{2}_{\nu}$ between each reduction and each forward model, with 57, 54,
or 55 degrees of freedom (dof) for the Tiberius, Tswift, and Eureka!
reductions, respectively. We include an offset in relative transit depth to
the mean of each reduction in our fitting. Our results for Eureka! are shown
in Figure 8, while those for Tiberius and Tswift are shown in Figures 14 and
15. We include the complete set of statistical constraints for each reduction
compared to each forward model case in Table 5.

Figure 8: Atmospheric forward models compared to the Eureka! reduction from
JWST/NIRCam. Because of the planet’s small size and low equilibrium
temperature, only hydrogen-dominated or methane atmospheres are strongly ruled
out or disfavored by the data. Atmospheres of water, carbon dioxide, or an
opaque haze layer – or no atmosphere – are all statistically consistent with
the data to within $\sim$1$\sigma$ across all three reductions.

Across all three reductions, hydrogen-dominated or methane atmospheric models
do not fit the data. For the Tswift reduction, hydrogen-dominated atmospheres
as metal-rich as 1000$\times$ solar can be ruled out to 3$\sigma$, while this
constraint is reduced to 350$\times$ solar for Tiberius and Eureka!. Such
hydrogen-rich atmospheres would be unexpected to persist against escape, so
this result is unsurprising. Methane atmospheres of 1 bar are disfavored at
3$\sigma$, 2.5$\sigma$, or 1$\sigma$ between Tswift, Tiberius, and Eureka!.

On the other hand, water atmospheres or a flat line – which could stem either
from an airless world or from a hazy world – are consistent with each
reduction, with $\chi^{2}_{\nu}$ $\lesssim$ 1\. Though water is a ubiquitous
molecule and common in warm planetary atmospheres, we are unable to rule it
out here since water has no strong absorption bandheads and instead acts as a
continuum-like absorber from 4-5 $\mu$m. Both Eureka! and Tiberius also are
consistent ($\chi^{2}_{\nu}$ $\lesssim$ 1) with a carbon dioxide atmosphere,
while Tswift disfavors this scenario to 2$\sigma$.

Despite the excellent agreement between each data reduction (Figure 6), there
are obvious differences in goodness of fit between atmospheric models. This
occurs due to 1) the planet’s small size and temperature, which results in a
very small atmospheric scale height, and 2) the ratio of the sub-Earth radius
(RP=0.92 R⊕) to the relatively large radius of the star for an M-dwarf
(R∗=0.506 R⊙), resulting in a very shallow transit. Given these physical
parameters, our 350$\times$ solar, hydrogen-dominated atmosphere model has a
peak feature size of 40 ppm while a carbon dioxide atmosphere model has a peak
feature size of only 15 ppm in the NIRCam/F444W bandpass. Even with the
extremely high precision obtained by our program (average uncertainties of 15
ppm per 0.02 $\mu$m bin), such small feature sizes mean that minute
differences between reductions driven by random noise or data treatments drive
divergent atmospheric interpretations. Overall, either no atmosphere, a hazy
atmosphere, or an atmosphere containing a species that does not have prominent
molecular bands across the NIRCam/F444W bandpass – e.g., a water-dominated
atmosphere – provide the best explanation of the data across the three
reductions compared to our forward models.

Table 5: Atmospheric Constraints

Model | Tiberius | Eureka! | Tswift | Interpretation
---|---|---|---|---
350$\times$ solar | 1.78 | 1.71 | 2.91 | rejected at $\geq$ 3$\sigma$
H2O | 1.18 | 0.79 | 1.11 |  $\sim$ consistent with data
CO2 | 1.02 | 0.88 | 1.41 | inconclusive
CH4 | 1.48 | 1.10 | 1.70 | ruled out to $\geq$1$\sigma$
Flat line | 1.15 | 0.82 | 1.22 |  $\sim$ consistent with data

Note. — The $\chi^{2}_{\nu}$ of each model is reported in the columns for each
reduction name. All atmospheric models are computed for 1 bar and an
equilibrium temperature of 540 K. The “flat line” model can indicate no
atmosphere or an opaque aerosol layer.

###  4.4 Atmosphere retrievals

To examine a broader range of atmospheric parameter space that is
(in)consistent with our measurements, we ran a set of retrievals using the
Spectral Mapping Atmospheric Radiative Transfer Exoplanet Retrieval code
(Lustig-Yaeger et al., 2022, 2023a, smarter;; Lustig-Yaeger & Fu et al.
2023b). We use the same smarter modeling approach and assumptions as described
in Lustig-Yaeger & Fu et al. (2023b), except slightly modified to account for
the different planet parameters. We run separate retrievals on each of the
final coadded spectra from the Eureka!, Tiberius, and Tswift reductions.

Our nominal retrieval setup for GJ 341b uses 11 free parameters. We include
the following six spectroscopically active gases that have (or partially have)
absorption bands in the NIRCam F444W wavelength range and are plausible
candidate gases for secondary atmospheres at temperatures at or near GJ 341b’s
equilibrium temperature of 540 K: H2O, CH4, CO, CO2, N2O, and O3. Each gas is
fitted in terms of the log10 evenly-mixed volume mixing ratio (VMR) with a
flat prior, VMR $\sim\mathcal{U}(-12,0)$. The bulk atmospheric constituent is
left agnostic by fitting for the mean molecular weight of the primary
constituent using a flat prior, $\mu\sim\mathcal{U}(2.5,50)$ g/mol, which
ranges from a low-mass mix of H2+He to heavier than CO2 (44 g/mol) and O3 (48
g/mol). Although aerosols are not explicitly included in our retrievals, we
fit for the opaque planet radius using a uniform prior
$R_{p}\sim\mathcal{U}(0.83,1.01)R_{\oplus}$, which could be due to a cloud
deck, the solid planetary surface (Benneke & Seager, 2012), additional
continuum opacity, or the critical refraction limit (Bétrémieux & Kaltenegger,
2013; Misra et al., 2014). We fit for the pressure at the opaque surface using
a uniform prior, $P_{0}\sim\mathcal{U}(-1,6)$ Pa. An isothermal temperature
structure is assumed with a uniform prior, $T\sim\mathcal{U}(300,700)$ K.
Finally, we allow the planet mass to vary uniformly between
$M_{p}\sim\mathcal{U}(0.3,2.2)M_{\oplus}$.

Our smarter retrievals use the dynesty (Speagle, 2020) nested sampling
algorithm (Skilling, 2004) with 600 live points. We run the model until the
contribution to the total evidence from the remaining prior volume drops below
dlogz=0.075.

Figure 9: Comparison of retrieval results on data reduced by the Eureka!,
Tiberius, and Tswift pipelines, respectively (rows of subplots). For each
reduction’s retrieval, the median spectral fit is shown in the left panel, the
1D marginalized CO2 posterior is shown in the middle panel, and the 1D
marginalized O3 posterior is shown in the right panel. Although either CO2 or
O3 are weakly detected in retrievals of the Tiberius and Tswift spectra,
respectively, (the Eureka! result is unconstrained) these findings are
incompatible with one another. Thus, relatively small differences in the
reductions can cause comparatively large changes in the subsequent physical
interpretation, including spurious molecular detections.

Figure 9 shows a summary of our retrieval results, comparing top-level
findings for each data reduction pipeline. Figures 16, 17, and 18 (in Appendix
E) show results for the Eureka!, Tiberius, and Tswift pipelines, respectively.
Based on the $\chi^{2}_{\nu}$ values from the median retrieved spectrum, the
retrievals are overfitting the spectrum, and therefore likely fitting noise.
Nonetheless, they reveal important aspects of our observational sensitivity to
the possible atmosphere of GJ 341b as well as the agreement/disagreement
between the reduction pipelines. As previously discussed, the Eureka!
reduction appears the most flat and as a result the retrieval rules out a
limited slice of parameter space that corresponds to only the thickest
(reference pressures $P_{0}>10^{4}$ Pa) and most extended atmospheres, with
high abundances of absorbing gases, that we explored.

On the other hand, the Tiberius and Tswift reductions tell a different story
due to their low significance detection of a spectroscopic feature (Table 3),
but they do not agree on the wavelength/molecule responsible for the feature,
implying that these are not real astrophysical signals. The retrieval results
for Tiberius favor high abundances of CO2 ($-2.1^{+1.2}_{-1.9}$ dex) with
moderate abundances of N2O ($-4.3^{+2.1}_{-3.9}$ dex) and low abundances of O3
($-8.3^{+2.5}_{-2.7}$ dex; 2$\sigma$ upper limit of $-3.6$ dex). The retrieval
results for Tswift show the opposite, favoring low abundances of CO2
($-8.7^{+2.6}_{-2.1}$ dex; 2$\sigma$ upper limit of $-3.6$ dex) and N2O
($-8.2^{+2.8}_{-2.4}$ dex; 2$\sigma$ upper limit of $-2.7$ dex) and high
abundances of O3 ($-1.8^{+1.0}_{-1.8}$ dex). These O3 and CO2 constraints
differ at ${\sim}1.5{\sigma}$. Our Bayesian model comparisons (using the
Bayesian evidence ratios computed by dynesty) find that the Tiberius reduction
implies a weak detection of CO2 (3.1$\sigma$), but no O3, while the Tswift
reduction instead leads to a weak detection of O3 (2.9$\sigma$), but no CO2.
Conversely, the Eureka! reduction favors no detected gases at ${>}1\sigma$
significance. However, both Tiberius and Tswift retrievals levy similar
constraints on the reference pressure, isothermal temperature, reference
radius, planet mass, and mean molecular weight which are all in agreement at
better than 1$\sigma$. These findings suggest that different noise instances
with otherwise similar noise properties dominate the character of each reduced
spectrum, as opposed to planetary characteristics. Furthermore, the
statistical significance of results from any one pipeline must be interpreted
as lower than the assessment of that spectrum in isolation would suggest.

##  5 Discussion and Conclusions

We have measured the 3.9–5.0 m transmission spectrum of the sub-Earth sized
exoplanet GJ 341b using JWST’s NIRCam instrument. Our transmission spectrum is
the result of co-adding the transmission spectra measured from three separate
transit observations and represents the first observation of this planet’s
atmosphere. We reduced the data using three independent pipelines (Eureka!,
Tiberius and Tswift) and performed several tests to determine the significance
of an atmosphere detection.

###  5.1 The NIRCam/F444W data quality

Our NIRCam/F444W data is dominated by an exponential ramp that varies between
the three visits (Figure 1). Despite fitting for this ramp, there is residual
red noise in the white light curve data (Figure 2). This red noise becomes
apparent at timescales between $\sim 1-30$ minutes. We attempted to remove
this red noise by using up to ten engineering telemetry timeseries
(‘mnemonics’) in our systematics model (Figure 3). Despite encouraging signs,
we were unable to significantly improve the residual red noise. We believe
that the JWST engineering mnemonics do present an opportunity to detrend JWST
spectrophotometry in general, however, we leave further attempts to do so to
future work. Fortunately, this red noise is not significant in the
spectroscopic light curves with our light curves reaching between 1.03 and
1.14$\times$ the noise limit (Fig. 7). As a result, we were able to obtain a
very precise transmission spectrum from our three coadded visits, with an
average uncertainty of 15 ppm in 0.02 µm bins. For these reasons, we recommend
the use of NIRCam/F444W for studies of bright exoplanet systems. While NIRCam
also provides us with short wavelength photometry (as used successfully in the
ERS program, Ahrer et al., 2023), our extraction of this here shows that is
not always beneficial.

###  5.2 GJ 341b’s transmission spectrum

The main conclusion we draw from our high precision transmission spectrum of
GJ 341b is that the planet has not retained a low mean molecular weight
atmosphere. This is consistent with expectations from the cosmic shoreline
(Zahnle & Catling, 2017) that small planets like GJ 341b should be stripped of
any primordial H/He-dominated atmosphere. Digging deeper, we find that the
very high precision of our data, combined with the small expected feature
amplitudes, leads our statistical tests, forward models and retrievals to
slightly different conclusions.

In our null hypothesis tests, the Eureka! spectrum was consistent with being
flat to within $0.2\sigma$ while the Tiberius and Tswift spectra were less
consistent at $1.3\sigma$ and $1.4\sigma$, respectively. However, we believe
that the relative consistency with being flat is driven by the relative size
of the uncertainties for each transmission spectrum. Repeating the null
hypothesis tests but switching the uncertainties between the reductions, we
find that Tiberius is inherently the flattest spectrum, followed by Eureka!
and then Tswift.

Our Gaussian feature tests, using the default uncertainties for each
reduction, find $2.2\sigma$ evidence against a feature for Eureka! and
$2.3\sigma$ and $2.9\sigma$ evidence for features in the Tiberius and Tswift
spectra, respectively. However, our atmospheric retrievals find that different
species are responsible for these candidate features: CO2 in the Tiberius
spectrum at $3.1\sigma$ and O3 in the Tswift spectrum at $2.9\sigma$, with a
relative difference between the CO2 and O3 abundance constraints of
${\sim}1.5\sigma$. These differences suggest that these ‘features’ are not
real astrophysical signals.

Our forward model analysis finds that each reduction’s spectrum rejects a
$<285\times$ solar metallicity atmosphere to $\geq 3\sigma$, with the Tswift
spectrum rejecting a $<1000\times$ solar metallicity atmosphere. Instead, the
forward models find our transmission spectra are consistent with no
atmosphere, a hazy atmosphere, or an atmosphere containing a species that does
not have prominent molecular bands across the NIRCam/F444W bandpass, e.g. H2O.
Furthermore, CH4-dominated atmospheres are ruled out to $1\sigma$–$3\sigma$
depending on the reduction, while CO2 atmospheres are consistent with the
Eureka! and Tiberius spectra to $<1\sigma$ and inconsistent with the Tswift
spectrum to $>2\sigma$.

At the radius and equilibrium temperature of GJ 341b, atmospheres dominated by
CO2, H2O or O2/O3, resulting from post-runaway atmospheres, are all plausible
and even predicted (Luger & Barnes, 2015; Gao et al., 2015; Schaefer et al.,
2016; Lincowski et al., 2018; Kite & Schaefer, 2021; Krissansen-Totton &
Fortney, 2022). Therefore, the single gas end member atmospheres that result
from our forward modelling, and that are consistent with our retrievals, are
in line with expectations. However, since a rocky exoplanet atmosphere has yet
to be detected and characterized, the possibilities are largely restricted to
theoretical predictions and the range of possible temperatures, pressures and
gas abundance combinations that are consistent with the observations is vast.
Furthermore the fact that each pipeline prefers a different gas demonstrates
that the evidence for these gases is not due to real astrophysical signals.

Our study of this shallow transit (240 ppm) demonstrates that small
differences between pipelines’ spectra can lead to differences in the
conclusions drawn from statistical tests, atmospheric forward models and
retrievals. Therefore, it is imperative that such studies use multiple
independent pipelines, two or more visits, and several statistical tests to
determine the robustness of future claims of atmospheric absorption for small
exoplanets.

###  5.3 The Cosmic Shoreline for M-dwarfs

GJ 341b is the fourth planet to be studied as part of our reconnaissance
program (PID:1981) to test the existence of an M-dwarf Cosmic Shoreline
(Zahnle & Catling, 2017). Our results so far are consistent with expectations
(Figure 10). GJ 341b (this work) and LHS 475b (Lustig-Yaeger & Fu et al.
2023b) sit on the dry/airless side of the shoreline, which is consistent with
our transmission spectra. GJ 486b (Moran & Stevenson et al. 2023) and GJ 1132b
(May & MacDonald et al. 2023) sit at the edge of the shoreline, where
atmospheres and atmosphere-less are both possible outcomes. GJ 486b shows
water in its transmission spectrum, although this could be stellar
contamination. GJ 1132b paints an inconsistent picture between the two visits’
transits, making inferences more challenging. Our remaining target
(TRAPPIST-1h) is likely to have an atmosphere based on expectations from the
shoreline. Our upcoming transit observations will shed light on the existence
of its atmosphere. Taken together, our program has demonstrated the tight
constraints that can be placed on Earth-sized exoplanet atmospheres.

Figure 10: The M-dwarf Cosmic Shoreline (Zahnle & Catling, 2017) plotted in
escape velocity–instellation space. Planets that sit below the shoreline (gray
line) are more likely to have atmospheres and those above are less likely. The
targets in our program are: GJ 341b (this paper), LHS 475b (Lustig-Yaeger & Fu
et al. 2023b), GJ 486b (Moran & Stevenson et al. 2023), GJ 1132b (May &
MacDonald et al. 2023) and TRAPPIST-1h (Rustamkulov & Peacock et al., in
prep.). We include the other TRAPPIST-1 planets for reference.

Looking beyond our program to other JWST studies of Earth-sized exoplanets
around M-dwarfs, TRAPPIST-1b and c are both consistent with having no
atmosphere (Greene et al., 2023; Lim et al., 2023; Zieba et al., 2023).
TRAPPIST-1b is close to the Cosmic Shoreline and thus the lack of an
atmosphere is not in tension with expectations. TRAPPIST-1c, on the other
hand, sits on the atmosphere side. However, the shoreline presented in Figure
10 is plotted in terms of the relative instellation. The XUV cosmic shoreline
would sit below this (Zahnle & Catling, 2017) and may be consistent with the
TRAPPIST-1c result. If the results of no atmospheres for these planets hold,
then this could be emerging evidence that XUV dominates the cosmic shoreline
in these late M-dwarf systems, consistent with past escape predictions
(Wheatley et al., 2017; Dong et al., 2018). Only by combining results from
wider and deeper surveys, including ongoing JWST programs (e.g., PID: 2512),
will we be able to confidently determine the boundaries of a cosmic shoreline
for M-dwarfs, if indeed one exists.

We thank the anonymous reviewer for their constructive comments that improved
the manuscript. We also thank Victoria DiTomasso for sharing the system
parameters of GJ 341b from her joint analysis of the radial velocity and TESS
data. JK acknowledges financial support from Imperial College London through
an Imperial College Research Fellowship grant. 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/kt9e-h986 (catalog DOI).

##  Appendix A The reduction and fitting approaches taken by Eureka!, Tiberius
and Tswift

Table 6: Comparison of the different reduction and analysis choices used by the three pipelines, along with the resulting noise properties. | Eureka! | Tiberius | Tswift
---|---|---|---
jwst pipeline version | 1.8.2 | 1.8.2 | 1.8.1
Group-level | col.-by-col. mean | None | None
background subtraction?a | & row-by-row mean |  |
Outlier detection | Full-frame: 4$\sigma$ | Full-frame: 5$\sigma$ | Full-frame: 5$\sigma$
/ sigma clipping | Light curves: 3.5$\sigma$ | Light curves: 4$\sigma$ | Light curves: 5$\sigma$
Pixel supersampling? | No | Yes, $10\times$ | No
Standard / optimal extraction | Optimal | Standard | Standard
Aperture full-width | 6 | 5 | 6
(pixels) |  |  |
Integration-level | None | col.-by-col. median | row-by-row median
background subtraction? |  |  | & col.-by-col. median
Spectral range extracted | 3.90–5.00 | 3.89–5.02 | 3.92–4.99
(m) |  |  |
Wavelength bin width | 0.02 | 0.02 | 0.02
(m) |  |  |
Number of integrations trimmed | First 300 from WLCb | 0 from WLC | None
| First 300 from SLCc | First 600 from SLC |
WLC model | Exponential $\times$ linear ramp | Exponential $\times$ linear ramp | Exponential $\times$ linear ramp
|  |  | (amplifier-dependent)
SLC model | Exponential $\times$ linear ramp | Linear ramp | Exponential $\times$ linear ramp
Common mode correction used? | No | Yes | Yes
|  | (visit dependent) | (visit & amplifier dependent)
WLC RMS | 223 | 292 | 303
(ppm) |  |  |
SLC RMS | 716 | 645 | 661
(ppm) |  |  |
Transit depth precision | 17 | 13 | 15
(ppm) |  |  |
aThe goal of group-level background subtraction is to remove 1/f noise.
bWLC corresponds to the white light curves.
cSLC corresponds to the spectroscopic light curves.

##  Appendix B The light curves and residuals from each pipeline

Figure 11: The light curves and best fit models resulting from the Tiberius
pipeline. The columns show, from left to right: transit 1, transit 2, transit
3. The rows show, from top to bottom: 1) the white light curve and residuals,
2) the spectroscopic light curves, 3) the best-fit models to the spectroscopic
light curves, 4) the residuals from the spectroscopic light curve fitting. The
vertical band at the left hand edge of the panels in rows 2–4 correspond to
the 600 integrations that were masked during the spectroscopic light curve
fitting stage. Figure 12: The light curves and best fit models resulting from
the Tswift pipeline. The columns show, from left to right: transit 1, transit
2, transit 3. The rows show, from top to bottom: 1) the white light curve and
residuals, 2) the spectroscopic light curves, 3) the best-fit models to the
spectroscopic light curves, 4) the residuals from the spectroscopic light
curve fitting.

##  Appendix C Transmission spectra measured for each visit

The transmission spectra measured for each visit by each reduction are shown
in Figure 13.

Figure 13: The transmission spectra for visit 1 (blue), visit 2 (orange) and
visit 3 (green) as measured by Tiberius (top panel), Eureka! (middle panel)
and Tswift (bottom panel).

##  Appendix D Forward models for all reductions

Figure 14: Atmospheric forward models compared to the Tiberius reduction from
JWST/NIRCam. Figure 15: Atmospheric forward models compared to the Tswift
reduction from JWST/NIRCam.

##  Appendix E Retrieval Posteriors for all Reductions

Figure 16: Atmospheric retrieval posteriors (lower left corner pairs) and
range of spectral fits (upper right) from our coadded NIRCam observations of
GJ 341b reduced with Eureka!. A subset of model parameters and covariances are
shown, with H2O, CH4, and CO not displayed (nor constrained) for clarity.
Analogous results are shown for the Tiberius and Tswift data reduction
pipelines in Figures 17 and 18, respectively. Figure 17: Atmospheric retrieval
posteriors (lower left corner pairs) and range of spectral fits (upper right)
from our coadded NIRCam observations of GJ 341b reduced with Tiberius. A
subset of model parameters and covariances are shown, with H2O, CH4, and CO
not displayed (nor constrained) for clarity.  Figure 18: Atmospheric retrieval
posteriors (lower left corner pairs) and range of spectral fits (upper right)
from our coadded NIRCam observations of GJ 341b reduced with Tswift. A subset
of model parameters and covariances are shown, with H2O, CH4, and CO not
displayed (nor constrained) for clarity.

## References

  * Ahrer et al. (2023) Ahrer, E.-M., Stevenson, K. B., Mansfield, M., et al. 2023, Nature, 614, 653, doi: 10.1038/s41586-022-05590-4
  * Alderson et al. (2023) Alderson, L., Wakeford, H. R., Alam, M. K., et al. 2023, Nature, 614, 664, doi: 10.1038/s41586-022-05591-3
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765, doi: 10.1098/rsta.2011.0269
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Batalha et al. (2020) Batalha, N., Freedman, R., Gharib-Nezhad, E., & Lupu, R. 2020, Resampled Opacity Database for PICASO, 2.0, Zenodo, doi: 10.5281/zenodo.6928501
  * Batalha et al. (2019) Batalha, N. E., Lewis, T., Fortney, J. J., et al. 2019, ApJ, 885, L25, doi: 10.3847/2041-8213/ab4909
  * Batalha et al. (2017) Batalha, N. E., Mandell, A., Pontoppidan, K., et al. 2017, ArXiv e-prints.  https://arxiv.org/abs/1702.01820
  * Bell et al. (2022) Bell, T. J., Ahrer, E.-M., Brande, J., et al. 2022, arXiv e-prints, arXiv:2207.03585.  https://arxiv.org/abs/2207.03585
  * Benneke & Seager (2012) Benneke, B., & Seager, S. 2012, ApJ, 753, 100, doi: 10.1088/0004-637X/753/2/100
  * Bétrémieux & Kaltenegger (2013) Bétrémieux, Y., & Kaltenegger, L. 2013, ApJ, 772, L31, doi: 10.1088/2041-8205/772/2/L31
  * Bourque et al. (2021) Bourque, M., Espinoza, N., Filippazzo, J., et al. 2021, The Exoplanet Characterization Toolkit (ExoCTK), 1.0.0, Zenodo, doi: 10.5281/zenodo.4556063
  * Bushouse et al. (2022) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022, JWST Calibration Pipeline, 1.8.2, Zenodo, doi: 10.5281/zenodo.7325378
  * Charbonneau & Deming (2007) Charbonneau, D., & Deming, D. 2007, arXiv e-prints, arXiv:0706.1047, doi: 10.48550/arXiv.0706.1047
  * Chen & Kipping (2017) Chen, J., & Kipping, D. 2017, ApJ, 834, 17, doi: 10.3847/1538-4357/834/1/17
  * de Wit et al. (2016) de Wit, J., Wakeford, H. R., Gillon, M., et al. 2016, Nature, 537, 69, doi: 10.1038/nature18641
  * de Wit et al. (2018) de Wit, J., Wakeford, H. R., Lewis, N. K., et al. 2018, Nature Astronomy, 2, 214, doi: 10.1038/s41550-017-0374-z
  * Diamond-Lowe et al. (2018) Diamond-Lowe, H., Berta-Thompson, Z., Charbonneau, D., & Kempton, E. M. R. 2018, AJ, 156, 42, doi: 10.3847/1538-3881/aac6dd
  * Diamond-Lowe et al. (2022) Diamond-Lowe, H., Mendonça, J. M., Charbonneau, D., & Buchhave, L. A. 2022, arXiv e-prints, arXiv:2210.11809, doi: 10.48550/arXiv.2210.11809
  * DiTomasso et al. (2023) DiTomasso, V., López-Morales, M., Peacock, S., & Malavolta, L. 2023, ApJL, submitted
  * Dong et al. (2018) Dong, C., Jin, M., Lingam, M., et al. 2018, Proceedings of the National Academy of Science, 115, 260, doi: 10.1073/pnas.1708010115
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Gao et al. (2015) Gao, P., Hu, R., Robinson, T. D., Li, C., & Yung, Y. L. 2015, ApJ, 806, 249, doi: 10.1088/0004-637X/806/2/249
  * Gardner et al. (2023) Gardner, J. P., Mather, J. C., Abbott, R., et al. 2023, PASP, 135, 068001, doi: 10.1088/1538-3873/acd1b5
  * 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
  * Greene et al. (2023) Greene, T. P., Bell, T. J., Ducrot, E., et al. 2023, Nature, 618, 39, doi: 10.1038/s41586-023-05951-7
  * Guillot (2010) Guillot, T. 2010, A&A, 520, A27, doi: 10.1051/0004-6361/200913396
  * 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
  * Horne (1986) Horne, K. 1986, Publ. Astron. Soc. Pac., 98, 609, doi: 10.1086/131801
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * Kirk et al. (2019) Kirk, J., López-Morales, M., Wheatley, P. J., et al. 2019, AJ, 158, 144, doi: 10.3847/1538-3881/ab397d
  * Kirk et al. (2021) Kirk, J., Rackham, B. V., MacDonald, R. J., et al. 2021, AJ, 162, 34, doi: 10.3847/1538-3881/abfcd2
  * Kite & Schaefer (2021) Kite, E. S., & Schaefer, L. 2021, ApJ, 909, L22, doi: 10.3847/2041-8213/abe7dc
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 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. (2019) Kreidberg, L., Koll, D. D. B., Morley, C., et al. 2019, Nature, 573, 87, doi: 10.1038/s41586-019-1497-4
  * Krissansen-Totton & Fortney (2022) Krissansen-Totton, J., & Fortney, J. J. 2022, ApJ, 933, 115, doi: 10.3847/1538-4357/ac69cb
  * Lammer et al. (2003) Lammer, H., Selsis, F., Ribas, I., et al. 2003, ApJ, 598, L121, doi: 10.1086/380815
  * Lim et al. (2023) Lim, O., Benneke, B., Doyon, R., et al. 2023, arXiv e-prints, arXiv:2309.07047, doi: 10.48550/arXiv.2309.07047
  * Lincowski et al. (2018) Lincowski, A. P., Meadows, V. S., Crisp, D., et al. 2018, ApJ, 867, 76, doi: 10.3847/1538-4357/aae36a
  * Line et al. (2014) Line, M. R., Knutson, H., Wolf, A. S., & Yung, Y. L. 2014, ApJ, 783, 70, doi: 10.1088/0004-637X/783/2/70
  * Line & Yung (2013) Line, M. R., & Yung, Y. L. 2013, ApJ, 779, 3, doi: 10.1088/0004-637X/779/1/3
  * Line et al. (2013) Line, M. R., Wolf, A. S., Zhang, X., et al. 2013, ApJ, 775, 137, doi: 10.1088/0004-637X/775/2/137
  * Luger & Barnes (2015) Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119, doi: 10.1089/ast.2014.1231
  * Luque & Pallé (2022) Luque, R., & Pallé, E. 2022, Science, 377, 1211, doi: 10.1126/science.abl7164
  * Lustig-Yaeger et al. (2023a) Lustig-Yaeger, J., Meadows, V. S., Crisp, D., Line, M. R., & Robinson, T. D. 2023a, arXiv e-prints, arXiv:2308.14804, doi: 10.48550/arXiv.2308.14804
  * Lustig-Yaeger et al. (2022) Lustig-Yaeger, J., Sotzen, K. S., Stevenson, K. B., et al. 2022, AJ, 163, 140, doi: 10.3847/1538-3881/ac5034
  * Lustig-Yaeger et al. (2023b) Lustig-Yaeger, J., Fu, G., May, E. M., et al. 2023b, arXiv e-prints, arXiv:2301.04191, doi: 10.48550/arXiv.2301.04191
  * Magic et al. (2015) Magic, Z., Chiavassa, A., Collet, R., & Asplund, M. 2015, A&A, 573, A90, doi: 10.1051/0004-6361/201423804
  * May et al. (2023) May, E. M., MacDonald, R. J., Bennett, K. A., et al. 2023, ApJL, submitted
  * McElwain et al. (2023) McElwain, M. W., Feinberg, L. D., Perrin, M. D., et al. 2023, PASP, 135, 058001, doi: 10.1088/1538-3873/acada0
  * Ment et al. (2023) Ment, K., Charbonneau, D., Irwin, J., et al. 2023, arXiv e-prints, arXiv:2304.01920, doi: 10.48550/arXiv.2304.01920
  * Misra et al. (2014) Misra, A., Meadows, V., & Crisp, D. 2014, ApJ, 792, 61, doi: 10.1088/0004-637X/792/1/61
  * Moran et al. (2023) Moran, S. E., Stevenson, K. B., Sing, D. K., et al. 2023, ApJ, 948, L11, doi: 10.3847/2041-8213/accb9c
  * Owen & Jackson (2012) Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931, doi: 10.1111/j.1365-2966.2012.21481.x
  * Owen & Mohanty (2016) Owen, J. E., & Mohanty, S. 2016, MNRAS, 459, 4088, doi: 10.1093/mnras/stw959
  * Peacock et al. (2019) Peacock, S., Barman, T., Shkolnik, E. L., Hauschildt, P. H., & Baron, E. 2019, ApJ, 871, 235, doi: 10.3847/1538-4357/aaf891
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Pizzolato et al. (2003) Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 397, 147, doi: 10.1051/0004-6361:20021560
  * Preibisch & Feigelson (2005) Preibisch, T., & Feigelson, E. D. 2005, ApJS, 160, 390, doi: 10.1086/432094
  * Rieke et al. (2023) Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001, doi: 10.1088/1538-3873/acac53
  * Rigby et al. (2023) Rigby, J., Perrin, M., McElwain, M., et al. 2023, PASP, 135, 048001, doi: 10.1088/1538-3873/acb293
  * Rustamkulov et al. (2022) Rustamkulov, Z., Sing, D. K., Liu, R., & Wang, A. 2022, ApJ, 928, L7, doi: 10.3847/2041-8213/ac5b6f
  * Rustamkulov et al. (2023) Rustamkulov, Z., Sing, D. K., Mukherjee, S., et al. 2023, Nature, 614, 659, doi: 10.1038/s41586-022-05677-y
  * Salvatier et al. (2016) Salvatier, J., Wieckiâ, T. V., & Fonnesbeck, C. 2016, PyMC3: Python probabilistic programming framework, Astrophysics Source Code Library, record ascl:1610.016.  http://ascl.net/1610.016
  * Schaefer et al. (2016) Schaefer, L., Wordsworth, R. D., Berta-Thompson, Z., & Sasselov, D. 2016, ApJ, 829, 63, doi: 10.3847/0004-637X/829/2/63
  * Schlawin et al. (2020) Schlawin, E., Leisenring, J., Misselt, K., et al. 2020, AJ, 160, 231, doi: 10.3847/1538-3881/abb811
  * Shkolnik & Barman (2014) Shkolnik, E. L., & Barman, T. S. 2014, AJ, 148, 64, doi: 10.1088/0004-6256/148/4/64
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405, doi: 10.1063/1.1835238
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * STScI Development Team (2013) STScI Development Team. 2013, pysynphot: Synthetic photometry software package, Astrophysics Source Code Library, record ascl:1303.023.  http://ascl.net/1303.023
  * Swain et al. (2021) Swain, M. R., Estrela, R., Roudier, G. M., et al. 2021, AJ, 161, 213, doi: 10.3847/1538-3881/abe879
  * The JWST Transiting Exoplanet Community Early Release Science Team et al. (2022) The JWST Transiting Exoplanet Community Early Release Science Team, Ahrer, E.-M., Alderson, L., et al. 2022, arXiv e-prints, arXiv:2208.11692.  https://arxiv.org/abs/2208.11692
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Wakeford et al. (2019) Wakeford, H. R., Lewis, N. K., Fowler, J., et al. 2019, AJ, 157, 11, doi: 10.3847/1538-3881/aaf04d
  * Wheatley et al. (2017) Wheatley, P. J., Louden, T., Bourrier, V., Ehrenreich, D., & Gillon, M. 2017, MNRAS, 465, L74, doi: 10.1093/mnrasl/slw192
  * Zahnle & Catling (2017) Zahnle, K. J., & Catling, D. C. 2017, The Astrophysical Journal, 843, 122, doi: 10.3847/1538-4357/aa7846
  * Zieba et al. (2023) Zieba, S., Kreidberg, L., Ducrot, E., et al. 2023, Nature, 620, 746, doi: 10.1038/s41586-023-06232-z