Draft version April 2, 2024

Typeset using LATEX twocolumn style in AASTeX631

JWST COMPASS: NIRSpec/G395H Transmission Observations of the Super-Earth TOI-836b

Lili Alderson ,1 Natasha E. Batalha ,2 Hannah R. Wakeford ,1 Nicole L. Wallack ,3

Artyom Aguichine ,4 Johanna Teske ,3 Jea Adams Redai ,5 Munazza K. Alam ,6 Natalie M. Batalha ,4
Peter Gao ,3 James Kirk ,7 Mercedes López-Morales ,5 Sarah E. Moran ,8 Nicholas Scarsdale ,4
Nicholas F. Wogan ,2 and Angie Wolfgang 9
1 School of Physics, University of Bristol, HH Wills Physics Laboratory, Tyndall Avenue, Bristol BS8 1TL, UK
2 NASA Ames Research Center, Moffett Field, CA 94035, USA
arXiv:2404.00093v1 [astro-ph.EP] 29 Mar 2024

3 Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road, NW, Washington, DC 20015, USA
4 Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA
5 Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA
6 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
7 Department ofPhysics, Imperial College London, London, UK
8 Department of Planetary Sciences and Lunar and Planetary Laboratory, University of Arizona, Tuscon, AZ, USA
9 Eureka Scientific Inc., 2452 Delmer Street Suite 100, Oakland, CA 94602-3017

ABSTRACT
We present two transit observations of the ∼870 K, 1.7 R⊕ super-Earth TOI-836b with JWST NIR-
Spec/G395H, resulting in a 2.8–5.2 µm transmission spectrum. Using two different reduction pipelines,
we obtain a median transit depth precision of 34 ppm for Visit 1 and 36 ppm for Visit 2, leading to
a combined precision of 25 ppm in spectroscopic channels 30 pixels wide (∼ 0.02 µm). We find that
the transmission spectrum from both visits is well fit by a zero-sloped line by fitting zero-sloped and
sloped lines, as well as step functions to our data. Combining both visits, we are able to rule out atmo-
spheres with metallicities < 250×Solar for an opaque pressure level of 0.1 bar, corresponding to mean
molecular weights of ≲ 6 g mol−1 . We therefore conclude that TOI-836b does not have a H2 -dominated
atmosphere, in possible contrast with its larger, exterior sibling planet, TOI-836c. We recommend that
future proposals to observe small planets exercise caution when requiring specific numbers of transits
to rule out physical scenarios, particularly for high metallicities and planets around bright host stars,
as PandExo predictions appear to be more optimistic than that suggested by the gains from additional
transits implied by our data.

Keywords: Exoplanet atmospheric composition (2021); Exoplanet atmospheres (487); Exoplanets

(498); Infrared spectroscopy (2285)

1. INTRODUCTION planets were not massive enough to retain their primor-

Despite their ubiquity (Batalha et al. 2013), 1–4 R⊕ dial atmospheres (e.g., Lopez et al. 2012). Both photo-
exoplanets present some of the most complex challenges evaporation (e.g., Owen & Jackson 2012; Owen & Wu
for observation and interpretation, hampered further by 2013) and core-powered mass loss (e.g., Ginzburg et al.
the lack of solar system counterparts. In this radius 2018) can generate heat-driven hydrodynamic outflows
regime, exoplanets are typically split between two cate- from the upper atmosphere, but it is not clear which
gories: the larger sub-Neptunes (> 1.8 R⊕ ), which may mechanism is dominant (e.g., Rogers et al. 2021; Owen
have hydrogen-rich envelopes (Lopez & Fortney 2013; & Schlichting 2023).
Buchhave et al. 2014), and the smaller super-Earths, Regardless of which category an individual exoplanet
with likely more tenuous (if any) atmospheres (Rogers may lie in, the densities of planets in this parameter
2015; Rogers et al. 2021), while a dearth of planets ex- space typically allow for many possible interior compo-
ists in between (known as the radius valley, e.g., Fulton sitions (e.g., Rogers & Seager 2010; Zeng et al. 2019).
et al. 2017). Currently, one of the leading theories pro- Without the ability to directly probe the interiors of
posed to explain the radius valley is that the smaller these exoplanets, observing their atmospheres remains
2

the only way to understand this population in detail. spectrum. In this case, May & MacDonald et al. 2023
Indeed, studies suggest that determining atmospheric find that these discrepancies are most likely due to an
mean molecular weights or metallicities can help break unlucky random noise draw.
some of the degeneracies presented by interior structure In order to draw conclusions about the broader small
modelling and constrain bulk compositions (Figueira exoplanet population as a whole, the JWST COMPASS
et al. 2009; Rogers et al. 2011; Fortney et al. 2013). (Compositions of Mini-Planet Atmospheres for Statisti-
The advent of JWST and its exquisite precision across cal Study) Program (GO-2512, PIs N. E. Batalha & J.
a wide wavelength range (Ahrer et al. 2023b; Alderson Teske) is focusing on observing a statistically motivated
et al. 2023; Feinstein et al. 2023; Rustamkulov et al. sample of 1–3 R⊕ planets. The program will obtain NIR-
2023) unlocks the ability to explore small exoplanets in Spec/G395H transmission spectra of eleven exoplanets,
detail. Broad coverage of the infra-red (IR) enables the while the full statistical sample includes 12 planets, with
detection of a variety of molecular species that were pre- four pairs in the same systems.1 The targets were se-
viously inaccessible (Batalha et al. 2017), and critically lected from a subset of the ≤ 3 R⊕ planets observed as
offers the opportunity to explore wavelengths that are part of the Magellan-TESS Survey (MTS, Teske et al.
less prone to muting from clouds and hazes that have 2021), in order to understand to what extent small plan-
plagued the observation of small exoplanet atmospheres ets have detectable atmospheres, and explore the com-
with the Hubble Space Telescope and ground-based in- positional diversity of the population as a whole. Sim-
struments (e.g., Crossfield et al. 2013; Kreidberg et al. ilarly to the MTS targets, the COMPASS targets were
2014; Louden et al. 2017; Kirk et al. 2018; Ahrer et al. chosen using a quantitatively selected sample using a
2023a). Of particular interest to the study of super- merit function based on RP , insolation flux, stellar ef-
Earth atmospheres is JWST’s NIRSpec/G395H mode fective temperature, and exposure time required to ob-
(Jakobsen et al. 2022; Birkmann et al. 2022). Span- tain 30 ppm precision in an R∼100 NIRSpec/G395H bin
ning 2.87–5.14 µm at R∼2700, the G395H grating cov- at 4 µm (see Batalha et al. 2023). Specifically, Batalha
ers spectral features from the major absorption bands et al. (2023) showed that a quantitatively chosen sam-
of CO2 , CO and CH4 as well as a partial band of ple was shown to be a useful method for enabling con-
H2 O – molecules expected to sculpt super-Earth trans- straints on population-level parameters, which is the ul-
mission spectra across a variety of metallicities (e.g., timate goal of the COMPASS Program. By observing
Wordsworth & Kreidberg 2022). Furthermore, since multi-planet systems, the COMPASS Program also has
super-Earths have predominately been found around the ability to test a variety of formation and evolution
bright stars, G395H’s brightness limit affords the ability theories that are heavily dependent on insolation flux.
to observe these atmospheres without saturating. As- The first multi-planet system observed in the COM-
sessing the presence or absence of these four molecules PASS Program is that of TOI-836 (HIP 73427), in which
provides a zeroth-order assessment of the carbon-to- two planets are orbiting a bright (J mag ∼ 7.58) K-
oxygen ratio (C/O) of the atmosphere (Batalha et al. dwarf (Hawthorn et al. 2023). The larger of the planets,
2023). the sub-Neptune TOI-836.01 (planet c), has a bulk den-
Historically, Hubble and ground-based observations sity consistent with a gaseous envelope, with a radius of
of super-Earth (R<1.8 R⊕ ) atmospheres have typically 2.59±0.09 R⊕ and a mass of 9.6±2.7 M⊕ , and orbits on
yielded featureless transmission spectra (e.g., Diamond- a period of 8.59 days. The smaller TOI-836.02 (planet
Lowe et al. 2020; Libby-Roberts et al. 2022; Diamond- b) is interior, on a period of 3.81 days and Teq ∼870 K.
Lowe et al. 2023), but even with the power of JWST, With a radius of 1.70±0.06 R⊕ and mass of 4.5±0.9 M⊕ ,
exploring these exoplanets has not been without chal- TOI-836b is a super-Earth at the lower edge of the ra-
lenges. NIRSpec/G395H observations of GJ 486b show dius valley, and likely has a much smaller gas fraction.
evidence of a deviation from a flat line consistent with Given their positions respective to the radius valley,
either a water-rich atmosphere or with contamination the TOI-836 system presents an excellent opportunity
from unocculted starspots (Moran & Stevenson et al. to directly compare and contrast the atmospheres of
2023). The transmission spectra also showed a consis- differently-sized exoplanets that formed within the same
tent offset between the two NIRSpec detectors, poten- stellar environment. Here, we focus on TOI-836b, pre-
tially due to the superbias subtraction step in the data senting the 3–5 µm transmission spectrum before taking
reduction. Similar conflicting inferences have been seen
for GJ 1132b, where the transmission spectrum obtained 1 Our full sample includes L 98-59c and L 98-59d, the latter of
during one visit is consistent with a water-dominated at- which is being observed by GTO-1224, PI S. Birkmann.
mosphere, while the second visit presents a featureless
 3

Table 1. System Properties for TOI-836b. 3. DATA REDUCTION

Values used in the light curve fitting are shown
in Table 2. To check for consistency and ensure robust con-
clusions, we reduced the data using two independent
Property Value pipelines: ExoTiC-JEDI (Alderson et al. 2022, 2023) and
Eureka! (Bell et al. 2022). Each reduction process is
K (mag) 6.804 ± 0.018
described in detail below and follows similar procedures
R∗ (R⊙ ) 0.665±0.010
to other NIRSpec/G395H transmission spectra analyses.
T∗ (K) 4552±154
log(g) 4.743±0.105
3.1. ExoTiC-JEDI
[Fe/H]∗ -0.284±-0.067
Period (days) 3.81673 ± 0.00001
The Exoplanet Timeseries Characterisation - JWST
Extraction and Diagnostic Investigator (ExoTiC-JEDI)
MP (M⊕ ) 4.5+0.92
−0.86
package2 performs end-to-end extraction, reduction, and
RP (R⊕ ) 1.704 ± 0.067
analysis of JWST time-series data from uncal files
Teq (K) 871±36
through to light curve fitting to produce planetary spec-
e 0.053±0.042
tra. Throughout, NRS1 and NRS2 data are reduced
ω (°) 9 ± 92
independently, and each visit is treated separately. In
Semi-major axis (AU) 0.04220 ± 0.00093
all cases, we tried a variety of values for each reduction
All values from Hawthorn et al. (2023) parameter, and determined the value which resulted in
the smallest out-of-transit scatter in the resulting white
light curve.
We begin with a modified version of Stage 1 of the
a broader look at the system as a whole, including the jwst pipeline (v.1.8.6, context map 1078; Bushouse
observations of TOI-836c presented in Wallack & COM- et al. 2022), performing linearity, dark current and sat-
PASS et al. (2024). uration corrections, and using a jump detection thresh-
In §2, we describe our observations and detail our re- old of 15. We next perform a custom destriping routine
duction procedures in §3. We present the transmission to remove group level 1/f noise, masking the spectral
spectrum of TOI-836b in §4, and interpret the trans- trace 15σ from the dispersion axis for each integration,
mission spectrum using 1D radiative-convective atmo- subtracting the median pixel value of non-masked pixels
spheric models in §5. Finally, we discuss the implica- from each detector column in each group. We also per-
tions of our results on the TOI-836 system and for fu- form a custom bias subtraction, building a pseudo-bias
ture observations in §6 and summarise our conclusions image by computing the median of each detector pixel in
in §7. the first group across every integration in the time series.
This median image is then used in place of a bias and
subtracted from every group, and was found to improve
2. OBSERVATIONS
the out-of-transit scatter for both detector white light
We observed two transits of TOI-836b with JWST curves for both visits (see also Alderson et al. 2023).
NIRSpec using the high-resolution (R∼2700) G395H We then proceed with the standard ramp fitting step.
mode, which commenced on March 4 2023 at 18:09 UTC ExoTiC-JEDI also utilises Stage 2 of the jwst pipeline
and March 8 2023 at 13:45 UTC, respectively. These to produce the 2D wavelength map needed to obtain the
visits were coincidentally separated by one orbital pe- wavelength solution.
riod, considerably less than the 22-day rotation period In Stage 3 of ExoTiC-JEDI, we extract our 1D stellar
of TOI-836 (Hawthorn et al. 2023). NIRSpec/G395H spectra, performing additional cleaning steps and 1/f
provides spectroscopy between 2.87–5.14 µm across the correction. Using the standard data quality flags pro-
NRS1 and NRS2 detectors (with a gap between 3.72– duced by the jwst pipeline, replacing any pixels flagged
3.82 µm). Both observations were taken in NIRSpec as do not use, saturated, dead, hot, low quantum effi-
Bright Object Time Series (BOTS) mode using the ciency or no gain value with the median of the neigh-
SUB2048 subarray, F290LP filter, S1600A1 slit, and bouring 4 pixels in each row. To replace any spurious
NRSRAPID readout pattern. Each 5.3-hour observa- pixels that have not yet been corrected (such as cos-
tion consisted of 5259 integrations with 3 groups per in- mic rays), we identify any that are outliers from their
tegration, and was designed to cover the 1.8-hour transit
and sufficient pre- and post-transit baseline.
2 https://github.com/Exo-TiC/ExoTiC-JEDI
4

nearest neighbours on the detector, or throughout the light curves, we fit for Rp /R∗ , holding T0 , i and a/R∗
time series. We use a 20σ threshold in time and a 6σ fixed to the respective white light curve fit value, as
threshold spatially, replacing the problem pixel with the shown in Table 2.
median of that pixel in the surrounding 10 integrations For both the white and spectroscopic light curves, we
or 20 pixels in the row, respectively. Any remaining 1/f removed any data points that were greater than 4σ out-
noise and background are removed by subtracting the liers in the residuals, and refit the light curves until no
median unilluminated pixel value from each column in such points remained. We also rescaled the flux time
each integration. To extract the 1D stellar spectra, we series errors using the beta value (Pont et al. 2006) as
fit a Gaussian to each column to obtain the centre and measured from the white and red noise values calcu-
width of the spectral trace across the detector, fitting lated using the extra functions.noise calculator()
a fourth-order polynomial to each. The spectral trace in ExoTiC-JEDI to account for any remaining red noise
centres and widths are then smoothed with a median fil- in the data. We removed the first 370 integrations (∼
ter of window size 11 to determine the aperture region. 22 minutes) of visit 1 and the first 440 integrations (∼
For both visits and both detectors, we used an aperture 26 minutes) of visit 2 to remove settling ramps at the
five times the FWHM of the trace, approximately 8 pix- start of the observations. We additionally removed 259
els wide from edge to edge. An intrapixel extraction is integrations from the end of visit 1 (∼ 15 minutes), and
used to obtain the 1D stellar spectra, where intrapixel is 508 integrations from the end of visit 2 (∼ 30 minutes),
defined as the fraction of the FWHM which falls on each which removed a slight linear slope in the residuals of
pixel, such that at the edge of the aperture, the flux in- NRS1 fits and removed an ∼10 ppm offset between the
cluded from any intersected pixel is equal to the fraction transit depths of NRS1 and NRS2. The ExoTiC-JEDI
of the pixel within the aperture, multiplied by the total fitted white light curves and residuals for each visit are
flux value of that pixel. The 1D stellar spectra are then shown in Figure 1.
cross-correlated to obtain the x- and y-positional shifts
throughout the observation for use in systematic light 3.2. Eureka!
curve detrending.
For our second independent reduction, we utilise
Finally, we fit white light curves for both NRS1 and
Eureka!, an end-to-end pipeline for analysing JWST
NRS2, as well as spectroscopic light curves across the full
transiting planet data. We used the default proce-
NIRSpec/G395H wavelength range at a variety of reso-
dures for Stage 1 and Stage 2 from the Eureka! wrap-
lutions for both visits. For the white light curves (span-
per of the jwst pipeline, following the same proce-
ning 2.814–3.717 µm for NRS1 and 3.824–5.111 µm for
dures as ExoTiC-JEDI but using the standard super-
NRS2), we fit for the system inclination, i, ratio of semi-
bias subtraction. We also used the aforementioned
major axis to stellar radius, a/R∗ , centre of transit time,
ExoTiC-JEDI group-level background subtraction to ac-
T0 , and the ratio of planet to stellar radii, Rp /R∗ , hold-
count for 1/f noise. Following this, we use Eureka!
ing the period and eccentricity, e, and argument of pe-
Stage 3 to extract the stellar spectra and produce the
riastron, ω, fixed to values presented in Hawthorn et al.
broadband and spectroscopic light curves. Note that
(2023). The stellar limb darkening coefficients are held
while ExoTiC-JEDI maintains the slightly curved shape
fixed to values calculated using the ExoTiC-LD package
of the NIRSpec/G395H spectral trace throughout spec-
(Grant & Wakeford 2022) based on the stellar T∗ , log(g),
tral extraction, the Eureka! pipeline flattens this trace
and [Fe/H]∗ presented in Hawthorn et al. (2023) (see Ta-
by bringing the centre of mass of each column to the
ble 1), with Set One of the MPS-ATLAS stellar models
centre of the subarray, allowing for a straight box ex-
(Kostogryz et al. 2022, 2023) and the non-linear limb
traction to be used. Eureka! allows for the customisa-
darkening law (Claret 2000). We used a least-squares
tion of a variety of reduction parameters in Stage 3,
optimiser to fit for a batman (Kreidberg 2015) transit
including the trace extraction width, the region and
model simultaneously with our systematic model S(λ),
method for the background extraction, and trace ex-
which took the form
traction parameters. To find the best combination of
S(λ) = s0 + (s1 × t) + (s2 × xs |ys |), values, we tested extraction apertures consisting of com-
binations of 4-8 pixels from the centre of the flattened
where xs is the x-positional shift of the spectral trace, trace, background apertures of 8-11 pixels, sigma thresh-
|ys | is the absolute magnitude of the y-positional shift of olds for optimal extraction outlier rejection of 10 and
the spectral trace, t is the time and s0 , s1 , s2 are coeffi- 60 (which approximates standard extraction), and two
cient terms, as previously used for ExoTiC-JEDI analysis different methods of background subtraction (an addi-
in May & MacDonald et al. 2023. For the spectroscopic tional column-by-column mean subtraction and a full
 5

ExoTiC-JEDI White Light Curves ExoTiC-JEDI Residuals

1.003 Visit 1 0.002
1.002
Relative Flux

NRS1 0.001
1.001
1.000 0.000
0.999 NRS2
0.001
2 1 0 1 2 2 1 0 1 2
1.003
Visit 2 0.002
1.002
Relative Flux

1.001 NRS1 0.001

1.000 0.000
0.999 NRS2
0.001
2 1 0 1 2 2 1 0 1 2

Eureka! White Light Curves Eureka! Residuals

1.003 Visit 1 0.002
Relative Flux

1.002 NRS1 0.001

1.001
1.000 0.000
0.999 NRS2 0.001
2 1 0 1 2 2 1 0 1 2

1.003 Visit 2 0.002

1.002
Relative Flux

NRS1 0.001
1.001
1.000 0.000
0.999 NRS2
0.001
2 1 0 1 2 2 1 0 1 2

Time From Mid-Transit (Hours)

Figure 1. White light curves for each reduction, the best-fit models to those white light curves, and the associated residuals
for ExoTiC-JEDI and Eureka! for each visit and detector. Binned light curves and residuals are also shown in the alternate
colours. Histograms of the residuals are shown in the rightmost column.
6

Table 2. Best fit values for the four individual white light curve fits for ExoTiC-JEDI and Eureka! as shown in Figure 1 and
the Eureka! joint fit.
T0 (MJD) a/R∗ i (◦ ) Rp /R∗
Hawthorn et al. (2023) Value 581599.9953±2e-3 – 87.57±0.44 0.0235 ± 0.0013
NRS1 60007.86562 ± 7e−5 15.42 ± 1.01 88.10 ± 0.49 0.02458 ± 0.00016
Visit 1
NRS2 60007.86552 ± 8e−5 15.65 ± 1.20 88.19 ± 0.58 0.02489 ± 0.00018
ExoTiC-JEDI
NRS1 60011.68218 ± 7e−5 14.12 ± 1.17 87.45 ± 0.59 0.02478 ± 0.00015
Visit 2
NRS2 60011.68213 ± 9e−5 15.01 ± 1.38 87.91 ± 0.78 0.02448 ± 0.00017
NRS1 60007.86570 ± 6e−5 17.50 ± 0.77 89.70 ± 0.82 0.02451 ± 0.00013
Visit 1
NRS2 60007.86560 ± 8e−5 15.75 ± 1.25 88.24 ± 0.74 0.02489 ± 0.00018
NRS1 60011.68205 ± 7e−5 14.12 ± 1.18 87.45 ± 0.60 0.02503 ± 0.00020
Visit 2
NRS2 60011.68207 ± 9e−5 15.08 ± 1.49 87.94 ± 1.06 0.02450 ± 0.00019
Eureka!
60007.86552 ± 5e−5
NRS1 14.99 ± 0.84 87.89 ± 0.42 0.02482 ± 0.00013
60011.68225 ± 5e−5
Joint
60007.86548 ± 6e−5
NRS2 15.05 ± 0.89 87.92 ± 0.44 0.02469 ± 0.00014
60011.68221 ± 6e−5
 7

frame median subtraction). We select the final reduc- During our spectroscopic light curve fits, we once
tion parameters as the combination that minimises the again fit for i, a/R∗ , T0 , and Rp /R∗ , which results in pa-
scatter in the resulting white light curves, doing this sep- rameters that are consistent with the fitted white light
arately for each detector and each visit. We find that curve values to within 2σ. We utilise a prior when fit-
for both visits, both NRS1 and NRS2 favoured an addi- ting the i, a/R∗ , T0 for the spectroscopic light curves.
tional column-by-column background subtraction using To obtain priors for these fits, we utilise the posterior
a sigma threshold of 60 (which approximates a standard distribution for each free parameter that resulted from
box extraction). The optimal aperture half-widths for the MCMC chains of the white light curve fits for each
the trace extraction for NRS1 and NRS2 was 6 pixels visit. We use the median 3σ values from combining the
for visit 1 and 4 pixels for visit 2. The background sub- chains from NRS1 and NRS2 as Gaussian priors for each
traction region spans from the upper and lower edge of astrophysical parameter, meaning that the prior repre-
the detector subarray to an inner bound defined by a sents the combined constraints from NRS1 and NRS2.
number of pixels away from the flattened trace. This For Rp /R∗ , we use a flat uninformed prior to not bias
inner bound was found to be 8 pixels in both detectors our transmission spectrum.
for Visit 1 and 9 pixels for NRS1 and 8 pixels for NRS2
3.2.1. Joint Fit of the Eureka! Light Curves
for Visit 2. We then extract white light and 30-pixel
binned (∼ 0.02 µm, R∼200) spectroscopic light curves In order to evaluate the power of combining multiple
for both visits for both detectors. visits, we also produce a joint fit of the Eureka! reduc-
During the light curve fitting stage, we move away tion light curves. Here, we fit both visits simultaneously
from the Eureka! pipeline and utilise a custom light but continue to separate NRS1 and NRS2 to account for
curve fitting code, but refer to this reduction as the any offsets between the two detectors and the differing
“Eureka!” reduction for simplicity. We fit each white systematic effects. We follow the same procedures as
and spectroscopic light curve separately using emcee for the individual fits (§3.2), but now obtain a universal
(Foreman-Mackey et al. 2013), fitting for i, a/R∗ , T0 , value for i, a/R∗ , T0 , and Rp /R∗ for each detector in
and Rp /R∗ and fixing the other orbital parameters to both the white and spectroscopic light curves. In the
those from Table 1 using the batman package (Kreidberg case of T0 , we assumed a centre of transit time for each
2015). We utilise quadratic limb-darkening coefficients visit and fit for a common offset from this value which
calculated using ExoTiC-LD and the stellar parameters applies to both visits for each detector (note that the
from Table 1. We fit our transit model and a systematic visits are separated by a single orbital period, see Table
model simultaneously, which took the form 2). We again use an MCMC fit, initialising our walkers
using a Levenberg-Marquardt least-squares minimizer,
S(λ) = s0 + s1 × T + s2 × X + s3 × Y , with three times the number of walkers as free parame-
ters in our fit (resulting in 42 walkers). For the MCMC
where X and Y are the normalised positions of the trace we discard a burn-in of 50,000 steps before utilising a
on the detector and si are the free parameters in our production run of 50,000 steps. Our best-fit joint white
systematic noise model. We use an iterative rolling me- light curve parameters are shown in Table 2 in compar-
dian outlier rejection with a 50 data point-wide window ison to those from the individual fits. For the spectro-
three times on both the white and spectroscopic light scopic light curves, we once again fit a Gaussian to the
curves, removing outliers more than 3σ from the me- posterior distribution for each free parameter which re-
dian. We initialise our MCMC walkers using the best fit sulted from the MCMC chains, and apply this as our
results from a Levenberg–Marquardt least-squares min- prior for the respective parameters. As with the indi-
imisation. We utilise three times the number of free pa- vidual fits, for Rp /R∗ we use a flat uninformed prior
rameters as the number of walkers (resulting in 27 walk- to not bias our transmission spectrum. For this joint
ers) and run a burn-in of 50,000 steps which is discarded fit, the fitted spectroscopic light curve parameters agree
followed by a production run of 50,000 steps, with un- with the white light curve values to within 1σ in all
informed priors on all of the parameters. We trim the wavelength channels.
initial 444 points (∼20 minutes) from all the light curves
to remove any initial ramp that may be present, and re- 4. TRANSMISSION SPECTRUM
moved the last 259 points of Visit 1 and 508 points of The 3–5 µm transmission spectra of TOI-836b using
Visit 2 (see Section 3.1). The fitted white light curves a 30-pixel binning scheme for each of the two visits are
and residuals resulting from the Eureka! reduction for shown in the upper panels of Figure 2, where no offsets
each visit are shown in Figure 1, while the fitted white have been applied between NRS1 and NRS2 or between
light curve parameters are shown in Table 2. the visits. In general, each visit appears to be consistent
8

900 ExoTiC-JEDI Visit 1 Eureka! Visit 1

800 ExoTiC-JEDI Visit 2 Eureka! Visit 2
Transit Depth

700
(ppm)

600
500
400
200
Difference
(ppm)

0
200
3.0 3.5 4.0 4.5 5.0 3.0 3.5 4.0 4.5 5.0
Wavelength ( m) Wavelength ( m)

ExoTiC-JEDI Weighted Average

800 Eureka! Weighted Average
Eureka! Joint Fit
Transit Depth

700
(ppm)

600

500

Eureka! Weighted - Eureka! Joint ExoTiC JEDI - Eureka! Weighted ExoTiC JEDI - Eureka! Joint
50
Difference
(ppm)

0
50
3.0 3.5 4.0 4.5 5.0
Wavelength ( m)
Figure 2. Top, Upper Panels: Individual visit transmission spectra for ExoTiC-JEDI (left, purples) and Eureka! (right, greens).
Lower Panel: Difference between individual visit transmission spectra for ExoTiC-JEDI (left) and Eureka! (right). On average,
the ExoTiC-JEDI reductions for visit 1 and visit 2 are consistent to within 39 ppm, while the Eureka! reductions agree to
within 48 ppm. The ExoTiC-JEDI and Eureka! reductions are consistent with each other to within the median transit depth
uncertainty for both visit 1 and visit 2. Bottom, Upper Panel: Weighted average transmission spectrum from the two visits
from the ExoTiC-JEDI (purple) and Eureka! (light green) reductions and joint fit transmission spectrum from the Eureka!
reduction (dark green). Lower Panel: Difference between each of the combined ExoTiC-JEDI and Eureka! transmission spectra
in ppm. As the difference between the two Eureka! methods is less the 5 ppm (black line), the two ExoTiC-JEDI– Eureka! lines
are difficult to distinguish (coloured lines). On average, the combined visit Eureka! and ExoTiC-JEDI spectra are consistent to
within 10 ppm.

between each reduction method, with a median differ- Visit 1 and 17 ppm for Visit 2, both less than their re-
ence in transit depth value of 39 ppm for ExoTiC-JEDI, spective median transit depth precisions. None of the
and 48 ppm for Eureka!, compared to the median tran- four transmission spectra (two reductions for two vis-
sit depth uncertainty for a single visit of 34 ppm for its) shown in the upper panels of Figure 2 show any
Visit 1 and 36 ppm for Visit 2 regardless of reduction obvious features immediately identifiable by eye as ab-
method. Each reduction across a single visit are similar, sorption from any expected chemical species in this at-
with the median difference between the ExoTiC-JEDI mosphere (see §5.2). In the lower panels of Figure 2,
and Eureka! transmission spectra equal to 11 ppm for we also show the weighted average transmission spec-
 9

trum for both ExoTiC-JEDI and Eureka!, as well as the 100

joint fit transmission spectrum for Eureka!. These com- Visit 1
bined transmission spectra similarly show no obvious 75
spectral features, and demonstrate consistency between
the two reduction pipelines. In particular, the Eureka! 50
weighted average and joint transmission spectra have
a median transit depth difference of less than 5 ppm,
25
and produce transit depth precisions within 0.5 ppm of 100 3.0 3.5 4.0 4.5 5.0

Transit Depth Precision (ppm)

each other. Given this similarity, we compare the tran- Visit 2
sit depth precisions from ExoTiC-JEDI and Eureka! for 75
each visit, and in the weighted average case, to the preci-
sions predicted by PandExo (Batalha et al. 2017), shown
50
in Figure 3. While both reductions are comparable, nei- 25
ther achieves the precision predicted by PandExo across
all wavelengths, with the maximum ExoTiC-JEDI pre- 3.0 3.5 4.0 4.5 5.0
cision 1.6× and the maximum Eureka! precision 2× 60 Pandexo Prediction
the PandExo value. On average, both reductions are
ExoTiC-JEDI
Eureka!
within 1.3× the PandExo value. This is likely due to 50
the more complicated noise presented by the low num- Combined
ber of groups required for this extremely bright target
(Ngroups =3 for TOI-836), as seen in similar NIRSpec 40
programs (e.g., Lustig-Yaeger & Fu et al. 2023; Moran
& Stevenson et al. 2023; Wallack & COMPASS et al. 30
(2024)).

5. INTERPRETATION OF TOI-836B’S
20
ATMOSPHERE
To interpret TOI-836b’s atmosphere, we must first 3.0 3.5 4.0 4.5 5.0
quantify how well the transmission spectra presented in Wavelength ( m)
Figure 2 agree with each other across a variety of met-
rics. In §5.1, we perform simple synthetic fits to the Figure 3. Comparison between the transit depth precisions
data to understand: 1) how well the data is fit by a zero achieved by each reduction and the predicted values from
(i.e., flat) or non-zero sloped line, 2) potential offsets PandExo simulations. We obtain a median transit depth
between NRS1 and NRS2, and 3) whether or not these precision in wavelength bins 30 pixels wide (∼ 0.02 µm,
fits are dependent on the visit and/or reduction method. R ∼ 200) of 34 ppm for Visit 1 and 36 ppm for Visit 2, with
both the ExoTiC-JEDI weighted and Eureka! joint transmis-
These choices are driven by structure that has been seen
sion spectra resulting in a median transit depth precision of
in other small exoplanet atmospheric observations that 25 ppm. The Eureka! weighted transit depth precisions are
have complicated the overall interpretation of the atmo- indistinguishable from the joint fit precisions, with a median
sphere (e.g., May & MacDonald et al. 2023; Moran & different of less than 0.5 ppm in each wavelength bin, and
Stevenson et al. 2023), and by the fact that the pres- therefore are not visible in this plot. In all cases, the average
ence of a slope in this wavelength region can be evidence achieved precisions are within 1.3× the PandExo prediction.
of stellar activity (e.g., Rackham et al. 2018). Once we
are satisfied that we understand whether these concerns For the synthetic model fits, we use the MLFriends
may impact our conclusions regarding TOI-836b’s atmo- statistic sampler (Buchner 2016, 2019) implemented in
sphere, we can then move to more physically motivated the open source code UltraNest (Buchner 2021). For
models. In §5.2, we use PICASO models (Batalha et al. each visit and each data reduction, we fit: 1) a one-
2019; Mukherjee et al. 2023) to enable us to understand parameter, zero-slope line, 2) a two-parameter step func-
what region of parameter space we are able to effectively tion composed of two zero-sloped lines, one each for
rule out in mean molecular weight and opaque pressure NRS1 and NRS2, and 3) a two-parameter sloped line.
level. The best-fit results are shown in Table 3 and Figure 4.
The zero-slope fit for the ExoTiC-JEDI Visit 1 and
5.1. Synthetic Fits to the Data 2 data are consistent to within 1σ, resulting in a
10

Figure 4. Three synthetic fits to the data for both reductions and visits, as well as the combined visit spectra. We use three
simple models to fit the data in order to demonstrate the agreement between reductions and visits: 1) a zero-sloped line, 2) a
step function to account for offsets between NRS1 and NRS2, and 3) a non-zero-sloped line. Shaded regions illustrate the 1
and 3σ bands derived from sampling the posteriors, whereas the line represents the median best-fit model. Overall, the final
combined visit spectra are well-fit by a zero-sloped line for both reductions.
 11

transit depth baseline of (Rp /R∗ )2 = 608±3 ppm and in this context (see Table 1 Trotta 2008). This is an
609±3 ppm for Visit 1 and 2, respectively. The understandable result given that the Eureka! Visit 2
Eureka! Visit 1 data is also consistent within the 1σ offset is larger than the transit depth uncertainty near
ExoTiC-JEDI range at (Rp /R∗ )2 = 605±3, however the the detector gap.
Eureka! Visit 2 data has a somewhat higher (∼2σ) As both ExoTiC-JEDI reductions produce consistent
baseline of (Rp /R∗ )2 = 620±4. results for our synthetic modelling, we proceed with a
The step function fit for both ExoTiC-JEDI and weighted average transmission spectrum of the two visits
Eureka! Visit 1 and 2 are not consistent - Visit 1 has for the ExoTiC-JEDI reduction. As there is “moderate”
a positive step function, while Visit 2 has a negative preference for a step function offset in the Eureka! re-
step function relative to NRS2. However, comparing duction of Visit 2, we choose to leverage the joint fit
across pipelines, the ExoTiC-JEDI and Eureka! reduc- transmission spectrum from the Eureka! reduction for
tions give consistent steps within 1σ for each visit. For the rest of our analysis.
Visit 1 ExoTiC-JEDI and Eureka! produce an offset We also ran our synthetic modelling on the weighted
of +13±7 ppm and +11±7 ppm, respectively, while for average ExoTiC-JEDI and joint fit Eureka! transmis-
Visit 2, ExoTiC-JEDI and Eureka! produce an offset sion spectra to confirm that they are in agreement. In
of -18±7 ppm and -31±7 ppm, respectively. This leads the case of the joint fit, we find that the Eureka! data
to a similar discrepancy when fitting the sloped line, now obtains an offset smaller than the median tran-
where there is agreement across the reduction methods sit depth uncertainty of 25 ppm, resulting in the step
but not from visit to visit. Within 1σ, both Visit 1 and slope models no longer being preferred over the
reductions produce a positive slope, while Visit 2 pro- zero-slope model. In the case of the weighted aver-
duces a negative slope. This largely suggests that the age, the ExoTiC-JEDI data continues not to prefer the
slope and step function are not astrophysical in nature, step or slope models, with the calculated offset between
unlike those that have been seen in other observations NRS1 and NRS2 now consistent with 0 ppm. The com-
of small exoplanet atmospheres with NIRSpec/G395H bined visit transmission spectra of TOI-836b is therefore
(e.g., Moran & Stevenson et al. 2023). Regardless, the well-described by a flat line regardless of the reduction
size of the offsets obtained for both ExoTiC-JEDI visits pipeline used, and we can proceed with our physically
and for Eureka! Visit 1 are smaller than the corre- motivated modelling.
sponding median transit depth uncertainty.
To confirm that the apparent offsets between NRS1
and NRS2 need not be a major consideration in our
final interpretation of TOI-836b’s atmosphere, we can 5.2. Ruling out Physical Parameter Space
also assess which of the zero-slope, step and slope model
The spectral feature sizes of transmission spectra are
is statistically preferred by the data. Table 3 lists the
largely driven by the scale height (=kT/µg), and poten-
likelihoods for each of these fits. For both ExoTiC-JEDI
tial muting by aerosols (e.g. Sing et al. 2016). In order to
and Eureka! Visit 1 reductions, and the ExoTiC-JEDI
understand what region of parameter space we can rule
Visit 2 reduction, the step function and slope model
out for this system, we create a grid of spectral models
are not preferred or are only weakly preferred over the
as a function of metallicity and “opaque pressure level”.
zero-slope model - i.e., the data are well described by a
For the former parameter, metallicity, it is unlikely
flat line given their comparative Bayes factors (lnB12 =
that this system (R = 1.7 R⊕ ) has a large hydrogen-
lnZ1 [Model 1] − lnZ2 [Model 2]). For the Eureka! re-
helium envelope. However, similar to the analysis of
duction of the Visit 2 data, comparisons of the Bayes
the TRAPPIST-1 system (Moran et al. 2018) and of
factors suggest that both the step function and slope
other small planets (e.g. Moran & Stevenson et al. 2023;
model are at least moderately preferred over the zero
Lustig-Yaeger & Fu et al. 2023), metallicity offers a suit-
slope with lnB12 =6.5 and 2.73 , respectively. Though
able proxy for the mean molecular weight. For example,
Bayes factors do not directly map to σ-significance for
for our given assumption in temperature-pressure pro-
non-nested models (Trotta 2008), for the step function
file, 100×Solar corresponds to a mean molecular weight
and slope model these roughly translate to a strong and
of 4.3 g mol−1 , while 1000×Solar corresponds to a mean
moderate preference, respectively, over the zero slope
molecular weight of 15.7 g mol−1 . The other parameter,
“opaque pressure level”, is a term adapted from Lustig-
3 The rounded integers for Eureka! V2 in Table 3 are -70.45, - Yaeger & Fu et al. 2023 and represents a pressure below
63.94, and -67.75 for the zero-slope, step function, and slope which the atmosphere is opaque (e.g., Seager & Sas-
models respectively
selov 2000; Charbonneau et al. 2002; Berta et al. 2012;
12

Table 3. Results of synthetic fits to Visit 1 & 2 of both ExoTiC-JEDI and Eureka! data reductions as well as the combined final
transmission spectra. Each column here signifies: for the zero-slope model, the (Rp /R∗ )2 baseline intercept in ppm units, 2) for
the step function model, the offset between NRS1 and NRS2 in (Rp /R∗ )2 ppm units, and 3) for the slope case, the gradient of
the slope (ppm/µm).

Exo-TiC-JEDI (v1/v2/weighted)
2
Model Type log Z χ /N Baseline Intercept NRS1/NRS2 Offset Slope Gradient
Zero slope -62/-66/-68 1.10/1.16/1.21 608±3.3/609±3.5/608±2.5 N/A N/A
Step Function -64/-65/-72 1.07/1.10/1.21 N/A +13±7/-18±7/-2±5 N/A
Slope -64/-67/-72 1.06/1.14/1.21 N/A N/A 13±6/-10±6/2±5
Eureka! (v1/v2/joint)
2
Model Type log Z χ /N Baseline Intercept NRS1/NRS2 Offset Slope Gradient
Zero-slope -64/-70/-71 1.13/1.25/1.28 605±3.4/620±4/614±2.4 N/A N/A
Step Function -66/-64/-73 1.11/1.07/1.25 N/A +11±7/-31±7/-6±5 N/A
Slope -67/-68/-75 1.10/1.13/1.25 N/A N/A 10±6/-23±6/-4±4

(a) (b)

Figure 5. (a) For a single choice in opaque pressure level (0.1 bar) we show the parameter space that can be ruled out in
metallicity. Blue lines show the reductions for ExoTiC-JEDI (Visit 1, 2, and weighted) and orange lines show the reductions for
Eureka! (Visit 1, 2, and joint). The black-dashed line indicates the 3σ level, below which we are unable to confidently rule out
models. Ultimately our data rules out metallicities < 250×Solar, corresponding to a mean molecular weight of ∼ 6 g mol−1 . (b)
For four of the metallicity cases shown in (a), we show the spectra relative to the weighted data from ExoTiC-JEDI. We also
indicate the χ2 /N and σ for reference.

Kreidberg et al. 2014; Knutson et al. 2014)4 . Other els (100–10−4 bar, log-spaced with 5 grid points), we
manuscripts have referred to this as a “cloud top pres- compute a grid of transmission spectra using the open
sure” (e.g., Kreidberg et al. 2014; Moran et al. 2018). source PICASO package (Batalha et al. 2019). For the
Here we use a more general term as we cannot differen- pressure-temperature profile, we use a 1D 5-parameter
tiate between a cloud top pressure from a surface pres- double-grey analytic formula (Guillot 2010). We also
sure. test whether or not the results are sensitive to our choice
Using a range of metallicities (1–1000×Solar, log of pressure-temperature profile by computing a similar
spaced with 26 grid points) and opaque pressure lev- grid with simple isothermal pressure-temperature pro-
files. Our conclusions do not change depending on this
4
choice. Given the pressure-temperature profile, we fix
Lustig-Yaeger & Fu et al. 2023 used the term “apparent surface
pressure” but it has the same meaning the elemental ratio C/O to solar (=0.55 Asplund et al.
2009) and obtain the chemistry by interpolating on a
 13

pre-computed chemical equilibrium grid. The chemistry sphere, in which the atmosphere only becomes opaque
grid was computed by Line et al. (2013) with NASA’s below 100 bar, the combined ExoTiC-JEDI spectrum is
CEA code (Gordon & McBride 1994). This grid is pub- able to rule out the 300×Solar case. For cases where the
licly available on GitHub as part of CHIMERA’s open opaque pressure level is 10−4 bar, comparable to a highly
source code5 . Of note to this analysis, the molecules lofted cloud, we can rule out cases ≤100×Solar. Both
which absorb from 3–5 µm are H2 O, CH4 , CO2 , and the 100 bar and 10−4 bar cases represent unlikely phys-
CO, which are all included in the grid. ical scenarios, as clouds are expected to form in super-
Figure 5 shows the results of the grid analysis, where Earth atmospheres (Mbarek & Kempton 2016) (i.e., at-
we show how our confidence level, expressed as a σ- mospheres are highly unlikely to be effectively cloud
level, changes as a function of metallicity for an in- free), and clouds lofted to high altitudes are unlikely
termediate opaque pressure level of 0.1 bar. The gen- to be completely opaque (Robinson & Catling 2014),
eral behaviour of the significance curves in Figure 5a, however as end-member cases they demonstrate that
which peak toward 10-50×Solar and decrease toward solar-like composition atmospheres are not plausible for
1000×Solar, is well-documented (e.g. Moran et al. 2018). TOI-836b regardless of the height of any opaque pres-
Toward 10×Solar the magnitude of spectral features in- sure level. Additionally, we tested whether or not our
creases because the added molecular opacity is better conclusions are affected by the choice of binning scheme
able to surpass the contribution from H2 /He continuum and found that the conclusions are unchanged.
without affecting the mean molecular weight, increas- Figure 5b shows four of the spectra used to compute
ing the significance to which the spectral features can the σ-significance curves in 5a, for reference, along with
be ruled out. Beyond ∼50×Solar the contribution from the weighted average spectra from ExoTiC-JEDI. The
the heavier metals starts to increase the mean molecular main features shown are that of CH4 and CO2 , in the
weight of the atmosphere and results in overall smaller 1×Solar case, and primarily H2 O and CO2 in the other
spectral features which are harder to more confidently cases. Figure 5b also lists the χ2 /N and sigma rejection
rule out. thresholds for each of the metallicity cases, demonstrat-
We choose to show 0.1 bar for reference as it is syn- ing, for example, that we cannot confidently distinguish
onymous with the tropopause of all Solar System objects between atmospheres with 250×Solar and 1000×Solar
(Robinson & Catling 2014). Here, σ is computed by con- metallicities given the constraints we achieve with two
verting χ2 /N to a p-value, and then to a σ-significance. combined transit observations.
For each individual visit, we are able to rule out metal-
licities lower than 100–160×Solar depending on the visit 5.3. Theoretical predictions of possible interiors of
and the reduction. For example, in the case of the first TOI-836b
visit, the individual ExoTiC-JEDI and Eureka! spectra
TOI-836b’s mass and radius place it at a very in-
enable metallicities to be ruled out at the <130×Solar
triguing position in the mass-radius diagram towards
and <160×Solar-level, respectively. With the combined
the lower edge of the radius valley. Furthermore, de-
spectra, we are further able to rule out metallicities
spite its low density, its parameters are compatible with
< 250×Solar and < 380×Solar, for ExoTiC-JEDI and
a pure rock composition (no iron core) at the 1σ level.
Eureka!, respectively, corresponding to mean molecu-
To infer the bulk properties of TOI-836b, we use the
lar weights of ∼6–9 g mol−1 . Combining both visits en-
SMINT (Structure Model INTerpolator) package from Pi-
ables us to rule out nearly double the parameter space
aulet et al. (2021), which performs an MCMC retrieval
in metallicity, demonstrating the potential of multi-visit
of planetary bulk compositions from pre-computed grids
observing strategies. As shown in Figure 5a, the re-
of theoretical interior structure models. We consider two
ported 3σ lower limit on metallicity is dependent on the
possible compositions for the interior: 1) an Earth-like
reduction method. However, the overall scientific con-
core with a H2 -He envelope of solar metallicity (Lopez
clusions are agnostic to the data reduction, as in all cases
& Fortney 2014), and 2) a refractory core with a vari-
we are able to rule out H2 -dominated atmospheres with
able core mass fraction and a pure H2 O envelope and
mean molecular weights less than ∼6 g mol−1 .
atmosphere on top (Aguichine et al. 2021). These com-
Figure 5 shows the results for an intermediate cloud
positions represent end-member cases between an en-
case, though we ran a full grid of both cloud-free and
velope that would form by accreting nebular gas with a
highly cloudy cases. For an effective cloud-free atmo-
Sun-like composition, and a high mean molecular weight
envelope where water is used as a proxy for all volatiles.
5 https://github.com/mrline/CHIMERA Based on its bulk properties and these scenarios, we find
that TOI-836b could have an envelope mass fraction of
14

into very different bulk compositions when inferred from

interior structure models. The interior of TOI-836b is
compatible with a pure rock (no iron core) composition
at 1σ, meaning that the possibility that TOI-836b is a
terrestrial planet cannot be excluded. From our interior
modelling, we also find that the planet could be made
of at most 0.1% solar metallicity gas or 9 ± 5% pure
water. In contrast, the possible bulk compositions of
TOI-836c are 1.74+0.55
−0.48 % in the case of solar metallic-
ity gas, or 52+15
−14 % in the pure water case (Wallack &
COMPASS et al. 2024). Figure 6 shows these possible
interior compositions for both planets, which represent
end-member cases for hydrogen-dominated and pure wa-
ter atmospheres such that intermediate compositions are
also possible.
Given TOI-836b’s high equilibrium temperature and
the stellar insolation flux received, photoevaporation
could likely be responsible for the observed lack of a
hydrogen-dominated atmosphere. Applying the photo-
evaporation model of Rogers et al. (2021) to a planet
with the properties of TOI-836b, and including a core
mass of 4.5 M⊕ , we find that any initial envelope mass
fraction in the range 2–30% is blown away in <400 Myr.
Applying the model to a planet with the properties of
TOI-836c, and including a core mass of 9.4 M⊕ , we
find that hydrogen envelopes of up to 10% can be re-
tained. Given the reported age of TOI-836 of 5.4+6.3
−5.0 Gyr
(Hawthorn et al. 2023), this analysis strongly suggests
the absence of a hydrogen-dominated atmosphere for
TOI-836b, which is in line with the apparent >6 g mol−1
mean molecular weight derived in §5.2. The bulk compo-
Figure 6. Mass-radius plots for the population of small ex- sition of TOI-836c, however, is still degenerate, with low
oplanets demonstrating possible interior compositions as cal- mean molecular weight atmospheres still plausible, par-
culated by SMINT, where the colour of each marker represents ticularly in the presence of clouds and hazes (Wallack &
either the bulk H2 -He (top) or bulk water (bottom) mass
COMPASS et al. 2024). These models provide the first
fractions. The planets in the TOI-836 system are denoted by
star-shaped markers. Density curves for the Earth-like and clues as to this system’s possible evolution, although we
50 wt% liquid water compositions (Zeng & Jacobsen 2016); caution that transmission spectra alone are unable to de-
50% and 100% steam atmospheres assuming Teq = 600 K termine whether the TOI-836 planets formed with their
and an Earth-like core (Aguichine et al. 2021); and 0.1%, current masses, or if the present-day difference in bulk
1%, 2% and 5% H2 -He composition assuming an age of 5 Gyr compositions is the consequence of photoevaporation.
(Lopez & Fortney 2014) are also plotted for reference. The
background planet population is obtained from the NASA
Exoplanet Archive. 6.2. Implications for Future JWST Observations
Here we examine how our results can inform the plan-
at most 0.1% in the case of gas of solar composition, or ning of future observations of small planets with high
a water mass fraction of 9 ± 5% in the pure H2 O case as atmospheric metallicities that orbit bright stars, com-
shown in Figure 6. paring our measured data to PandExo JWST simula-
tions, which are used by the community for planning
6. DISCUSSION
observations. Using the grid of models described in §5.2
we determine how many additional transits would be
6.1. The TOI-836 System needed to rule out a certain metallicity model. This is an
The two planets of the TOI-836 system are located identical exercise to that performed in Figure 5, except
on opposite sides of the radius valley, which translates here our “data” is a PandExo simulation of a featureless
 15

future observation planning of planets that are expected

to be heavily enriched in metals, particularly for those
around bright stars with more complex noise properties.

7. CONCLUSIONS
We have presented two JWST NIRSpec/G395H obser-
vations of the transmission spectrum of the super-Earth
TOI-836b. We produce two reductions of the data with
independent pipelines, ExoTiC-JEDI and Eureka!, re-
sulting in a median transit depth uncertainty for both
methods of 34 ppm for Visit 1 and 36 ppm for Visit 2
in 30 pixel wide bins. We combine our two visits using
a weighted average for ExoTiC-JEDI and a joint fit for
Eureka!, and find a combined median transit depth pre-
cision of 25 ppm in both cases. We also find sub-ppm dif-
ferences in the precision obtained by the Eureka! joint
fit and a weighted average of the individual Eureka!
visits at all wavelengths.
Figure 7. The number of transits needed to rule out a zero-
When modelling our transmission spectra, we find
sloped line at 3σ as a function of metallicity for an opaque
pressure level of 0.1 bar. Here, “real data” uses the preci-
that transmission spectra that appear to be flat “by-
sion derived from the first visit of the Eureka! reduction eye” can have different retrieved transit depth baselines
and assumes that additional transits improve the precision and detector offsets. We caution that these model pa-
based on the measured improvement from a single visit to rameterisations are a simple and basic test to determine
the Eureka! joint fit. The PandExo data curve follows simu- first-order structures in the data, but are not necessarily
lations computed with the JWST simulation tool PandExo. an accurate representation of the intrinsic scatter across
the transmission spectrum. Future work will be needed
spectrum. We compute noise simulations with PandExo as more data is collected to better characterise the noise
(Batalha et al. 2017) using an identical observational properties being seen in JWST observations. Careful
setup to our program here, assuming that √ each addi- analyses of each visit and each data reduction method
tional visit provides a gain in precision of n transit. should therefore be done individually and assessed col-
For the “real data” case, we use the noise budget mea- lectively, even if the data appear to be overall consistent
sured in this program from Visit 1, assuming that the to within 1σ.
increase in precision is equivalent to what we have mea- Our final combined transmission spectrum from each
sured moving from individual visits to the combined reduction method is well described by a flat line, with
joint Eureka! reduction. Doing √ so results in a mea- no obvious atmospheric features. PICASO modelling en-
sured precision gain of ∼98% 2 for TOI-836b. Using ables us to rule out atmospheres of at least <100×Solar
each of the estimates for noise, we compute simulated metallicity regardless of the height of an opaque pres-
observations of each modelled spectrum and then com- sure level (equivalent to either a cloud deck or surface).
pute the number of transits needed to rule out a zero- With our combined two visit spectra for the 0.1 bar
sloped line with 3σ confidence. Additionally, we include case, we specifically rule out <250×Solar metallicities
random noise and repeat the test for 1000 different ran- for the ExoTiC-JEDI spectrum and <380×Solar metal-
dom noise instances. Figure 7 shows the median result licities for the Eureka! spectrum. These constraints
for the case of an opaque pressure level of 0.1 bar. allow us to rule out atmospheres with mean molecu-
Overall, ruling out cases up to 1000×Solar metallic- lar weights less than ∼ 6 g mol−1 . Given the modelling
ity for TOI-836b would require up to eight additional setup considered in this work, combining both visits en-
transits, assuming that the√ data continued to result in ables us to rule out nearly double the parameter space
a precision gain of ∼98% n transit. For lower metallic- in metallicity when compared to each visit individually.
ity cases (<100×Solar), PandExo predictions are in line Comparing our mean molecular weight for TOI-836b
with those based on the real data. However, for higher to interior and photoevaporation evolution modelling
metallicities, PandExo data appears to be somewhat op- strongly supports our overall conclusion that this super-
timistic, resulting in estimates requiring 1-2 fewer tran- Earth does not possess a H2 -dominated atmosphere, in
sits than the real data. This result should be noted for possible contrast to the larger, exterior TOI-836c.
16

As JWST continues to observe small planets, we rec- work supported by NASA’S Interdisciplinary Consor-
ommend that care is taken when using simulation tools tia for Astrobiology Research (NNH19ZDA001N-ICAR)
to determine how many transits may be needed to rule under award number 19-ICAR19 2-0041. This work
out certain physical scenarios, particularly in the case benefited from the 2023 Exoplanet Summer Program
of observations that require a small number of groups. in the Other Worlds Laboratory (OWL) at the Univer-
For high-metallicity atmospheres (>100×Solar), we sity of California, Santa Cruz, a program funded by the
find that PandExo predictions are optimistic com- Heising-Simons Foundation. This research also made
pared to the precision gains from our measured data, use of the NASA Exoplanet Archive, which is operated
and yield estimates with 1–2 transits less than may by the California Institute of Technology, under contract
be required. This should be accounted for in future with the National Aeronautics and Space Administra-
observation proposals of small planets with JWST. tion under the Exoplanet Exploration Program.
Co-Author contributions are as follows: LA led the
We thank the anonymous referee for their comments data analysis and write-up of this study. NEB led the
that helped improve the quality and clarity of this pa- modelling efforts. NLW and JIAR provided reductions
per. The data products for this manuscript can be and analyses of the data. AA performed the interior
found at the following Zenodo repository: 10.5281/zen- modelling. HRW advised throughout the analysis and
odo.10658637. L.A. would like to thank D. Grant and manuscript preparation. All authors read and provided
M. Radica for useful discussions regarding JWST data comments and conversations that greatly improved the
analysis. This work is based on observations made with quality of the manuscript.
the NASA/ESA/CSA James Webb Space Telescope.
Software: astropy (Astropy Collaboration et al.
The data were obtained from the Mikulski Archive for
2013, 2018, 2022), batman (Kreidberg 2015), emcee
Space Telescopes at the Space Telescope Science Insti-
(Foreman-Mackey et al. 2013), Eureka! (Bell et al.
tute, which is operated by the Association of Universi-
2022), ExoTiC-Jedi (Alderson et al. 2022), ExoTiC-LD
ties for Research in Astronomy, Inc., under NASA con-
(Grant & Wakeford 2022), Matplotlib (Caswell et al.
tract NAS 5-03127 for JWST. These observations are
2019), NumPy (Oliphant 2006), pandas (McKinney et al.
associated with program #2512. Support for program
2011), PICASO (Batalha et al. 2018; Mukherjee et al.
#2512 was provided by NASA through a grant from the
2023), PandExo (Batalha et al. 2017), scipy (Virtanen
Space Telescope Science Institute, which is operated by
et al. 2019), SMINT (Piaulet et al. 2021), STScI JWST
the Association of Universities for Research in Astron-
Calibration Pipeline (Bushouse et al. 2022), ultranest
omy, Inc., under NASA contract NAS 5-03127. L.A.
(Buchner 2021), xarray (Hoyer & Hamman 2017)
acknowledges funding from STFC grant ST/W507337/1
and from the University of Bristol School of Physics PhD
Scholarship Fund. This work is funded in part by the Al- Facilities: JWST (NIRSpec)
fred P. Sloan Foundation under grant G202114194. Sup- The JWST data presented in this paper were ob-
port for this work was provided by NASA through grant tained from the Mikulski Archive for Space Telescopes
80NSSC19K0290 to J.T. and N.W. H.R.W. was funded (MAST) at the Space Telescope Science Institute. The
by UK Research and Innovation (UKRI) under the UK specific observations analysed can be accessed via DOI:
government’s Horizon Europe funding guarantee [grant 10.17909/fpwt-rn60.
number EP/Y006313/1]. This material is based upon

REFERENCES
Aguichine, A., Mousis, O., Deleuil, M., & Marcq, E. 2021, Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009,
ApJ, 914, 84, doi: 10.3847/1538-4357/abfa99 ARA&A, 47, 481,
Ahrer, E., Wheatley, P. J., Gandhi, S., et al. 2023a, doi: 10.1146/annurev.astro.46.060407.145222
MNRAS, 521, 5636, doi: 10.1093/mnras/stad779 Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,
Ahrer, E.-M., Stevenson, K. B., Mansfield, M., et al. 2023b, et al. 2013, A&A, 558, A33,
Nature, 614, 653, doi: 10.1038/s41586-022-05590-4 doi: 10.1051/0004-6361/201322068
Alderson, L., Wakeford, H. R., MacDonald, R. J., et al. Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M.,
2022, MNRAS, 512, 4185, doi: 10.1093/mnras/stac661 et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
Alderson, L., Wakeford, H. R., Alam, M. K., et al. 2023, Astropy Collaboration, Price-Whelan, A. M., Lim, P. L.,
Nature, 614, 664, doi: 10.1038/s41586-022-05591-3 et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
 17

Batalha, N. E., Lewis, N. K., Line, M. R., Valenti, J., & Figueira, P., Pont, F., Mordasini, C., et al. 2009, A&A, 493,
Stevenson, K. 2018, ApJL, 856, L34, 671, doi: 10.1051/0004-6361:20078951
doi: 10.3847/2041-8213/aab896 Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman,
Batalha, N. E., Marley, M. S., Lewis, N. K., & Fortney, J. 2013, Publications of the Astronomical Society of the
J. J. 2019, ApJ, 878, 70, doi: 10.3847/1538-4357/ab1b51 Pacific, 125, 306, doi: 10.1086/670067
Batalha, N. E., Wolfgang, A., Teske, J., et al. 2023, AJ, Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013,
165, 14, doi: 10.3847/1538-3881/ac9f45 ApJ, 775, 80, doi: 10.1088/0004-637X/775/1/80
Batalha, N. E., Mandell, A., Pontoppidan, K., et al. 2017, Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017,
PASP, 129, 064501, doi: 10.1088/1538-3873/aa65b0 AJ, 154, 109, doi: 10.3847/1538-3881/aa80eb
Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, Ginzburg, S., Schlichting, H. E., & Sari, R. 2018, MNRAS,
ApJS, 204, 24, doi: 10.1088/0067-0049/204/2/24 476, 759, doi: 10.1093/mnras/sty290
Bell, T., Ahrer, E.-M., Brande, J., et al. 2022, The Journal Gordon, S., & McBride, B. J. 1994, Computer program for
of Open Source Software, 7, 4503, calculation of complex chemical equilibrium compositions
doi: 10.21105/joss.04503 and applications. Part 1: Analysis, Tech. rep.
Berta, Z. K., Charbonneau, D., Désert, J.-M., et al. 2012, Grant, D., & Wakeford, H. R. 2022, Exo-TiC/ExoTiC-LD:
ApJ, 747, 35, doi: 10.1088/0004-637X/747/1/35 ExoTiC-LD v3.0.0, v3.0.0, Zenodo, Zenodo,
Birkmann, S. M., Giardino, G., Sirianni, M., et al. 2022, in doi: 10.5281/zenodo.7437681
Society of Photo-Optical Instrumentation Engineers Guillot, T. 2010, A&A, 520, A27,
(SPIE) Conference Series, Vol. 12180, Space Telescopes
doi: 10.1051/0004-6361/200913396
and Instrumentation 2022: Optical, Infrared, and
Hawthorn, F., Bayliss, D., Wilson, T. G., et al. 2023,
Millimeter Wave, ed. L. E. Coyle, S. Matsuura, & M. D.
MNRAS, 520, 3649, doi: 10.1093/mnras/stad306
Perrin, 121802P, doi: 10.1117/12.2629545
Hoyer, S., & Hamman, J. 2017, Journal of Open Research
Buchhave, L. A., Bizzarro, M., Latham, D. W., et al. 2014,
Software, 5
Nature, 509, 593, doi: 10.1038/nature13254
Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022,
Buchner, J. 2016, Statistics and Computing, 26, 383,
A&A, 661, A80, doi: 10.1051/0004-6361/202142663
doi: 10.1007/s11222-014-9512-y
Kirk, J., Wheatley, P. J., Louden, T., et al. 2018, MNRAS,
—. 2019, PASP, 131, 108005,
474, 876, doi: 10.1093/mnras/stx2826
doi: 10.1088/1538-3873/aae7fc
Knutson, H. A., Dragomir, D., Kreidberg, L., et al. 2014,
—. 2021, The Journal of Open Source Software, 6, 3001,
ApJ, 794, 155, doi: 10.1088/0004-637X/794/2/155
doi: 10.21105/joss.03001
Kostogryz, N., Shapiro, A., Witzke, V., et al. 2023,
Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022,
Research Notes of the AAS, 7, 39
JWST Calibration Pipeline, 1.8.2, Zenodo, Zenodo,
Kostogryz, N., Witzke, V., Shapiro, A., et al. 2022, arXiv
doi: 10.5281/zenodo.7229890
preprint arXiv:2206.06641
Caswell, T. A., Droettboom, M., Hunter, J., et al. 2019,
matplotlib/matplotlib v3.1.0, v3.1.0, Zenodo, Kreidberg, L. 2015, PASP, 127, 1161, doi: 10.1086/683602
doi: 10.5281/zenodo.2893252 Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014,
Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, Nature, 505, 69–72, doi: 10.1038/nature12888
R. L. 2002, ApJ, 568, 377 Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, ApJL,
Claret, A. 2000, A&A, 363, 1081 793, L27, doi: 10.1088/2041-8205/793/2/L27
Crossfield, I. J. M., Barman, T., Hansen, B. M. S., & Libby-Roberts, J. E., Berta-Thompson, Z. K.,
Howard, A. W. 2013, A&A, 559, A33, Diamond-Lowe, H., et al. 2022, AJ, 164, 59,
doi: 10.1051/0004-6361/201322278 doi: 10.3847/1538-3881/ac75de
Diamond-Lowe, H., Charbonneau, D., Malik, M., Kempton, Line, M. R., Knutson, H., Deming, D., Wilkins, A., &
E. M. R., & Beletsky, Y. 2020, AJ, 160, 188, Desert, J.-M. 2013, ApJ, 778, 183,
doi: 10.3847/1538-3881/abaf4f doi: 10.1088/0004-637X/778/2/183
Diamond-Lowe, H., Mendonça, J. M., Charbonneau, D., & Lopez, E. D., & Fortney, J. J. 2013, ApJ, 776, 2,
Buchhave, L. A. 2023, AJ, 165, 169, doi: 10.1088/0004-637X/776/1/2
doi: 10.3847/1538-3881/acbf39 —. 2014, ApJ, 792, 1, doi: 10.1088/0004-637X/792/1/1
Feinstein, A. D., Radica, M., Welbanks, L., et al. 2023, Lopez, E. D., Fortney, J. J., & Miller, N. 2012, ApJ, 761,
Nature, 614, 670, doi: 10.1038/s41586-022-05674-1 59, doi: 10.1088/0004-637X/761/1/59
18

Louden, T., Wheatley, P. J., Irwin, P. G., Kirk, J., & Rackham, B. V., Apai, D., & Giampapa, M. S. 2018, ApJ,
Skillen, I. 2017, MNRAS, 470, 742 853, 122
Lustig-Yaeger, J., Fu, G., May, E. M., et al. 2023, arXiv Robinson, T. D., & Catling, D. C. 2014, Nature Geoscience,
e-prints, arXiv:2301.04191, 7, 12, doi: 10.1038/ngeo2020
doi: 10.48550/arXiv.2301.04191 Rogers, J. G., Gupta, A., Owen, J. E., & Schlichting, H. E.
2021, MNRAS, 508, 5886, doi: 10.1093/mnras/stab2897
May, E. M., MacDonald, R. J., Bennett, K. A., et al. 2023,
Rogers, L. A. 2015, ApJ, 801, 41,
arXiv e-prints, arXiv:2310.10711,
doi: 10.1088/0004-637X/801/1/41
doi: 10.48550/arXiv.2310.10711
Rogers, L. A., Bodenheimer, P., Lissauer, J. J., & Seager,
Mbarek, R., & Kempton, E. M. R. 2016, ApJ, 827, 121,
S. 2011, ApJ, 738, 59, doi: 10.1088/0004-637X/738/1/59
doi: 10.3847/0004-637X/827/2/121 Rogers, L. A., & Seager, S. 2010, ApJ, 712, 974,
McKinney, W., et al. 2011, Python for high performance doi: 10.1088/0004-637X/712/2/974
and scientific computing, 14, 1 Rustamkulov, Z., Sing, D. K., Mukherjee, S., et al. 2023,
Moran, S. E., Hörst, S. M., Batalha, N. E., Lewis, N. K., & Nature, 614, 659, doi: 10.1038/s41586-022-05677-y
Wakeford, H. R. 2018, AJ, 156, 252, Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916
doi: 10.3847/1538-3881/aae83a Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature,
Moran, S. E., Stevenson, K. B., Sing, D. K., et al. 2023, 529, 59
ApJL, 948, L11, doi: 10.3847/2041-8213/accb9c Teske, J., Wang, S. X., Wolfgang, A., et al. 2021, ApJS,
256, 33, doi: 10.3847/1538-4365/ac0f0a
Mukherjee, S., Batalha, N. E., Fortney, J. J., & Marley,
M. S. 2023, ApJ, 942, 71, doi: 10.3847/1538-4357/ac9f48 Trotta, R. 2008, Contemporary Physics, 49, 71,
doi: 10.1080/00107510802066753
Oliphant, T. 2006, NumPy: A guide to NumPy, USA:
Virtanen, P., Gommers, R., Burovski, E., et al. 2019,
Trelgol Publishing. http://www.numpy.org/
scipy/scipy: SciPy 1.2.1, v1.2.1, Zenodo,
Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931,
doi: 10.5281/zenodo.2560881
doi: 10.1111/j.1365-2966.2012.21481.x Wallack, N., & COMPASS et al. 2024, AAS Journals,
Owen, J. E., & Schlichting, H. E. 2023, arXiv e-prints, submitted
arXiv:2308.00020, doi: 10.48550/arXiv.2308.00020 Wordsworth, R., & Kreidberg, L. 2022, ARA&A, 60, 159,
Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105, doi: 10.1146/annurev-astro-052920-125632
doi: 10.1088/0004-637X/775/2/105 Zeng, L., & Jacobsen, S. B. 2016, ApJ, 829, 18,
Piaulet, C., Benneke, B., Rubenzahl, R. A., et al. 2021, AJ, doi: 10.3847/0004-637X/829/1/18
161, 70, doi: 10.3847/1538-3881/abcd3c Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019,
Pont, F., Zucker, S., & Queloz, D. 2006, MNRAS, 373, 231, Proceedings of the National Academy of Science, 116,
9723, doi: 10.1073/pnas.1812905116
doi: 10.1111/j.1365-2966.2006.11012.x