Atmospheric Reconnaissance of TRAPPIST-1 b with JWST/NIRISS: Evidence for Strong Stellar Contamination in the Transmission Spectra

The stellar-contamination model assumes that the visible stellar disk has two unocculted regions, a spot and a facula, for a total of three components, including the quiescent photosphere. We opted for three components as it is the simplest model that allows one component to be hotter and one colder than the quiescent photosphere. The covering fractions of the spot and facula can go to zero, hence allowing a two- and a single-component scenario. We did not compare our results with models with more components, as the reduced χ² values were already close to 1, i.e., the three-component model already provided a reasonable fit to the observations (Section 6.1). However, we do acknowledge that the transit depth precision may be underestimated because it does not fully account for the uncertainty on the number of components. More analysis, and likely more observations, will be needed to determine the number of components that best fit the data.

We compute a weighted average of three PHOENIX ACES stellar spectra (Husser et al. 2013) at the temperatures and gravities of each component to simulate the spectrum of a heterogeneous disk. The weight of each spectrum is determined by the size of the spot and of the facula. We then compute the expected transmission spectrum assuming an atmosphere-less planet transiting this heterogeneous disk (Equation (3) from Rackham et al. 2018).

Since the two visits were approximately 1.5 days apart (Gillon et al. 2017) and adopting a stellar rotation period of 3.3 days (Luger et al. 2017), the two visits were about half a stellar rotation apart, meaning that the spot and facula properties can be different between these two visits. We thus fitted the stellar-contamination model to the transit spectrum of each visit separately, inferring the properties of the unocculted stellar heterogeneities, i.e., their temperature differences with respect to the quiescent photosphere, surface gravities, and covering fractions, with the priors listed in Table 1. We forced the sum of both covering fractions to remain smaller than 1. Following Fournier-Tondreau et al. (2023), we set the surface gravity as a free parameter as a proxy of the pressure in the heterogeneity, which may be affected by the presence of magnetic fields (Solanki 2003; Bruno et al. 2021). For simplicity, we forced both heterogeneities to have the same surface gravity in this particular fit. We also set the planet-to-star radius ratio and the photospheric temperature as free parameters to account for their uncertainties.

Table 1. Priors and Posteriors of the Stellar-contamination Parameters from the Sequential and Simultaneous Fits of the Stellar Contamination and Planetary Atmosphere

Stellar-contamination Parameter	Priors		Posteriors
	Sequential	Simultaneous	Sequential		Simultaneous
	Visits 1 & 2	Visits 1 & 2	Visit 1	Visit 2	Visit 1	Visit 2
Spot covering fraction, fspot	${ \mathcal U }(0,1)$	${ \mathcal U }(0,0.5)$	${0.29}_{-0.10}^{+0.11}$	${0.16}_{-0.08}^{+0.11}$	${0.23}_{-0.08}^{+0.10}$	${0.24}_{-0.14}^{+0.15}$
Facula covering fraction, ffac	${ \mathcal U }(0,1)$	${ \mathcal U }(0,0.5)$	${0.05}_{-0.04}^{+0.08}$	${0.29}_{-0.13}^{+0.12}$	${0.07}_{-0.04}^{+0.09}$	${0.27}_{-0.18}^{+0.15}$
Photospheric temperature, Tphot (K)	${ \mathcal N }({T}_{\star ,\mathrm{eff}},{\sigma }_{{T}_{\star ,\mathrm{eff}}})$	${ \mathcal N }({T}_{\star ,\mathrm{eff}},{\sigma }_{{T}_{\star ,\mathrm{eff}}})$	${2567}_{-46}^{+52}$	${2557}_{-45}^{+46}$	${2571}_{-42}^{+43}$	${2576}_{-37}^{+34}$
Stellar photosphere gravity, $\mathrm{log}{g}_{\mathrm{phot}}$ (dex)	5.0	${ \mathcal N }(\mathrm{log}{g}_{\star },{\sigma }_{\mathrm{log}{g}_{\star }})$	5.0	5.0	${5.24}_{-0.01}^{+0.01}$	${5.24}_{-0.01}^{+0.01}$
Spot temperature, Tspot or ΔTspot (K)	${ \mathcal U }(-266,-100)$ a	${ \mathcal U }(2300,{T}_{\star ,\mathrm{eff}}+3\,{\sigma }_{{T}_{\star ,\mathrm{eff}}})$ b	$-{197}_{-62}^{+56}$ a	$-{178}_{-60}^{+54}$ a	${2371}_{-49}^{+70}$ b	${2431}_{-89}^{+83}$ b
Facula temperature, Tfac or ΔTfac (K)	${ \mathcal U }(100,1000)$ a	${ \mathcal U }({T}_{\star ,\mathrm{eff}}-3\,{\sigma }_{{T}_{\star ,\mathrm{eff}}},1.2\,{T}_{\star ,\mathrm{eff}})$ b	$+{214}_{-88}^{+257}$ a	$+{153}_{-33}^{+61}$ a	${2706}_{-86}^{+124}$ b	${2734}_{-53}^{+105}$ b
Spot surface gravity, $\mathrm{log}{g}_{\mathrm{spot}}$ (dex)	${ \mathcal U }(2.5,5.5)$	${ \mathcal U }(3.0,5.4)$	${4.3}_{-0.5}^{+0.2}$	${5.22}_{-0.12}^{+0.11}$	${4.55}_{-0.42}^{+0.57}$	${5.33}_{-0.15}^{+0.05}$
Facula surface gravity, $\mathrm{log}{g}_{\mathrm{fac}}$ (dex)		${ \mathcal U }(3.0,5.4)$			${4.49}_{-0.81}^{+0.57}$	${5.28}_{-0.25}^{+0.09}$
Planet radius Rp (R⊕)	...	${ \mathcal U }(0.85,1.15){R}_{{\rm{p}},\mathrm{obs}}/{R}_{\oplus }$	${1.09}_{-0.01}^{+0.01}$	${1.13}_{-0.01}^{+0.01}$	${1.10}_{-0.01}^{+0.01}$	${1.12}_{-0.01}^{+0.01}$
Planet-to-star radius ratio, squared ${({R}_{{\rm{p}}}/{R}_{\star })}^{2}$ (ppm)	${ \mathcal U }(0.8,1.2)\mathrm{med}({T}_{\mathrm{depth}})$	...	${6951}_{-142}^{+144}$	${7397}_{-150}^{+148}$	${7062}_{-128}^{+128}$	${7321}_{-131}^{+131}$

Notes. In the prior columns, T_⋆,eff = 2566 K, $\mathrm{log}{g}_{\star }=5.24$, and R_p,obs = 1.116 R_⊕ are adopted from the literature (Agol et al. 2021), with slightly more conservative uncertainties ${\sigma }_{{T}_{\star ,\mathrm{eff}}}=50\,{\rm{K}}$ and ${\sigma }_{\mathrm{log}{g}_{\star }}=0.01$, and med(T_depth) is the median of the measured transit depths of each visit. In the posterior columns, the leftmost value is the 50th percentile of the posterior distribution. Upper and lower uncertainties are computed from the 16th, 50th, and 84th percentiles.

^a Temperature difference with respect to the quiet photosphere. ^b Absolute temperature, not relative to the quiet photosphere.

Download table as:  ASCIITypeset image

We then divided out the median stellar-contamination model from each of the two transmission spectra and multiplied back the median squared planet-to-star radius ratio. We took an error-weighted average of the two stellar-contamination-corrected transmission spectra, yielding a stellar-contamination-corrected and combined transmission spectrum.

Before performing the planetary atmosphere retrieval, we first compared the stellar-contamination-corrected and combined transmission spectrum to a selection of SCARLET planetary atmosphere forward models (Benneke & Seager 2012, 2013; Benneke 2015, 2019; Benneke et al. 2019; Pelletier et al. 2021) of TRAPPIST-1 b. We computed the χ² statistic between each atmospheric model and the corrected and combined spectrum, and measured how far out in the distribution tail the statistic fell, in units of the standard deviation of the χ² distribution (Gregory 2005). This frequentist calculation provides us with rejection (or detection) confidence levels for specific planetary atmosphere scenarios, which can guide the atmospheric retrieval.

We used SCARLET to perform a planetary atmosphere retrieval on the stellar-contamination-corrected and combined transmission spectrum. Since the frequentist comparison between atmospheric forward models and the corrected and combined spectrum did not detect an atmosphere, this retrieval analysis explored the range of mean molecular masses and surface/cloud pressures that remain consistent with the data, similar to an analysis done with GJ 1214 b and GJ 436 b (Benneke & Seager 2013; Kreidberg et al. 2014; Knutson et al. 2014). We used a parameter space that efficiently samples the full range of models between hydrogen-dominated and volatile-dominated atmospheres. We adopted a mix of H₂/He at solar ratio, along with molecules of higher mean molecular mass, namely H₂O, CH₄, CO, CO₂, NH₃, and N₂. This is because TRAPPIST-1 planets could host a wide variety of secondary atmospheres with outgassing/volcanism and cometary delivery as their main sources (Turbet et al. 2020 and references therein), along with photochemistry shaping the composition of these atmospheres given the strong X-ray and UV radiation and flares emitted by the host star (Wheatley et al. 2017; Bourrier et al. 2017; Vida et al. 2017). Optically thick and gray clouds were also included in the model, such that the pressure parameter corresponds to the pressure either at the ground or at the top of the cloud. We parameterized the forward model to sample the pressure and H₂/He abundance relative to all gases, using the two-component centered-log-ratio transform (Benneke & Seager 2013), also illustrated by the horizontal axis stretch of Figure 3.

Zoom In Zoom Out Reset image size

For the simultaneous fit of the stellar contamination and planetary atmosphere, we used POSEIDON (MacDonald & Madhusudhan 2017; MacDonald 2023) to analyze the transmission spectrum of each visit individually. We parameterized stellar contamination via a two-heterogeneity model defined by eight free parameters (Fournier-Tondreau et al. 2023; priors listed in Table 1): the spot and facula coverage fractions, temperatures, and local surface gravities, as well as the stellar photosphere temperature and surface gravity. We calculated the impact of stellar contamination by interpolating PHOENIX models (Husser et al. 2013) using the MSG package (Townsend & Lopez 2023). We fitted for the planetary atmosphere via a six-parameter model: the atmosphere radius at the 1 bar reference pressure, R_p,ref (with prior ${ \mathcal U }(0.85,1.15){R}_{{\rm{p}},\mathrm{obs}}$, where R_p,obs = 1.116 R_⊕; Agol et al. 2021), the atmospheric temperature, T (${ \mathcal U }(200,900)$ K), and the H₂, H₂O, and CO₂ mixing ratios, X_i (centered-log-ratio priors from 10⁻¹² to 1 with the remainder of the atmosphere filled with N₂; Lustig-Yaeger et al. 2023). We calculated each model transmission spectrum on a wavelength grid at a resolving power of R = 20,000 from 0.58 to 3.0 μm, including both the planetary atmosphere and stellar-contamination contribution.

The transmission spectra of TRAPPIST-1 b notably differ between the two visits (Figures 2(a)–(b)), with different wavelength-dependent slopes, spanning a few hundreds of parts per million, and different absorption band-like features. These differences can be explained by stellar contamination and its variation between the two epochs.

The stellar-contamination fit of the sequential fit shows that the transmission spectrum of the first visit is explained by an unocculted spot (with no evidence for faculae), with a covering fraction of ${29}_{-10}^{+11}\, \%$ and a temperature difference with respect to the photosphere of $-{197}_{-62}^{+56}$ K (Figure 2(a) and Table 1). The second visit shows a preference for an unocculted facula (with no evidence for spots), causing the strong slope in the blue part of the spectrum (Figure 2(b)). The covering fraction of this facula is ${29}_{-13}^{+12}\, \%$ and its temperature difference is $+{153}_{-33}^{+61}\,$K. The stellar-contamination model cannot perfectly fit the observed spectra, which is explained by the mismatch between the stellar models and the observations (see Appendix D), but it can, to some extent, explain the slope blueward of 1.3 μm as well as some features near 1.4 μm, especially for the second visit.

The results from the simultaneous stellar-contamination and planetary atmosphere fit (Figure 4 and Appendix E) are consistent with the unocculted spot and faculae properties inferred from the stellar-contamination fit of the sequential fit (Table 1). The first visit is explained by unocculted spots about 200 K cooler than the stellar photosphere, while the second visit requires unocculted faculae about 160 K hotter than the photosphere (Figures 4(c)–(d)). The stellar-contamination model is moderately to strongly preferred over a flat spectrum model, with log Bayes factors of 5.3 and 2.9 for the first and second visits, respectively, following Table 1 from Trotta (2008). Comparing the transmission spectra from the two visits to a stellar-contamination-only model and a joint stellar-contamination and planet atmosphere model (Figures 4(a)–(b)) shows that the observations can be described by stellar contamination alone; there is no significant improvement from adding an atmosphere. The χ² from the fits are as follows: for visit 1, χ² = 44.0 (with 41 degrees of freedom) for a retrieval with only stellar contamination and 43.5 (with 37 degrees of freedom) for stellar contamination and a planetary atmosphere. For visit 2, χ² = 30.6 and 29.6 for the same models.

Zoom In Zoom Out Reset image size

The inferred stellar heterogeneity parameters and atmospheric constraints from the simultaneous fit are generally similar to those inferred from the sequential fit (Section 6.2), likely due to the lack of features from the planetary atmosphere. These two analyses of the NIRISS/SOSS observations strongly support previous evidence that stellar contamination is a critical consideration when interpreting any transmission spectrum of planets transiting active stars (e.g., Moran et al. 2023), in particular TRAPPIST-1 planets (Zhang et al. 2018; Wakeford et al. 2019; Ducrot et al. 2020; Garcia et al. 2022). It is intriguing though that the amplitude of the TLS effect is over an order of magnitude larger in the two TRAPPIST-1 b spectra (∼750–1000 ppm peak to peak) than in the TRAPPIST-1 g spectra obtained with NIRSpec/Prism over two visits (B. Benneke et al. 2023, in preparation), also part of the JWST TRAPPIST-1 atmospheric reconnaissance program. This may indicate that the amplitude of the TLS effect varies strongly with epoch and/or transit impact parameter (b ∼ 0.1 and b ∼ 0.4 for TRAPPIST-1 b and g, respectively; Agol et al. 2021), but this remains speculative given the small number of observations. All this underscores the need for more observations to quantify better the impact of stellar contamination at any given epoch and geometrical configuration.

The median error bar of the stellar-contamination-corrected and combined transmission spectrum at R ∼ 15 is 89 ppm, not accounting for the ∼1800 ppm per-pixel model-driven uncertainty associated with the lack of stellar model fidelity (see Appendix D). Signatures of a cloud-free atmosphere on TRAPPIST-1 b could have amplitudes of 100–200 ppm, depending on the exact atmospheric conditions (Lincowski et al. 2018). The precision achieved with two NIRISS/SOSS transits is thus insufficient to detect or reject such atmospheres confidently yet, but it shows that this objective could be within reach with a reasonable number of visits considering the expected lifetime of JWST and its instrument.

The frequentist comparison between atmospheric forward models and the stellar-contamination-corrected and combined transmission spectrum rejects cloud-free, hydrogen-rich atmospheres at 1× and 100× solar metallicity at 29 and 16σ, respectively, further confirming previous analyses of HST observations (de Wit et al. 2016). We also compared the transmission spectrum to a limited set of cloud-free secondary atmospheres, namely pure CH₄, H₂O, and CO₂ models, as well as a clear Titan-like scenario (Figure 2(d)), but none of these scenarios can be confidently confirmed or rejected at present (0.82, 0.48, 0.84, and 0.44σ, respectively). We conclude that the current NIRISS transmission spectrum (two visits), after stellar-contamination correction, does not show any signatures of an atmosphere.

The sequential fit of stellar contamination and planetary atmosphere finds that the TLS-corrected and combined transmission spectrum is only consistent at 3σ with a particular range of atmospheres (Figure 3), from clear atmospheres with a mean molecular mass greater than approximately 9 amu to low-mean-molecular-mass atmospheres with cloud surface pressures less than 0.1 bar. However, if TRAPPIST-1 b has retained a primordial atmosphere, cloud formation could be unlikely given the amount of radiation the planet receives (Morley et al. 2015; de Wit et al. 2016). Clouds are also unlikely, at least on the dayside of TRAPPIST-1 b, given the near-zero Bond albedo inferred from the JWST/MIRI secondary eclipse observations (Greene et al. 2023). These Bayesian results are consistent with the frequentist rejection of the clear hydrogen-rich forward models.

The simultaneous fit of stellar contamination and planetary atmosphere finds a 3σ upper limit on the H₂ abundance of 52% from the first visit (Figure 4(c) and Appendix E), consistent with the results from the sequential fit. The second visit alone is inconclusive on the H₂ abundance, as the faculae-dominated spectrum yields a bimodal posterior distribution, allowing an unphysical atmosphere with 100% H₂ and no trace volatiles (Figure 4(d) and Appendix E). The peaks in the posteriors for H₂O and CO₂ are manifestations of the uninformative centered-log-ratio prior that treats all the gases equally (i.e., if there is an atmosphere with a high mean molecular mass, either or both H₂O and CO₂ could be the background gas). The results from this simultaneous retrieval of the planetary atmosphere and stellar contamination remain consistent with a clear, high-mean-molecular-mass atmosphere dominated by N₂, H₂O, and/or CO₂ (see the H₂O and CO₂ histogram peaks in Figures 4(c)–(d)).

In addition to being consistent with pre-JWST observations, the conclusions from both the sequential and simultaneous fits of stellar contamination and planetary atmosphere are consistent with the conclusions from the JWST/MIRI secondary eclipse observations (Greene et al. 2023; Ih et al. 2023) in that the transits reveal no evidence for an atmosphere. Both the transit and secondary eclipse observations agree that TRAPPIST-1 b does not have a thick atmosphere, at least not one with a low mean molecular mass based on the transit observations. The MIRI observations provide additional constraints by rejecting certain atmospheres containing CO₂, as well as atmospheres containing no CO₂ for certain pressures and compositions (right panel of Figure 2 from Ih et al. 2023), further ruling out regions of the parameter space allowed by the transmission data. In this sense, the transmission spectra are less constraining in terms of atmospheric characterization for TRAPPIST-1 b, but they partially confirm the emission data and, more importantly, they bring insight into the stellar-contamination aspect of transmission spectroscopy, which is currently the preferred technique to probe the potential atmospheres of the outer, cooler TRAPPIST-1 planets.

Being closest to the host star, TRAPPIST-1 b is the planet in this system most likely to have lost its atmosphere through atmospheric escape processes. Because of the highly irradiated environment that TRAPPIST-1 b is subjected to, the potential absence of an atmosphere on TRAPPIST-1 b would not necessarily forbid atmospheres on the other planets in the system, especially the outer planets (Krissansen-Totton 2023). Moreover, the nondetection of an atmosphere on TRAPPIST-1 b can provide insight into the evolutionary history of the system. Results from this work and from the eclipse observations (Greene et al. 2023; Ih et al. 2023) should thus not discourage anyone from monitoring transits of the TRAPPIST-1 planets to search for atmospheric signatures.

The supreme-SPOON pipeline has already been presented in detail in Radica et al. (2023; see also Feinstein et al. 2022; Coulombe et al. 2023), and the steps we follow here closely mirror those other analyses. We process the data of both visits through the standard supreme-SPOON stage 1, which includes a group-level correction of 1/f noise, and stage 2. During background correction, we found that the SOSS SUBSTRIP256 background model provided by STScI could not completely remove the pick-off mirror jump in the zodiacal background near spectral pixel 700 in Figure 5(b). We therefore constructed a new background model from a high signal-to-noise ratio (S/N) stacked image of each TRAPPIST-1 time series. We first masked all traces (target and field contaminants) and bad pixels. Then, assuming each column sees a constant background level, we adopted the 10th percentile value of each column and low pass filtered these values to produce a smoothed background. We performed these steps to build the background model with the deep stack rotated such that the background discontinuity at x ≈ 700 was perfectly vertical. The model was then reprojected in 2D, derotated, and subtracted from the time series. The 1D spectral extraction was performed with a 25 pixel wide box aperture. Dilution introduced due to the overlap of the first and second spectral orders on the detector is predicted to be negligible for relative flux measurements (e.g., transmission spectra; Darveau-Bernier et al. 2022; Radica et al. 2022); we estimate that the amount of dilution introduced by the order overlap would be less than 5 ppm throughout the 0.6–2.8 μm wavelength range (following Equation (5) from Darveau-Bernier et al. 2022), significantly smaller than the transit depth precision.

Zoom In Zoom Out Reset image size

We also reduced the data with the NAMELESS pipeline (Feinstein et al. 2022; Coulombe et al. 2023; Radica et al. 2023; Figure 5), which first runs all the steps of the jwst pipeline stage 1 except for dark current subtraction, as it leaves residual 1/f noise. A parallel reduction was also performed with the dark current subtraction, but it had virtually no effect on the transmission spectrum. We then went through the assign_wcs, srctype, and flat_field steps of the jwst pipeline stage 2 before performing custom routines for the correction of bad pixels, background, as well as 1/f noise (Coulombe et al. 2023).

As with the supreme-SPOON pipeline, we found that the 2D background model from STScI could not reproduce the observed background jump near spectral pixel 700. To address this issue and to avoid introducing dilution through background correction, we split the background model into two separate sections along the line tracing the jump in background flux and fitted their amplitude separately. The background model made from the two distinct sections was then subtracted from all integrations. Spectral extraction of the first and second spectral orders was performed with a 24 pixel wide box aperture.

As a third reduction pipeline, we used the SOSS-Inspired SpectroScopic Extraction (SOSSISSE) tool. The SOSSISSE tool uses rateints files, i.e., the post–ramp-fitting time series of 2D images of the trace, to express any given trace image in the time series as the linear combination of a reference trace and its spatial derivatives. This approach assumes that the trace will only suffer achromatic transforms and provides a residual map from which spectroscopic information is then retrieved. Small signals from planetary transits are thus encoded into residual maps that have a much smaller dynamic range than the input data set, such that low-amplitude effects (e.g., 1/f noise) are more readily handled.

The SOSSISSE tool first constructs a model that consists of (1) a high-S/N, normalized, out-of-transit, median trace, M_i,j , where i and j are the indices of the pixels in the x- (dispersion) and y- (cross-dispersion) directions, respectively, and (2) time-dependent morphological changes in the trace. These morphological changes are represented by spatial derivatives of M_i,j : shifts in the x- and y-directions (∂M/∂x and ∂M/∂y), a trace rotation (∂M/∂θ), and changes in the contrast of the point-spread function (PSF) structures expressed as the second derivative along the y-direction (∂² M/∂y²). The model can be written as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{i,j}(t) & = & a(t)\,{M}_{i,j}+b(t)\,\displaystyle \frac{\partial {M}_{i,j}}{\partial x}+c(t)\,\displaystyle \frac{\partial {M}_{i,j}}{\partial y}\\ & & +d(t)\,\displaystyle \frac{\partial {M}_{i,j}}{\partial \theta }+e(t)\,\displaystyle \frac{{\partial }^{2}{M}_{i,j}}{\partial {y}^{2}},\end{array}\end{eqnarray} \tag{ A1 }$$

where F_i,j (t) is the model trace flux at time t on pixel (i, j), and a, b, c, d, and e are time-dependent scalars (Figures 1(e)–(h) and (m)–(p)), fitted in a least-squares fashion to each 2D trace image. Scalars b, c, d, and e can be interpreted as the amplitudes of their respective derivatives (e.g., b(t) can be seen as "dx," the amplitude of the trace displacement in the x-direction).

These morphological changes in the trace should be accounted for to measure the photometric signal accurately, and the SOSSISSE tool includes them as part of the photometric measurement instead of attempting to remove them through detrending in a subsequent step. Indeed, the fitted amplitude a(t) of the high-S/N trace consistently accounts for these morphological changes in the trace, ensuring that the photometric measurements are consistent across the entire time series. We also found that e(t) is particularly effective at detecting telescope tilt events (Rigby et al. 2023; Doyon et al. 2023). These events can affect the shape of the trace, but may only slightly shift its overall position. Including this term in the model accurately captures the impact of any potential tilt events on the data.

Once an optimal trace model has been found for each 2D image of the time series, SOSSISSE subtracts this best-fit trace model from each image, leaving only chromatic temporal variations, or "residuals." SOSSISSE then collapses each 2D residual image into a 1D spectrum by performing a least-squares adjustment of the model trace F at spectral pixel (i, j) onto the residual. The underlying assumption is that residuals resulting from spectroscopic features will have the same morphology as the underlying trace. This least-squares fitting is performed with the propagation of the uncertainty from the rateints files. We also account for the probability that the per-pixel residuals are either consistent with their uncertainties or with a cosmic-ray event by giving to each pixel a weight corresponding to 1 minus its likelihood of being affected by a cosmic ray. The time series a(t), which can be interpreted as the broadband light curve, is finally added back to each spectral channel, yielding a spectral time series that can be used to perform the transit light-curve fit (Section 4).

Similar structures are present across the spectral time series of the three pipelines. For example, the temporal slope in visit 1 and the same correlated noise pattern associated with stellar activity near the transit egress in visit 2 are present in all three reductions. There are some differences between reductions, typically localized in time and wavelength, but these are expected given the differences between the pipelines.

We fitted the broadband and spectroscopic light curves of all three reductions to compare the resulting transmission spectra (Figure 6(a)). Although the transit depths are generally consistent within their uncertainties, there are differences at the 100–200 ppm level in some wavelength bins. Differences of this order of magnitude between reductions have already been reported in JWST data (Feinstein et al. 2022; Coulombe et al. 2023; Rustamkulov et al. 2022; Radica et al. 2023), but these become more problematic for Earth-sized planets like TRAPPIST-1 b, whose expected atmospheric signatures have amplitudes of 100–200 ppm (Lincowski et al. 2018) at most.

Zoom In Zoom Out Reset image size

To assess the impact of the differences between reductions, we ran the simultaneous fit of the stellar contamination and planetary atmosphere (Section 5.2) on all three sets of TRAPPIST-1 b transmission spectra resulting from the three pipelines, generally finding consistent conclusions. Regardless of the reduction, we see evidence of stellar contamination and cloud-free, hydrogen- and helium-rich models can be ruled out at high confidence. Secondary atmospheres, especially those with a low surface pressure, and/or clouds cannot be confidently rejected from these SOSS observations. Since the conclusions remain the same across all three reductions, we consider that our reductions are reliable for the purpose of this work.

We started the transit light-curve fit with a set of approximately 2048 light curves, with each light curve corresponding to a spectral channel. For spectral order 1, we produced the broadband light curve by binning all light curves into a single bin, whereas the spectroscopic light curves were generated by binning the light curves down to 50 bins, equally spaced in wavelength. For spectral order 2, we first removed all light curves at wavelengths above 0.82 μm, as this spectral domain is already covered by spectral order 1 (Albert et al. 2023) and removing these redder wavelengths also prevents contamination from order 1. Given the lower S/N and the fewer native-resolution channels in order 2, we only performed a broadband light-curve fit for this order. We then normalized each light curve by its median and performed 10 iterations of 3.5σ clipping to remove outliers in flux.

We first fitted the broadband light curves from the two visits independently with ExoTEP (Benneke et al. 2017, 2019; Benneke 2019; Figures 1(a)–(c) and (i)–(k)), which fits a batman transit model (Kreidberg 2015) and a systematics model to the light curve, along with a photometric noise parameter, using the emcee implementation (Foreman-Mackey et al. 2013) of the Markov Chain Monte Carlo method, with the priors listed in Table 2. We applied a Gaussian prior on the impact parameter based on the value and uncertainty from Agol et al. (2021) given the low temporal resolution during the transit ingress and egress. For the GP length scale (Appendix B.3), we applied a lower bound corresponding to about five integrations to prevent overfitting, and an upper bound of 2.4 hr, approximately half the duration of the observation, to prevent degeneracies with the parameters of the linear function (Appendix B.3). These independent fits of the two broadband light curves provided a first estimate of the global transit parameters (a/R_⋆, b), transit ephemeris (T_mid), and transit depths in the SOSS bandpass. We then fitted the two broadband light curves simultaneously (posteriors listed in Table 2), allowing each visit to have different systematics and photometric noise parameters, initializing all parameters at their best-fit values from the individual fit to ease convergence.

Table 2. Priors and Posteriors of the Joint Broadband Light-curve Fit of the Two Visits

White Light Curve
Parameter	Priors	Posteriors
Orbital period, Porb (days)	1.51087081	1.51087081
Midtransit time, visit 1, Tmid,1 (BJD)	${ \mathcal U }(\mathrm{2,459,779.1},\mathrm{2,459,779.3})$	${\mathrm{2,459,779.210475}}_{-0.000025}^{+0.000025}$
Midtransit time, visit 2, Tmid,2 (BJD)	Tmid,1 + 1.51087081	Tmid,1 + 1.51087081
Planet-to-star radius ratio, (Rp /Rs )	...	${0.08501}_{-0.00049}^{+0.00047}$
Planet-to-star radius ratio squared, ${({R}_{p}/{R}_{s})}^{2}$ (ppm)	${ \mathcal U }(100,\mathrm{50,000})$	${7226}_{-83}^{+80}$
Impact parameter, b	${ \mathcal N }(0.095,0.65)$	${0.098}_{-0.045}^{+0.046}$
Quadratic LD coeff. 1, q1	${ \mathcal U }(0,1)$	${0.409}_{-0.081}^{+0.101}$
Quadratic LD coeff. 2, q2	${ \mathcal U }(0,1)$	${0.266}_{-0.072}^{+0.080}$
Eccentricity, e	0	0
Argument of periapsis, ω (deg)	90	90
Semimajor axis, a/R⋆	${ \mathcal U }(10,30)$	${20.66}_{-0.15}^{+0.12}$
White noise, visit 1, s1 (ppm)	${ \mathcal U }({10}^{-10},{10}^{10})$	${172.3}_{-11.2}^{+12.3}$
White noise, visit 2, s2 (ppm)	${ \mathcal U }({10}^{-10},{10}^{10})$	${164.10}_{-9.89}^{+10.86}$
GP amplitude, visit 1, ${\mathrm{log}}_{10}({a}_{\mathrm{GP},\ 1})$	${ \mathcal U }(-10,1)$	$-{3.79}_{-0.18}^{+0.26}$
GP amplitude, visit 2, ${\mathrm{log}}_{10}({a}_{\mathrm{GP},\ 2})$	${ \mathcal U }(-10,1)$	$-{3.68}_{-0.11}^{+0.14}$
GP length scale, visit 1, ${\mathrm{log}}_{10}({\lambda }_{\mathrm{GP},\ 1}(\mathrm{days}))$	${ \mathcal U }(-2.21,-1)$	$-{1.62}_{-0.31}^{+0.18}$
GP length scale, visit 2, ${\mathrm{log}}_{10}({\lambda }_{\mathrm{GP},\ 2}(\mathrm{days}))$	${ \mathcal U }(-2.21,-1)$	$-{2.01}_{-0.12}^{+0.28}$

Spectroscopic Light Curves
Parameters	Priors
Orbital period, Porb (days)	1.51087081
Midtransit time, Tmid (BJD)	WLC best fit
Planet-to-star radius ratio squared, ${({R}_{p}/{R}_{s})}^{2}$ (ppm)	${ \mathcal U }(100,\mathrm{50,000})$
Impact parameter, b	WLC best fit
Quadratic LD coeff. 1, q1	${ \mathcal U }(0,1)$
Quadratic LD coeff. 2, q2	${ \mathcal U }(0,1)$
Eccentricity, e	0
Argument of periapsis, ω (deg)	90
Semimajor axis, a/R⋆	WLC best fit
White noise (ppm)	${ \mathcal U }({10}^{-10},{10}^{10})$
GP amplitude, ${\mathrm{log}}_{10}({a}_{\mathrm{GP}})$	${ \mathcal U }(-10,1)$
GP length scale, ${\mathrm{log}}_{10}({\lambda }_{\mathrm{GP}}(\mathrm{days}))$	${ \mathcal U }(-2.21,-1)$

Note. For the posteriors, the leftmost value is the 50th percentile of the posterior distribution, and upper and lower uncertainties are computed from the 16th, 50th, and 84th percentiles. Posteriors are not listed for the spectroscopic light curves as they are different for each spectral bin.

We fitted the spectroscopic light curves of the two transits independently, with a transit and systematics model as well as a photometric noise parameter for each wavelength bin (Figure 7), with the priors listed in Table 2. We parameterized the quadratic limb darkening parameters following Kipping (2013). The transmission spectrum for each visit (Figures 2(a)–(b)) is given by the best-fit spectroscopic R_p/R_⋆ values squared.

Zoom In Zoom Out Reset image size

The broadband light curve of the second visit (Figure 1(i)) betrays the presence of 3–15 minute long correlated noise structures produced by stellar activity (Section 4). This correlated noise is most easily seen as flux variations before and after the transit egress of the second visit, but more subtle modulations are also present elsewhere in both visits. To mitigate the impact of this correlated noise on the measured transit depths, we tested two different systematics treatments for both transits, and three additional treatments for the second visit, specifically to leverage the apparent correlation between the Hα integrated flux and the flux in other wavelength bins. We tested each of these five systematics treatments in parallel during the transit light-curve fitting (Appendix B.2). In all five systematics treatments, we used a linear function as the basic systematics model S(t), such that the most basic flux model is written as

$$\begin{eqnarray}&&\begin{array}{l}f(t)={f}_{0}\cdot T(t,\theta )\cdot S(t),\quad \mathrm{with}\\ S(t)=v\cdot (t-{t}_{0})+c,\end{array}\end{eqnarray} \tag{ B1 }$$

where f₀ is the baseline flux, T(t, θ) is the transit model, v and c are the slope and normalization constant of the linear function, respectively, and t₀ is the time at the start of the visit.

The first systematics treatment consisted of trimming the baseline and masking temporally correlated points such that the most obviously correlated noise patterns would not affect the fit. For visit 1, we removed 2 hr and 0.6 hr from the start and the end of the visit, respectively. For visit 2, we masked six in-transit and 12 posttransit points, all of which seemed correlated in time and with Hα. We also removed 2 hr from the start, but did not remove any additional parts from the end, since 12 points were already masked, such that trimming would have left very little posttransit baseline.

The second systematics treatment consisted in using a GP (Rasmussen & Williams 2006). We used a squared-exponential kernel to define the covariance matrix of the likelihood function using the george Python package (Ambikasaran et al. 2015), keeping the white photometric noise parameter on the diagonal of the covariance matrix, added in quadrature to the kernel function. This added two free hyperparameters to the fit of each light curve: the amplitude a_GP and timescale λ_GP of the GP.

The subsequent systematics treatments were only applied to the second visit as they all involved the Hα luminosity, which was most obviously correlated with the flux at other wavelengths for this visit. The third systematics treatment consisted of detrending the light curves against the Hα integrated flux, that is

$$\begin{eqnarray}&&S(t)=v\cdot (t-{t}_{0})+c+w\cdot ({F}_{{\rm{H}}\alpha }(t)-{F}_{{\rm{H}}\alpha ,0}),\end{eqnarray} \tag{ B2 }$$

where w is a free parameter for each light curve, F_Hα (t) is the Hα luminosity time series, and F_Hα,0 is the Hα luminosity at the start of the visit.

In the fourth systematics treatment, we first trained a squared-exponential kernel GP on the Hα luminosity time series, yielding a Gaussian-like posterior distribution function (PDF) for both a_GP and λ_GP. We used the mean and standard deviation of the PDF of λ_GP to define a Gaussian prior applied to λ_GP as we fitted the broadband and spectroscopic light curves with a squared-exponential kernel GP.

The fifth systematics treatment was a combination of the second and third treatments, that is, detrending against the Hα luminosity and using a (untrained) squared-exponential kernel GP.

Transit depths are consistent within their uncertainties across all systematics treatments (Figure 6(b)). This shows that the GP does not bias the measured R_p/R_s, at least not in a way that is different from simply masking parts of the light curves or simply detrending against the Hα luminosity. The median error bars, in parentheses in the legend, are larger for the trimmed baseline treatment and any treatment involving a GP. The larger error bars for the trimmed baseline treatment is simply explained by the reduced amount of data, especially in visit 2 where in-transit points were masked. We explain the larger error bars in the GP treatments by the fact that the GP hyperparameters may be correlated with the transit depths, which inflates their error bars. However, a smaller median error bar on the transit depths does not necessarily mean a better fit to the light curves. A visual inspection of the light-curve residuals from the best-fit model of each wavelength bin for each systematics treatment indicated that the uninformed GP (treatment 2) provided a better fit, that is, it yielded residuals with the least correlated noise. This is expected given its more flexible nature compared to detrending against the fixed Hα time series. This is why we selected treatment two as our reference.