High Tide or Riptide on the Cosmic Shoreline? A Water-rich Atmosphere or Stellar Contamination for the Warm Super-Earth GJ 486b from JWST Observations

We use a modified version of the jwst Stage 1 pipeline, starting from the _uncal.fits files. We perform group-level background subtraction before determining the flux per integration. For each group, we exclude the region within 9 pixels of the trace before computing and subtracting a median background value per pixel column. We process the _rateints.fits files through the regular jwst Stage 2 pipeline, skipping the flat-fielding and absolute photometric calibration steps when our goal is to derive the planet's spectrum at later stages. Conversely, we include these steps when our goal is to compute the flux-calibrated stellar spectrum (see Section 4.4). Stage 3 of Eureka! converts the time series of 2D integrations into 1D spectra using optimal spectral extraction (Horne 1986) and an aperture within 5 pixels of the trace. We flag bad pixels at numerous points within this stage using thresholds optimized to minimize scatter in the white light curves.

For the NRS1 detector, we extract the flux from 2.777 to 3.717 μm and split the light into 47 spectroscopic light curves, each 20 nm (0.02 μm) in width. For the NRS2 detector, we adopt the same resolution in extracting 67 spectroscopic light curves spanning 3.825–5.165 μm. For each detector, we manually mask 9 pixel columns that exhibit significant scatter in their individual light curves. Doing so improves the quality of the spectroscopic light curves and yields more consistent transit depths.

With two NIRSpec detectors and two transit observations, we fit four white light curves and their systematics (see Figure 1). We determine the system parameters using batman (Kreidberg 2015) and fix the quadratic limb-darkening coefficients to those provided by ExoTiC-LD (Grant & Wakeford 2022), assuming the stellar parameters given by Trifonov et al. (2021) and the MPS-ATLAS set 1 models (Kostogryz et al. 2023). For the NRS1 detector, we find that a quadratic trend in time provides the best fit. For the NRS2 detector, a linear trend suffices to remove systematics. Table 1 in Appendix A lists our best-fit system parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images

When fitting the spectroscopic light curves (see Figure 1), we fix the planet's transit midpoint, inclination, and semimajor axis to the weighted mean values in Table 1. We fix the quadratic limb-darkening parameters to the values provided by ExoTiC-LD for each spectroscopic channel. For the NRS1 detector, we also fix the quadratic term in our time-dependent systematic model to that of the best-fit white light-curve value (Transit 1: c₂ = 0.0335, Transit 2: c₂ = 0.0248). For all spectroscopic light curves, we fit for the zeroth and first-order terms (c₀ and c₁) of our polynomial. Light curves from the NRS2 detector only require a linear model in time. Including the term that rescales the uncertainties, each spectroscopic light curve has four free parameters, of which only the planet-to-star radius ratio is a physical parameter.

For each light curve, we first perform a least-squares minimization using the Powell method (Powell 1964) and then initialize our Markov chain Monte Carlo routine using our best-fit values. We estimate the parameter uncertainties using emcee (Foreman-Mackey et al. 2013) and, at each iteration, we increase the uncertainties by an average factor of ∼1.5 to achieve a reduced χ² = 1. All of our posteriors are Gaussian distributed, and there are no parameter degeneracies.

We run the jwst pipeline through Stages 1 and 2 using the uncal.fits files. We utilize group-level 1/f subtraction and apply a scaled superbias to account for the vertical offset seen in NRS2 Transit 1 (see Appendix A.1). We correct for cosmic rays and bad and hot pixels in the Stage 2 output rateints.fits files and apply a second 1/f correction at the integration level by masking the spectral trace and then calculating the median of the background pixels in each column. This value is then subtracted from the cleaned 2D image.

We next cross-correlate each 2D image with the median-aligned image to determine the x- and y-shifts of the spectral trace, which are used to align all 2D images. A Gaussian profile is then cross-correlated to each column in the y-direction, and a fourth-order polynomial is fit in the x-direction to determine the spectral trace, which is used to extract the spectra.

The white light curves for Transits 1 and 2 are fit from the extracted spectra by summing the spectra in the wavelength direction over a detector. We fit a/R_⋆, limb-darkening parameters, and the impact parameter b using the weighted mean from both transits and both detectors. We then fix a/R_⋆, b, and the period, and fit for R_P /R_⋆, T₀, and limb darkening in the white light curve. A low-order polynomial in time (third-order in NRS1 and up to fourth-order for NRS2) was used to model the baseline, with additional detrending parameters of the x- and y-shifts and superbias scale factor. We then fix the system parameters (presented in Table 2 in Appendix A) and limb-darkening coefficients in each wavelength column to fit the spectroscopic light curves.

With Tiberius we started by running STScI's jwst stage 0 pipeline on the uncal.fits files from the group_scale step through gain_scale step. We set --odd_even_columns = True at the ref_pix step and ran our own 1/f correction step at the group level prior to running ramp_fit, which removes the median background flux for every column of every group's spectral image. We define the background as a 14 pixel wide region that avoided 18 pixels centered on the curved trace and mask bad pixels using our own custom bad-pixel map. We subsequently ran assign_wcs and extract_2d to obtain the wavelength solution and proceeded to run Tiberius's spectral extraction on the gainscalestep.fits files.

First we oversampled each pixel by a factor of 10 using a linear interpolation. This allows us to measure the stellar flux at the subpixel level, which reduces noise in the light curves (The JWST Transiting Exoplanet Community Early Release Science Team et al. 2022). We used a fourth-order polynomial to trace the NRS1 detector stellar spectrum and a sixth-order polynomial for NRS2. We performed standard aperture photometry at every pixel column, with a 4 pixel wide aperture. We performed an additional background subtraction step at this stage by calculating the background in 14 pixels on either side of the trace, excluding 7 pixels on each side. For NRS1 we fit these background pixels with a linear polynomial, while for NRS2 we used a median since our defined background regions were mostly above the stellar trace.

We remove cosmic rays and residual bad pixels manually and then correct for small shifts in the stellar spectra along the dispersion direction by cross-correlating all spectra in the time series with the first, resampling each spectrum onto a common pixel grid. Finally, we created a white light curve between 2.75 and 3.72 μm for NRS1 and 3.83–5.15 μm for NRS2. Our spectroscopic light curves were created at 1 pixel resolution over the same wavelength range.

We fit the four white light curves (two transits ×two detectors) with batman (Kreidberg et al. 2015), leaving a/R_*, R_P /R_*, the orbital inclination (i), and the time of mid-transit (T₀) as free parameters, and fixing the period to the value from Trifonov et al. (2021). For our white and spectroscopic light curves, we assumed quadratic limb darkening with coefficients fixed to values from 3D stellar atmosphere models (Magic et al. 2015) using ExoTiC-LD (Grant & Wakeford 2022). We adopted T_eff = 3340 K, [Fe/H] = 0.070, and $\mathrm{log}{g}_{* }=4.9155$ (Trifonov et al. 2021). For our systematics model we used a combination of polynomials: quadratic-in-time, linear-in-x-position, and linear-in-y-position, resulting in nine free parameters: four transit model parameters and four systematics model parameters.

To determine the best-fitting values and uncertainties, we used emcee (Foreman-Mackey et al. 2013) with 90 walkers for two runs of 20,000 steps. After the first run we inflated our photometric uncertainties to give a reduced χ² = 1 for our best-fitting model before the second run. Table 3 in Appendix A summarizes the results of our white light-curve fits. For our spectroscopic light-curve fits, we fixed a/R_*, i, and T₀ to the weighted mean values from our four white light-curve fits and only fitted for R_P /R_* and the five parameters defining our systematics model. Here we used a Levenberg–Marquadt sampler for computational speed as we had to fit 6876 spectroscopic light curves.

We performed a flat line hypothesis rejection test to determine the statistical significance of the slope in the transmission spectrum. We fitted the spectrum from each pipeline using two models: a flat, featureless model that uses one free parameter for the transit depth and a Gaussian spectral feature model with four free parameters: the flat transit depth and the central wavelength, amplitude, and width of a Gaussian feature added to the baseline featureless spectrum. We fitted both models to each data set using the dynesty nested sampling code for Bayesian inference (Speagle 2020) and then used the Bayesian evidence to calculate the Bayes factor of each model (e.g., Trotta 2008, 2017). We then converted the Bayes factors to more classical "sigma" detection significances using the relationship detailed by Benneke & Seager (2013).

Figure 2 demonstrates that each spectrum separately favors the Gaussian model and rejects a featureless spectrum. The strength of the signal detection is 3.20σ for Eureka!, 2.24σ for FIREFLy, and 3.29σ for Tiberius. The FIREFLy detection significance is lower due to slightly larger uncertainties associated with that reduction, which stem from FIREFLy's choice of spectroscopic binning to produce similar transit depth errors across the full wavelength range and wavelength-dependent baseline functions. Nevertheless, the same shape is seen in the spectra from the three pipelines. Thus, the flat line hypothesis is rejected by all three analyses with varying confidence. Each individual reduction hypothesis rejection test is available.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (4 panels)

Tarball of figure set images

We ran a suite of forward models using the stellar and planet parameters from Caballero et al. (2022) to compare to each transmission spectrum. We also generated forward models using an updated stellar log(g) = 4.91 ± 0.02, obtained from our updated a/R_s constraints (see Appendix A) (Seager & Mallén-Ornelas 2003; Sandford & Kipping 2017), finding consistent results.

We focus on higher mean molecular weight scenarios to explain the transmission spectrum. For completeness, however, we simulate a 1000× solar metallicity atmosphere with a parameterized pressure-temperature profile in thermochemical equilibrium with CHIMERA (Line & Yung 2013; Line et al. 2014) as in our previous work (Lustig-Yaeger et al. 2023). We include the species H₂O, CH₄, CO, CO₂, NH₃, HCN, H₂S, H₂, and He. The CHIMERA thermochemical equilibrium abundances result in a model spectrum that is primarily shaped by methane, carbon dioxide, and water. After generating the temperature-pressure profile and atmospheric abundances with CHIMERA, we use the radiative transfer suite of PICASO (Batalha et al. 2019), with opacities resampled to R = 10,000 from Batalha et al. (2020), to generate model spectra.

In each case, we bin the resulting model transmission spectrum to the resolution of the data before performing a reduced-χ² comparison.  . As with the Gaussian hypothesis tests, we exclude the data points in the gray shaded region of Figure 2 from our model fitting due to steeply falling instrument throughout at these wavelengths (<2.87).

As shown in Figure 3, the slight slope and flatness of the spectra from each reduction allow us to confidently disregard low mean molecular weight atmospheres dominated by hydrogen/helium—up to metallicities of 1000× solar—to greater than 3σ. This improves upon the previous high-resolution data obtained by Ridden-Harper et al. (2022) that could only strongly rule out atmospheres up to a few times solar. Our 1000× solar metallicity atmosphere has an average mean molecular weight of 13.86 g mol⁻¹ compared to the high resolution's 5 g mol⁻¹ limit though our constraint is less stringent for non-chemically consistent atmospheres (see Section 4.3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images

We also compare the data from each reduction to a set of end-member forward models from PICASO with single-gas 1 bar, isothermal atmospheres. For ease of interpretation, we focus here on the results from the Eureka! reduction, as we determined that it was the most representative data set, with the smallest weighted average deviation from the median of all three reductions. However, the trend in best fit agrees among all three reductions (for a complete description of each reduction's fit, see Appendix C). The slight slope on the blue end of NRS1 results in best-fitting (reduced χ² = 1.01) forward models that contain pure water vapor, as this molecule has a strong absorption feature from 2.2 to 3.7 μm, consistent with the slope we observe in NRS1.

Our data across all reductions also moderately to weakly rule out carbon-rich atmospheres of either CH₄ or CO₂ to 6.5σ and 2.3σ, respectively. A flat line model, representative of an airless body or a high-altitude (0.1 μbar) cloud deck, fit the data with reduced χ² = 1.11, which is statistically equivalent to the clear water atmosphere model within the forward modeling framework. However, between its equilibrium temperature and size, GJ 486b is not expected to support clouds to such low pressures, as there are few condensable species in this temperature range. Photochemical hazes could dampen the presence of any spectral features with a haze layer at this altitude and create a flat line spectrum (Gao et al. 2020; Pidhorodetska et al. 2021; Caballero et al. 2022); however, given the Bayesian evidence of the Gaussian absorption tests discussed above, the water atmosphere is the preferred explanation from the PICASO analysis for all reductions. We note that the FIREFLy reduction only weakly rejects the flat line hypothesis and, therefore, an airless planet or very hazy planet is still a possibility. In Figure 3, we show the results of our PICASO forward modeling compared to the Eureka! data. The full set of results for each reduction is available.

In addition to our forward model comparisons, we performed an atmospheric retrieval analysis to assess the robustness of our tentative evidence for a water-rich atmosphere and consider alternative astrophysical explanations. We apply two independent retrieval codes—POSEIDON (MacDonald & Madhusudhan 2017; MacDonald 2023) and rfast (Robinson & Salvador 2023)—to all three data reductions to ensure reliable inferences.

4.3.1. Water-rich Atmosphere Scenario

Our POSEIDON atmospheric retrieval considers six potential gases that can range in abundance from being trace volatiles to the dominant background gas: N₂, H₂, H₂O, CH₄, CO₂, and CO. The opacity contributions from these gases include line opacity (Tashkun & Perevalov 2011; Li et al. 2015; Yurchenko et al. 2017; Polyansky et al. 2018) and collision-induced absorption (CIA) from H₂–H₂, H₂–N₂, H₂–CH₄, H₂–CO₂, CO₂–CO₂, CO₂–CH₄, and N₂–N₂ (Karman et al. 2019). Since the mixing ratios must sum to unity, we have five free parameters describing their mixing ratios that each follow centered log-ratio (CLR) priors, ranging from 10⁻¹² to 1, as described by Lustig-Yaeger et al. (2023). The other free parameters are the isothermal temperature (${ \mathcal U }$ [200 K, 900 K]), the atmosphere radius at the 1 bar reference pressure (${ \mathcal U }$ [0.9 R_p,obs, 1.1 R_p,obs]), and the log-pressure of an opaque surface (${ \mathcal U }$ [−7, 2], in bar). We calculate transmission spectra via opacity sampling at a resolving power of R = 20,000 from 0.5 to 5.4 μm, with the lower wavelength limit set far below our shortest wavelength (2.8 μm) to later demonstrate how retrieval solutions diverge at optical wavelengths. These eight-parameter POSEIDON retrievals used the PyMultiNest (Feroz et al. 2009; Buchner et al. 2014) package to explore the parameter space with 2000 live points.

Figure 4 shows our POSEIDON retrieval results for this atmospheric model scenario (blue retrieved spectrum and histograms) for the Eureka! data reduction—see the figure set corresponding to Figure 4 for the other two reductions. For Eureka! and FIREFLy, the preferred explanation for the observed rise in the blue wavelengths of the transmission spectrum is H₂O opacity from the wing of the band centered on 2.8 μm. Bayesian model comparisons favor the presence of H₂O with Bayes factors of 133 and 8 (3.6σ and 2.6σ) for Eureka! and FIREFLy, respectively. The retrieved H₂O abundance posterior indicates that water is the most likely background gas (e.g., Eureka! requires an H₂O mixing ratio >10% to 2σ confidence), with an upper limit ruling out a H₂-dominated atmosphere. The Eureka! and Tiberius reductions also yield upper limits on the CH₄ and CO₂ abundances (see Appendix C). The Tiberius reduction, however, does not uniquely infer a water-rich atmosphere. Though a water-rich atmosphere remains the preferred solution for Tiberius, a secondary mode permits a clear, H₂-dominated atmosphere with no other gases contributing to the spectrum. This secondary mode reflects a solution where the wavelength dependence of H₂–H₂ CIA is used to fit the spectrum. This solution is unphysical since an H₂-dominated atmosphere will always contain other trace molecules with more prominent absorption features at these wavelengths. Upon further investigation, we found that the unphysical solution is driven by the upward rise at the longest wavelengths that are only present in the Tiberius reduction (see Figure 2). We, therefore, conclude that a consistent explanation for GJ 486b's transmission spectrum, assuming the observed non-flatness is caused by atmospheric absorption, can be readily explained (${\chi }_{\nu }^{2}\approx 1.0$) by a water-rich atmosphere—in agreement with the forward models in Section 4.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (3 panels)

Tarball of figure set images

We also conducted single-composition atmospheric retrievals with rfast for all three reductions. These retrievals consider atmospheres with a single absorbing gas alongside a spectrally inactive background gas with an agnostic mean molecular weight. Our rfast retrieval model has six free parameters: the log-gas mixing ratio, ${\mathrm{log}}_{10}{f}_{\mathrm{gas}}$ (${ \mathcal U }$ [−12, 0]); the log-surface pressure, ${\mathrm{log}}_{10}{P}_{0}$ (${ \mathcal U }$ [−1, 6], in Pa); the surface temperature, T₀ (${ \mathcal U }$ [300, 1100] K); the mean molecular weight of the background gas, m_b (${ \mathcal U }$ [2, 50] amu); the planet radius, R_p (${ \mathcal U }$ [1.1, 1.4] R_⊕); and the planet mass, M_p (${ \mathcal N }$ [2.28, 0.12] M_⊕). For the single gases, we consider, in separate retrievals, H₂O, CO₂, CO, and CH₄. The rfast retrievals use emcee (Foreman-Mackey et al. 2013) with 100 walkers for 15,000 steps, where the first 5000 are discarded for burn-in.

We show our rfast 1D posteriors in Appendix C. Our rfast retrievals also identify an H₂O-rich atmosphere as a consistent explanation for the Eureka! and FIREFLy reductions (though the lower limits on H₂O are weaker compared with POSEIDON due to the combination of a free mean molecular weight, planet mass, and log-uniform versus CLR priors). rfast also finds that the Tiberius reduction permits lower mean molecular weight atmospheres for similar reasons to POSEIDON.

4.3.2. Unocculted Starspot Scenario

To further investigate the possibility of stellar contamination, we return to the JWST/NIRSpec G395H data to probe the Stage 3 stellar spectra and examine whether the star is consistent with a particular stellar model. Upon completing Stage 2 of the jwst pipeline with the flat-fielding and absolute photometric calibration steps enabled, we noticed that only the region within 8 pixels of the trace is converted to units of MJy. The remaining pixel regions are in data numbers per second (DN/s), so we manually mask them before running Stage 3 of Eureka!. Due to the lack of unmasked background pixels, we disable Stage 3 background subtraction for this flux-calibrated reduction. This change does not skew the final calibrated spectrum since we already performed group-level background subtraction in Stage 1.

To compute the stellar baseline spectrum, we exclude 1040 integrations during transit (1560–2599) and then compute median values along the time axis. We manually mask a few obvious outliers before estimating the baseline spectrum uncertainties by computing the standard deviation in flux along the time axis. Typical uncertainties are 3–5 mJy but can be as large as 55 mJy for some spectral channels. The typical uncertainty values are consistent with the uncertainties derived from our standard spectral extraction routine. We do not use the standard error calculation for our uncertainties. That is, we do not divide our uncertainties by the square root of the number of integrations because, as demonstrated below, the standard deviation in flux better represents the true uncertainty in our flux-calibrated spectrum. We note that the derived baseline spectrum is remarkably consistent between both transits (see Figure 5).

Zoom In Zoom Out Reset image size

We used PHOENIX stellar models produced by Allard et al. (2012) to analyze whether the observed stellar baseline spectrum is best explained by a spotless or spotted star. We utilized the Allard et al. (2012) models, as in Section 4.3.2, because they account for the formation of molecular bands including H₂O, CH₄, and TiO₂ and have higher (Δλ = 2 Å) resolution than the observations. This grid of models also has sufficient temperature and gravity coverage to model the photospheres of M-dwarf stars and their spots and faculae (T_eff ≥ 2000 K, log(g) = 0–6 cm s⁻²).

We employed single PHOENIX models to represent spotless (or one-component) stars. We used weighted linear combinations of PHOENIX models to create inhomogeneous models. Two-component models include one model with T_eff ≥ 3000 K to represent the background photosphere and a second, cooler model with ${T}_{\mathrm{eff}}\leqslant {T}_{\mathrm{eff},\mathrm{photosphere}}-100$ K to represent spots. Three-component models include an additional ${T}_{\mathrm{eff}}\,\geqslant {T}_{\mathrm{eff},\mathrm{photosphere}}+100$ K model to represent faculae. In the two- and three-component models, all spots have the same T_eff and log(g), as do the faculae. Linear combinations were computed by interpolating the spot and faculae models onto the photosphere wavelength grid before summing the fluxes in a weighted fraction where the photosphere was required to be ≥50% of the total.

To compare the models to the observed baseline spectra, we converted the native wavelengths from ångstrom to micron and the flux densities from ergs s⁻¹ cm⁻² cm⁻¹ to fluxes in units of mJy. We then scaled the models by ${R}_{* }^{2}$/dist² using literature values for GJ 486: R_* = 0.33 R_⊙ (Trifonov et al. 2021) and dist = 8.07 pc (Gaia Collaboration et al. 2021). We smoothed and interpolated the models to be the same resolution as the observations before calculating a reduced χ². In our reduced-χ² calculations, we considered 3187 wavelength points for Transit 1 (3180 for Transit 2) and three fitted parameters (T_eff, log(g), and a scaling factor). The multicomponent models included additional fit parameters for determining the percent coverage for the spots and faculae. The scaling factor was multiplied by the ${R}_{* }^{2}$/dist² term to account for uncertainty in either measured quantity and varied from 0.9 to 1.1. To get the final reduced-χ² value for each model, we computed reduced χ² individually for Transits 1 and 2 and then took the average.

Considering each type of one-, two-, and three-component model individually, we find that the models with smallest reduced-χ² values are fairly consistent with the existing literature values though no model is a particularly good fit with a reduced χ² near 1 (for numerical details, see Table 4 in Appendix B). A 100% T_eff = 3300 K, log(g) = 4.5 cgs model with ${\chi }_{\nu }^{2}$ = 72.0 is the preferred one-component photosphere model (scale factor = 1.05), yielding a lower surface gravity than expected for a field-age mid-M dwarf like GJ 486. In agreement with our updated log(g) = 4.91 ± 0.02 cgs, we disfavor the low stellar surface gravity of the best-matched photosphere-only model when taking into account inhomogeneities on the stellar surface. A 75% T_eff = 3400 K, log(g) = 5 cgs background photosphere with 25% spot coverage at T_eff = 3000 K, log(g) = 5 cgs is the preferred two-component model (${\chi }_{\nu }^{2}$ = 53.4; scale factor = 1.05). The model most preferred overall is a three-component model with ${\chi }_{\nu }^{2}$ = 49.0 that has a background photosphere with T_eff = 3200 K, log(g) = 5 cgs, 20% spot coverage at T_eff = 3000 K, log(g) = 5 cgs, and 25% faculae coverage at T_eff = 3400 K, log(g) = 5 cgs (scale factor = 1.1). These three models are shown in Figure 5 compared to the baseline GJ 486 spectra from Transits 1 and 2. There is decent general agreement for each model throughout the full ∼2.9–5 μm range, with slightly better agreement for the three-component photosphere+spot+faculae model, indicating that we cannot rule out starspots as a source for the presumed water detection.

As stated in the main text, all initial reductions showed a consistent offset in measured transit depth for the Transit 1, NRS2 detector relative to the other white light-curve depths. Since this shift is not seen in the NRS1 detector, we can confidently rule out all astrophysical effects (e.g., stellar variability) as a source of the discrepancy. For the FIREFLy reduction, we altered our application of the superbias in the bias subtraction step and light-curve fitting stages, which we found produced more consistent transit depths for NRS1 and NRS2.

In our FIREFLy reduction, we measured the superbias level by rescaling the superbias image to match the level in the trace-masked groups of each integration. We note that a full readout of the detector mitigates bias drifts using reference pixels, but the subarray readouts used here do not have such pixels. We find that the superbias level changes by hundreds of parts per million (ppm) throughout the time series, with typical values of the scaling factor about 1.003. We use the standard-deviation-normalized time series of the superbias scaling coefficient as a detrending vector at the light-curve fitting stage, added linearly to our usual systematics model. We find that the superbias decorrelation coefficient is statistically preferred in the systematics model, with some residual structure in the photometry well explained by this term. The addition of superbias detrending reduced the transit depth tension between NRS1 and NRS2, with the white light-curve transit depths agreeing within the uncertainties.

For the Eureka! reduction, we also investigated time-dependent variations in the NRS2 detector bias level. We found that applying a scale factor correction to the superbias frame for each integration in Stage 1 marginally improved the consistency in measured transit depths (by ∼20 ppm) but also led to increased scatter. Applying a single-scale factor correction for all integrations yielded a similar improvement but without the increased scatter. We continue to investigate different methods of scaling the superbias frame. In the meantime, we elect to adopt the standard bias correction in our final Eureka! analysis and apply a manual offset of 78 ppm in transit depth to NRS2, Transit 1.

To account for NRS2 transit visit discrepancy for the final Tiberius reduction, we also manually offset the transmission spectrum for NRS2, Transit 1, by 63 ppm, such that the median transit depth was equal to NRS2, Transit 2.

After this superbias detrending in FIREFLy and manual offsets in Eureka! and Tiberius, we saw excellent agreement between the Eureka!, FIREFLy, and Tiberius spectra across both NRS1 and NRS2 in both transits, as shown in Figure 2. Since the superbias correction alters FIREFLy's absolute transit depths, we elect to compare their relative transit depths. In Tables 1, 2, and 3 we show the best fit parameters obtained for the Eureka!, FIREFLy, and Tiberius reductions, respectively.