We employ a modified version of the JExoRES pipeline (M. Holmberg & N. Madhusudhan 2023; N. Madhusudhan et al. 2023b), adapted to MIRI LRS. First, we perform Stages 1–2, using the JWST Science Calibration Pipeline (H. Bushouse 2020) to perform calibrations and the detector ramp fitting. In Stage 1, we perform the data quality initialization, electromagnetic interference (EMI) correction, saturation flagging, first and last frame flagging, linearity correction, reset switch charge decay (RSCD) correction, dark-current subtraction, and ramp fitting. For this, we use both the standard linearity and RSCD corrections and a custom approach detailed in Appendix A.2, as pursued by several recent works (D. Grant et al. 2023; E. M. R. Kempton et al. 2023; A. Dyrek et al. 2024). Overall, we find that the approach to nonlinearity correction has a minimal effect on the transmission spectrum of K2-18 b, as discussed in Appendix A.2 and illustrated in Figure 7. Next, we perform the flat-field correction in Stage 2.

In Stage 3, we first use the gain reference file to convert the flux from DN s^–1 to e⁻ s^–1. We then search for cosmic-ray hits by performing sigma clipping on the time series of each pixel by first removing a running median of seven integrations and using a threshold of 7σ. In addition, we also mask neighboring pixels surrounding detected outliers as well as pixels that have been flagged with any issues during Stages 1–2. We then correct for the background by subtracting the mean of the flux outside the trace, for each detector column and integration. To estimate the background, we use columns with pixel numbers 11–30 and 44–63, while not including bad pixels or outliers. Finally, we extract the spectra with optimal extraction (K. Horne 1986), using the median of all integrations as a model for the point-spread function for each spectral channel. For this, we use an aperture of 9 pixels. A box extraction approach produces very similar results, as shown in Figure 8. We reject additional outliers during optimal extraction and then in the light curves themselves, using sigma clipping with a running median filter.

Next, we fit the light curves using the batman (L. Kreidberg 2015) transit model and perform nested sampling with MultiNest (F. Feroz et al. 2009), assuming a circular orbit with a period from B. Benneke et al. (2019). We first construct a white-light curve by integrating the flux between 4.8 and 10 μm. Using this wavelength range gave slightly lower (5%) scatter compared to the full 4.8–12 μm range. For fitting the light curve, we use normally distributed priors on the normalized semimajor axis a/R_* and the inclination i from the weighted average of the NIRISS and NIRSpec fits from N. Madhusudhan et al. (2023b). These orbital parameters are also consistent with pre-JWST measurements (B. Benneke et al. 2019; A. Tsiaras et al. 2019). For the limb darkening, we use the quadratic law with the parameterization and priors by D. M. Kipping (2013). We use wide uniform priors on all other parameters. For the baseline trend, we adopt an exponential and a linear component: F_obs(t) = F_out (1 + ατ + γe^−τ/) F_transit(t), where τ is the time since the start of the observation, F_transit is the batman transit model, and F_out, α, γ, and are trend parameters. Furthermore, we mask the first 250 integrations (1034 s) to remove the strongest effect of the detector settling. Figure 1 shows the MIRI white-light curve of K2-18 b together with our model fit. Moreover, we also fit for an error inflation parameter to account for additional white noise. The parameters measured from the white-light-curve fitting are shown in Table 1.

Table 1. Parameters Estimated from the White-light-curve Analysis of Our JWST MIRI LRS Observation of K2-18 b

Parameter	JExoRES	JexoPipe
Mid-transit time, T0 (BJD – 2400000.5)
Inclination, i (deg)
Normalized semimajor axis, a/R*
Planet-to-star radius ratio, Rp/R*
First limb-darkening coefficient, u1
Second limb-darkening coefficient, u2

Note. The results are shown from the two independent data analyses. In both cases, the period was fixed to 32.940045 days (B. Benneke et al. 2019) and wavelengths between 4.8 and 10 μm were used.

Download table as:  ASCIITypeset image

Finally, we perform the spectroscopic light-curve fitting, fixing a/R_*, i, the mid-transit time, and the two limb-darkening coefficients to the values obtained from the white-light curve. We bin the light curves in wavelength prior to fitting, nominally with a width of 0.2 μm or 5 pixels, whichever contains the most pixels, as described in Appendix A.4. A minimum bin width of 0.2 μm is in line with the bin widths used in previous works, ranging between 0.15 and 0.5 μm (e.g., J. Bouwman et al. 2023; D. Grant et al. 2023; T. J. Bell et al. 2024; D. Powell et al. 2024). Due to instability of the spectrum below 5.6 μm as a function of binning, we disregard this part from further analysis, to be conservative, as discussed in Appendix A.4. We again mask the first 250 integrations. We explore different choices of trends in Appendix A.5 and find that an exponential and a linear component produce a stable transmission spectrum when varying the number of integrations to mask at the start. Furthermore, we find evidence for time-correlated noise on the timescale of minutes, with increasing amplitude toward shorter wavelengths. Accounting for this additional noise using Gaussian processes (GPs), as described in Appendix A.7, we find that the transit depth uncertainty increases by around 30% at 6 μm, and less at longer wavelengths. We also test different assumptions for the limb darkening and find that it does not significantly alter the transmission spectrum, as shown in Appendix A.6. The resulting MIRI transmission spectrum of K2-18 b is shown in Figure 2.

We use a second pipeline, JexoPipe , a JWST data reduction framework previously applied to NIRSpec Prism and G395H data (S. Sarkar et al. 2024) and adapted here for MIRI LRS. Starting with the .uncal files, in Stage 1 we utilize the official JWST Science Calibration Pipeline steps in the following sequence: group scale, data quality array initialization, EMI correction, saturation flagging, first frame, last frame, reset, linearity correction, RSCD correction, dark subtraction, jump detection, ramp fitting, and gain scale. Stage 1 ends with integration-level .rateints files. After combining all Stage 1 segments into a single .rateints file, we begin Stage 2 by applying the JWST calibration pipeline steps: assign World Coordinate System (WCS) and flat field. We then apply a custom bad-pixel flagging step (S. Sarkar et al. 2024). This step applies NaN values to all pixels that have abnormal Data Quality (DQ) flags as well as to 3σ outliers found on a row-by-row basis in each integration image. Next, we perform background subtraction per integration, using pixel columns (10)–(14) and (60)–(66). We apply an outlier mask, then subtract the mean of each row from all pixel values in that row.

We then apply a custom bad-pixel correction step using a combination of temporal and spatial interpolation, as described in S. Sarkar et al. (2024). Next, we detect the remaining outliers for each integration by comparison with a rolling median (of 20 contiguous integrations), replacing pixel values ±5σ from the rolling median with the median value. For Stage 3, 1D spectral extraction, for each integration, we apply an aperture of 9 pixel columns centered on the spectrum maximum and then perform optimal extraction (K. Horne 1986). We obtain the wavelength for each pixel row from the mean of the wavelengths in that row after application of the aperture.

We use the 1D spectra time series to construct a white-light curve (between 4.8 and 10 μm). We exclude the first 250 integrations, where the systematic trend is most extreme. We identify outliers on the white-light curve, ±2.5σ from a rolling median, and replace the 1D spectra corresponding to these outliers with linearly interpolated spectra from adjacent integrations. We scale the error bars on the light-curve points such that the average error bar equals the observed standard deviation of the scatter in the out-of-transit residuals. We use emcee (D. Foreman-Mackey et al. 2013) to perform a Markov Chain Monte Carlo parameter estimation of the white-light curve, fitting for a transit model with quadratic limb darkening generated by pylightcurve (A. Tsiaras et al. 2016) multiplied by a systematic trend consisting of an exponential term and a linear term (as in Section 2.1). In the white-light curve, we fit for R_p/R_*, the mid-transit time, a/R_*, i, quadratic limb-darkening coefficients, and four parameters for the trend. Uniform priors are used except for a/R_* and i, where we apply Gaussian priors based on values in N. Madhusudhan et al. (2023b) and use the Kipping parameterization (D. M. Kipping 2013) for limb-darkening priors. We fix the period to 32.940045 days (B. Benneke et al. 2019), the argument of periastron to 90°, and the eccentricity to 0. The white-light-curve parameter estimates are given in Table 1.

We bin the spectral light curves following the same prescription as in Section 2.1. We fit the spectral light curves in a similar manner to the white-light curve, with the exception of the mid-transit time, a/R_*, and i, which we fix to the white-light values. We also fix the limb-darkening coefficients to those from the white-light fit. To check the sensitivity of the transmission spectrum to the treatment of limb darkening, we also adopt fixed wavelength-dependent limb-darkening coefficients from ExoTic-LD (the Kurucz model; D. Grant & H. R. Wakeford 2024). In doing so, we find that the spectra from the two treatments are consistent to well within the 1σ uncertainties.

For each spectral light-curve fit, the free parameters are R_p/R_* and the four systematic trend parameters. To account for time-correlated noise, we obtain the β factor for each spectral light curve, using the “time-averaging” method to produce Allan deviation plots (J. N. Winn et al. 2007; P. Cubillos et al. 2017). The spectral light-curve fits are then rerun after inflating the uncertainties on the original light curves by the corresponding β value. The final transmission spectrum is shown in Figure 3.

The transmission spectra obtained from the two pipelines show good agreement, as shown in Figure 3, with all the data points agreeing within the 1σ uncertainties. While cross-checks between pipelines are useful to ensure consistency, it is also important to check the robustness and sensitivity of the transmission spectrum obtained using a given pipeline (M. Holmberg & N. Madhusudhan 2024). In the present work, we conduct robustness tests with both pipelines across various considerations in the data reduction and analyses. We discuss some of these results in Appendix A. We demonstrate the robustness of the resultant transmission spectra to a wide range of parameters and assumptions in data reduction and light-curve analyses. These include different treatments for the nonlinearity correction, spectral extraction, spectral binning, and limb darkening, and different trends in the light-curve fitting, as well accounting for correlated noise.

The resulting spectra are consistent within the uncertainties across all these cases, as shown in Appendix A. The corresponding atmospheric constraints are also consistent across these cases, as shown in Table 2. We also conduct atmospheric retrievals with spectra from both JExoRES and JexoPipe , to ensure consistency in the overall conclusions. The abundances and detection significances obtained from each reduction are summarized in Table 2. This is further discussed below.

Table 2. Retrieved Abundance Estimates and Detection Significances from the MIRI Transmission Spectrum of K2-18 b

Canonical Model	Data	DMS	DMDS	()	(B)	Detection Significance
DMS Only	JExoRES		⋯	215.45	2.86	DMS (2.9σ)
DMDS Only	JExoRES	⋯		216.40	3.81	DMDS (3.2σ)
DMS + DMDS	JExoRESa	<−2.44		216.40	3.80	DMS+DMDS (3.2σ)
DMS Only	JexoPipe		⋯	211.18	3.18	DMS (3.0σ)
DMDS Only	JexoPipe	⋯		211.18	3.17	DMDS (3.0σ)
DMS + DMDS	JexoPipe	−4.68		211.59	3.59	DMS+DMDS (3.2σ)
Effect of Trendb
DMS + DMDS	JExoRES (Exp+Linear1)a	<−2.44		216.40	3.81	DMS+DMDS (3.2σ)
DMS + DMDS	JExoRES (Exp+Linear2)	<−1.81		216.40	2.89	DMS+DMDS (2.9σ)
DMS + DMDS	JExoRES (Exp+Quadratic)	<−2.42		212.92	2.85	DMS+DMDS (2.9σ)
DMS + DMDS	JExoRES (Quadratic)			214.50	4.22	DMS+DMDS (3.4σ)
Masking of Potential Shadow Regionc
DMS + DMDS	JExoRES (λ < 10.5 μm)	<−2.10		168.44	3.51	DMS+DMDS (3.1σ)
DMS + DMDS	JExoRES (λ < 10 μm)	<−2.15		153.20	3.09	DMS+DMDS (3.0σ)
Accounting for Time-correlated Noise
DMS + DMDS	JExoRES (no GP)	<−2.16		215.94	4.08	DMS+DMDS (3.3σ)
DMS + DMDS	JExoRES (GP)a	<−2.44		216.40	3.81	DMS+DMDS (3.2σ)
Effect of Binningd
DMS + DMDS	JExoRES (0.2 μm)	<−2.16		215.94	4.08	DMS+DMDS (3.3σ)
DMS + DMDS	JExoRES (0.4 μm)	<−2.44		112.72	3.98	DMS+DMDS (3.3σ)
DMS + DMDS	JExoRES (0.8 μm)	<−2.87		63.71	3.79	DMS+DMDS (3.2σ)

Notes. The first column shows the canonical model used with DMS and/or DMDS included. The second column shows the spectrum used from one of the two pipelines, JExoRES or JexoPipe , along with any specific treatments in the data reduction or light-curve analysis. The DMS and DMDS abundances in the third and fourth columns are shown as of the volume mixing ratios. The retrieved abundances are reported as median values with 1σ uncertainties, for cases with a well-defined peak in the posterior, or as 95% upper limits otherwise. The model evidence is shown for each case as (), where is the Bayesian evidence. (B) refers to the natural logarithm of the Bayes factor, comparing the canonical model with a model with DMS and/or DMDS removed. The corresponding molecules are shown in the last column, with their detection significances in parentheses. We note that there is a typical uncertainty of ∼0.1σ in the calculation of the detection significance using MultiNest. GP indicates analysis using Gaussian processes to account for time-correlated noise, as outlined in Appendix A.7. ^aThese correspond to our canonical case presented in Figure 2. ^bThe trend refers to the systematic trend used in the light-curve analysis, as explored in Figure 10. For the two Exp+Linear cases, we remove the first 250 and 500 integrations for Linear1 and Linear2, respectively. ^cWe do not find a strong discontinuity in the detector behavior past 10.6 μm, called the shadow region, sometimes found in MIRI LRS data (T. J. Bell et al. 2024). Nevertheless, we conduct retrievals with different wavelength limits for robustness. ^dThese cases correspond to the spectra presented in Figure 9, which explore the effect of binning. These cases do not use the GP model. For the 0.2 μm case, we use a minimum bin width of 5 pixels, representing our nominal binning as discussed in Section 2.1.

Download table as:  ASCIITypeset image

The present transmission spectrum is the first mid-infrared (∼6–12 μm) transmission spectrum of a candidate hycean world ever observed. Given the wide range of molecules that could contribute spectral features in this range, we follow a hierarchical Bayesian approach to determine the preferred model for explaining the observations. We carry out the model selection in two stages, with the model parameters in the retrieval at the second stage being informed by the first stage, as discussed below. At the end, the preferred model given the data is chosen by comparing the Bayesian evidence between the two stages.

We first conduct a maximal retrieval using the JExoRES spectrum, with an atmospheric model that includes a substantially expanded set of molecules over previous atmospheric retrievals of K2-18 b. We consider 20 molecules, including 11 molecules considered in previous work (N. Madhusudhan et al. 2023b) and nine new ones. The previous molecules include H₂O, CH₄, NH₃, CO, CO₂, HCN, (CH₃)₂S (DMS), CH₃Cl, CS₂, OCS, and N₂O. The molecules H₂O, CH₄, NH₃, CO, CO₂, and HCN have been predicted to be the prominent C-, O-, and N-bearing species in temperate H₂-rich atmospheres, including hycean worlds (e.g., R. Hu 2021; N. Madhusudhan et al. 2021; X. Yu et al. 2021). The remaining molecules (CH₃)₂S (DMS), CH₃Cl, CS₂, OCS, and N₂O have been predicted as observable potential biosignatures in H₂-rich atmospheres of habitable super-Earths (A. Segura et al. 2005; S. D. Domagal-Goldman et al. 2011; S. Seager et al. 2013a) and hycean worlds (N. Madhusudhan et al. 2021).

The new molecules include other prominent sulfur-based molecules, hydrocarbons, and potential biomarkers (e.g., C. Sousa-Silva et al. 2020; E. W. Schwieterman & M. Leung 2024; S.-M. Tsai et al. 2024): H₂S, SO₂, C₂H₂, C₂H₄, C₂H₆, (CH₃)₂S₂ (dimethyl disulfide; DMDS), CH₃SH, PH₃, and CH₃OH. Of these molecules, DMDS, CH₃SH, C₂H₆, and PH₃ are particularly relevant as potential biomarkers both in hycean as well as terrestrial-like atmospheres, as suggested in recent studies (e.g., C. Sousa-Silva et al. 2020; E. W. Schwieterman & M. Leung 2024; S.-M. Tsai et al. 2024). Of all the molecules considered, DMS, DMDS, CH₃SH, CH₃Cl, and N₂O are expected to be the most reliable biosignatures in hycean conditions, whereas others may have alternate abiotic sources (A. Segura et al. 2005; S. D. Domagal-Goldman et al. 2011; S. Seager et al. 2013a; E. W. Schwieterman & M. Leung 2024; S.-M. Tsai et al. 2024).

The retrieval obtains a good fit to the spectrum, with the model being inconsistent with a featureless spectrum, i.e., a flat line, at 3σ significance. However, an inspection of the posteriors reveals only one well-constrained peak for the molecule DMDS, while all the other molecules are unconstrained, as shown in Appendix B. We also do not find significant constraints for cloud/haze parameters that could affect the mid-infrared spectrum.

We further conduct a retrieval without DMDS and find a well-constrained peak in the posterior for DMS. We find that DMS replaces DMDS, due to the strong degeneracy between the two molecules in the MIRI wavelength range where prominent peaks are evident in the data, as shown in Figure 2. The spectral contributions due to the two molecules are shown in Figure 4. Thus, the retrievals highlight the degeneracy between DMS and DMDS, with a higher preference for DMDS from this data set. On removing both molecules from the retrieval, we find no significant constraints on any of the remaining molecules. The model without DMS and DMDS does not provide a good fit to the data, with the maximal model preferred over this model at 2σ. Moreover, the model without DMS and DMDS is only marginally favored over a flat spectrum below 2σ significance. The data therefore show notable evidence for only two molecules, DMDS and DMS.

Informed by the above retrieval, we arrive at a canonical model that includes DMS and DMDS. We additionally include CH₄ and CO₂, which have been detected in K2-18 b previously in the near-infrared (N. Madhusudhan et al. 2023b). In total, our canonical model includes 16 free parameters: four molecular mixing ratios, six for the P–T profile, four for the cloud parameters, one for P_ref, and one for the offset on the MIRI data. The retrieval with this model provides a comparable fit to the data as the maximal model above, but with a higher Bayesian evidence, due to the lower number of unconstrained free parameters. The canonical model is preferred over the maximal model by 2.2σ and over a featureless spectrum, i.e., a flat line, by ≥3.4σ. For computing the Bayesian evidence for the one-parameter model of a featureless spectrum, we assume uniform priors on the model transit depth, with a range of ±N ppm relative to the measured white-light transit depth. Considering the span of transit depths in the observed spectrum of ∼600 ppm, we explore values of N between 600–2000 ppm. Over this range, the canonical model is preferred over a featureless spectrum at 3.4σ–3.8σ, with higher preferences obtained for wider priors, i.e., larger N. We adopt the minimum preference of 3.4σ as a conservative estimate.

As with the maximal model, we find that only DMDS is constrained by this data set, while other parameters remain unconstrained. We also recover the degeneracy between DMDS and DMS, such that when DMDS is excluded from the model, DMS is recovered. In what follows, we present the detection significances of the molecules and atmospheric constraints with the canonical model using the spectra obtained from both pipelines.

The spectrum provides important insights into the atmospheric composition of K2-18 b. The retrieved spectral fit along with the 1σ and 2σ contours using the JExoRES spectrum is shown in Figure 2. As noted above, the present MIRI spectrum spans a wavelength range of ∼5.8–12 μm, which encompasses strong spectral features of several prominent molecules expected in temperate H₂-rich atmospheres, as shown in Figures 4 and 14. However, as discussed above, we are unable to explain the observed features with most of the 20 prominent species considered, with the maximal retrieval finding no significant evidence for 18 of the 20 species. The only species with notable evidence are DMDS and DMS. Both molecules have a similar double-peak feature between 6.8 and 8 μm, with additional broad peaks around 9.8 μm and 10.5 μm for DMS and DMDS, respectively, which are consistent with the observed spectrum. Additionally, the amplitudes of the observed spectral features, of ∼300–400 ppm, in the present MIRI data are significantly larger compared to the ∼200 ppm amplitudes seen in previous near-infrared (∼1–5 μm) observations with the NIRISS and NIRSpec instruments (N. Madhusudhan et al. 2023b). While previous observations detected strong spectral features of CH₄ and CO₂ in the near-infrared, neither of them are detectable in the present spectrum, as can be seen in Figure 4. While CO₂ does not have strong spectral features in the MIRI LRS range, those of CH₄ are dwarfed by stronger contributions due to DMS and/or DMDS. As such, the contributions of both CO₂ and CH₄ at previously reported mixing ratios (∼1%) are insufficient to explain the present data in the MIRI band, as can be seen in Figure 4. Instead, contributions from DMDS and DMS at volume mixing ratios of ∼10⁻⁵–10⁻³ readily explain the data, due to their strong absorption cross sections in the MIRI band.

The spectrum provides notable constraints on the presence and abundance of DMDS and DMS in the atmosphere, as shown in Figure 5. Due to the degeneracy between the spectral features of DMDS and DMS in the MIRI band, as shown in Figure 4, it is difficult to robustly distinguish between them at the current data quality. Considering the canonical model, where both DMDS and DMS are present, DMDS is preferred at a somewhat higher significance (∼2σ) compared to DMS for the JExoRES spectrum. The presence of DMDS in the model is favored at 2σ significance compared to the model without DMDS. DMS is not well constrained, and the model with DMS is preferred over the model without DMS at only 1σ significance. The lack of either molecule in the model is compensated for by the other, due to the degeneracy in their spectral features, as discussed above. However, the combination of DMDS and DMS together is favored at 3.2σ over a model with neither molecule included (Table 2).

On the other hand, retrievals with models including only one of the two molecules at a time show relatively high detection significances for both molecules individually. For example, considering the canonical model with only DMDS included, i.e., DMS removed, we find that the detection significance for DMDS is 3.2σ. Similarly, the canonical model with only DMS included, i.e., DMDS removed, provides a detection significance for DMS of 2.9σ. Therefore, the present data provide evidence for the presence of DMDS and/or DMS at 2.9σ–3.2σ significance.

The retrieved posterior distributions for the volume mixing ratios of DMDS and DMS shown in Figure 5 reflect the degeneracy between the two species that is apparent from their spectral contributions, as discussed above. When both molecules are included in the model, the spectral contribution of DMDS dominates over that of DMS. The posterior of DMDS shows a clear peak, albeit with a low-abundance tail, due to the degeneracy with DMS. The abundance of DMDS is retrieved to be = −3.48, while DMS is unconstrained. Here, X_DMDS is the volume mixing ratio of DMDS. In retrievals with only one of the two molecules present, the posterior of that molecule shows a single well-constrained peak with no low-abundance tail. In such retrievals, the abundances are retrieved to be = −3.25 and = −3.42, as shown in Table 2. Last, we repeat the canonical retrieval using the UltraNest nested sampling implementation and find consistent abundance estimates and detection significances to our standard retrievals using the MultiNest implementation.

We also conduct another canonical retrieval with our other spectrum obtained with the JexoPipe pipeline. We find that the detection significances and abundance estimates of DMS and DMDS are consistent with those obtained using the JExoRES spectrum, as shown in Table 2. We find the combined detection significance of DMDS and DMS to be 3.2σ, similar to that obtained with the canonical retrieval with the JExoRES data. The abundances are slightly lower but still consistent. For this canonical model, where DMS and DMDS are both included, the posteriors for DMS and DMDS are comparable, as shown in Figure 5. Using the individual DMS-only or DMDS-only cases, we find evidence for each of the molecules at 3σ significance and similar abundance estimates to the corresponding JExoRES cases: = and = . Therefore, across all our retrievals, we detect DMDS and/or DMS at 3σ significance.

The observations do not provide significant constraints on any of the other molecules or atmospheric properties but are consistent with previous inferences from the near-infrared JWST spectrum (N. Madhusudhan et al. 2023b). Specifically, we do not obtain significant constraints on CH₄ or CO₂, which have previously been detected. However, the upper limits we retrieve for the abundances of both molecules are consistent with their previous abundance estimates. The contributions from both molecules at previously measured abundances of ∼1% each are insufficient to explain the observed spectral features in the present spectral range with JWST MIRI (6–12 μm). While CO₂ has no strong spectral features in this range, those of CH₄ are significantly weaker compared to those of DMDS and/or DMS, as shown in Figure 4. The abundances we retrieve for DMS and/or DMDS are also consistent with the constraints for DMS reported in previous work.

Similarly, we do not find significant evidence for spectral contributions due to clouds/hazes in the MIRI band and none of the cloud parameters are constrained. We find the photospheric pressure corresponding to the white-light radius to be (P_ref/bar) = , with no significant offset to the data retrieved. The offset is retrieved to be 12 ppm, which is consistent with zero within the 1σ uncertainties. The photospheric temperature at 1 mbar at the terminator is also weakly constrained to be K, which is somewhat higher but consistent with previous estimates in the near-infrared to within 1σ uncertainties.

For additional robustness, we conduct retrievals on the transmission spectra obtained using a wide range of assumptions for the data reduction and analyses. Using JExoRES , we explore differences in the treatment of systematic trends and time-correlated noise, as well as different wavelength limits and the effect of spectral binning, as discussed in Appendix A. Across all of these cases, we find consistent abundance constraints and detection significances of DMDS and/or DMS, ranging from 2.9σ to 3.4σ, as outlined in Table 2.

As a further robustness check for our findings, we carry out leave-one-out cross-validation using the JExoRES data. We follow A. Vehtari et al. (2017) and L. Welbanks et al. (2023), calculating the expected log pointwise predictive density (elpd) for each data point. For each data point d_i, we carry out an atmospheric retrieval with all other points considered and that data point removed. We do not make any approximations, e.g., Pareto-smoothed importance sampling, and instead directly compute the elpd of each data point by carrying out a retrieval considering all other data points. Using the retrieved posterior samples, we compute the probability density of the ith point d_i, given a retrieval using model on all other data d_−i, resulting in a total of S samples of the parameter space for each retrieval. The elpd metric provides an indication of how well each data point can be predicted by the model and therefore quantifies the performance of the model. The metric is given by

where the parameters corresponding to the jth sample in each retrieval are denoted as θ_j. We perform the above calculation for two atmospheric models: the canonical model described in Section 3 and one with DMS and DMDS absent. As such, we run a total of two retrievals per data point. The difference between the elpd values of the two models, Δelpd_i, indicates their relative performance. Conversely, in this case, Δelpd_i serves as an indication for the relative contributions of the data points toward the present detection of DMDS and DMS. The pointwise Δelpd_i is shown in Figure 6. We find that most points have Δelpd_i values above zero, indicating that the detection of DMDS and DMS combined is driven by multiple data points to various extents. Δelpd_i is highest for two data points near 7 μm and another two points just under 9 μm. This indicates that these points are better predicted by a model that includes DMDS and DMDS than a model without the two molecules. Conversely, the data point below 6 μm has the most negative Δelpd_i, indicating that it is better predicted by a model without DMS and DMDS. Overall, our leave-one-out analysis indicates that the present detection is unlikely to be the result of localized systematics or individual outlier points.

Our findings provide new independent evidence for the possibility of a biosphere on K2-18 b. As discussed in Section 1, the detection of carbon-bearing molecules, CH₄ and CO₂, and nondetections of NH₃ and CO (N. Madhusudhan et al. 2023b) are consistent with prior predictions for a hycean world and inconsistent with mini-Neptune or gas dwarf scenarios (G. J. Cooke & N. Madhusudhan 2024; F. E. Rigby et al. 2024). The previous observations also showed tentative hints of DMS, a possible biosignature molecule in the atmosphere (N. Madhusudhan et al. 2023b), but the evidence was weak. This new independent evidence for the presence of DMS and/or DMDS at 2.9σ–3.2σ adds to the tentative inference of DMS reported previously and bolsters the chances of a biosphere on K2-18 b.

Recent photochemical modeling of K2-18 b suggests that large quantities of DMS and DMDS, with mixing ratios up to 10⁻², can accumulate in the atmosphere for assumed biogenic oceanic fluxes of these gases of ≳20 times Earth levels (S.-M. Tsai et al. 2024). Such quantities are consistent with our current MIRI results for DMDS and DMS and with previous abundance constraints of DMS (N. Madhusudhan et al. 2023b). These otherwise photochemically fragile molecules survive in the S.-M. Tsai et al. (2024) models under high-surface-flux-emission conditions: (1) because of the favorable ultraviolet spectral energy distribution of M dwarf stars (see also A. Segura et al. 2005; S. D. Domagal-Goldman et al. 2011; S. Seager et al. 2013a); (2) because DMS and DMDS can eventually self-shield, as well as shield molecules such as CO₂ from photolysis, when their column density is high enough, significantly cutting down on their photochemical destruction pathways (S.-M. Tsai et al. 2024); and (3) because the models assume that biogenically emitted DMS and DMDS do not have a return sink to the ocean, allowing the abundance to build up over time. Our present MIRI abundance estimates of DMS and/or DMDS at mixing ratios of ∼10⁻⁵–10⁻³ could, therefore, imply stronger biological activity on the planet in comparison to Earth. We note however, that DMDS could potentially condense on K2-18 b (see C. Sagan & B. N. Khare 1971; B. N. Khare et al. 1978; D. M. VonNiederhausern et al. 2006), depending on tropopause temperatures, restricting its gas-phase abundance in the stratospheric region that is probed in the transit observations. Therefore, abundances significantly higher than ∼10⁻⁴–10⁻³ may be difficult to achieve if the stratosphere is too cold.

On the other hand, it is possible that the abundance estimates derived in our work are strongly influenced by the uncertainties in spectral parameters and cross sections of these molecules used in the models. The derived abundance and temperature are strongly dependent on the absorption cross section of a molecule detected. For both DMDS and DMS, only limited cross sections are available in the literature (I. Gordon et al. 2017), which are derived experimentally, assuming an Earth-like N₂-rich background atmosphere at nearly STP conditions. It is possible that the collisional broadening factors and cross sections may be different for an H₂-rich background gas at lower pressures probed in the present transmission spectroscopy. Therefore, our results highlight the acute need for laboratory and theoretical work to derive high-fidelity absorption cross sections for these and other biosignature molecules, to enable their robust detection and abundance estimates in habitable exoplanets.

As with any potential biosignature molecule, it is important to ask whether there can be abiotic sources of DMS and DMDS in a temperate H₂-rich atmosphere that could explain the observations. Both DMS and DMDS have been predicted to be robust biosignatures for Earth-like planets as well as planets with H₂-rich atmospheres, including super-Earths and hycean worlds (S. D. Domagal-Goldman et al. 2011; S. Seager et al. 2013a; D. C. Catling et al. 2018; E. W. Schwieterman et al. 2018; N. Madhusudhan et al. 2021). Their identification as robust biosignatures is due to the fact that on Earth, both molecules are uniquely produced by life (particularly marine biota) in small quantities of ≲1 ppb by volume and are not supplied by abiotic photochemistry. Nevertheless, here we explore some potential alternatives.

Experimental studies have demonstrated the feasibility of forming several organosulfur compounds abiotically, using ultraviolet irradiation or electric discharges of gaseous mixtures containing H₂S and CH₄ (e.g., F. Raulin & G. Toupance 1975; C. He et al. 2020; V. Vuitton et al. 2021; N. W. Reed et al. 2024). In particular, some of these studies have demonstrated the production of both DMS (F. Raulin & G. Toupance 1975; N. W. Reed et al. 2024) and DMDS (C. Sagan & B. N. Khare 1971; B. N. Khare et al. 1978) in gas mixtures containing both CH₄ and H₂S, arguing in favor of their possible abiotic production in reduced planetary atmospheres. However, both DMS and DMDS are highly reactive and have very short lifetimes in the above experiments (e.g., a few minutes) and in the Earth’s atmosphere (between a few hours to ∼1 day), due to various photochemical loss mechanisms (e.g., S. Seager et al. 2013b). Thus, the resulting DMS and DMDS mixing ratios in the current terrestrial atmosphere are quite small (typically ≲1 ppb), despite continual resupply by phytoplankton and other marine organisms.

Therefore, sustaining DMS and/or DMDS at over 10–1000 ppm concentrations in a steady state in the atmosphere of K2-18 b would be implausible without a significant biogenic flux. Moreover, the abiotic photochemical production of DMS in the above experiments requires an even greater abundance of H₂S as the ultimate source of sulfur—a molecule that we do not detect—and requires relatively low levels of CO₂ to curb DMS destruction (N. W. Reed et al. 2024), contrary to the high reported abundance of CO₂ on K2-18 b (N. Madhusudhan et al. 2023b).

Another recent study has reported evidence for the presence of DMS on the comet 67P/Churyumov–Gerasimenko (N. Hänni et al. 2024), motivating the possibility of a potential abiotic source through cometary delivery to the exoplanet’s atmosphere. Comets are also known to contain other ices and sublimated gases that might be considered biosignatures if present in an exoplanet atmosphere, including O₂, CH₄, and CH₃Cl (A. Bieler et al. 2015; E. C. Fayolle et al. 2017; M. Rubin et al. 2019). However, cometary delivery is implausible for explaining such molecules in planetary atmospheres, due to the insufficient contribution that trace species within a comparatively small comet would add to the much more massive planetary atmosphere (e.g., R. W. Court & M. A. Sephton 2012; R. C. Felton et al. 2022; M. Leung et al. 2022). Moreover, molecules such as DMS and DMDS would be shock-heated during a hypervelocity cometary impact and its related plume splashback phase, which would reset the bulk composition to simpler molecules that are more stable in thermochemical equilibrium at high temperatures and (for the splashback phase) low pressures (K. Zahnle 1996). Therefore, DMS and DMDS would not be delivered at any significant measurable levels through such a process. We also note that DMDS has not been reported in the comet 67P. In summary, cometary impacts cannot be responsible for the relatively large DMS and/or DMDS abundances inferred from the present MIRI retrievals.

A conclusive identification of a biosignature necessitates a robust assessment of various factors, including the robustness of the detections, the environmental context, and potential false positives (e.g., E. W. Schwieterman et al. 2018; V. Meadows et al. 2022). It is widely recognized that the detection of a biosignature is unlikely to be instantaneous or unambiguous in the first instance, rather relying on continued accumulation of evidence and addressing the above factors (V. Meadows et al. 2022). Our study is the first formal step in this direction, building on the first possible hints of DMS reported in our previous work (N. Madhusudhan et al. 2023b) and the further evidence of DMDS and/or DMS with a higher significance seen in the present observations. Finally, our present detection, based on multiple spectral features with a different instrument in a different spectral range from previous work, provides an important independent line of evidence in this direction.

Further work is needed to robustly verify the current findings. More observations are required to robustly demonstrate the repeatability of our present findings, rule out potentially unaccounted-for instrumental systematics, as well as increase the detection significances. While DMDS and DMS best explain the current observations, their combined detection significance is ∼3σ, which is at the lower end of the robustness typically required for scientific evidence. The significance can be readily increased to a 4σ–5σ level by a modest amount of additional JWST time, e.g., between one and three additional transits with MIRI, i.e., only ∼8–24 hr. Second, while we have explored 20 prominent molecules in fitting the spectrum, our search may still not be fully exhaustive. Therefore, future studies could investigate other potential molecules that could explain the data. At the same time, as discussed above, new experimental and theoretical studies are needed to determine accurate absorption cross sections for DMDS and DMS, and other potential biosignature molecules, for conditions relevant to candidate hycean worlds like K2-18 b. Future laboratory experiments and/or theoretical modeling are also needed to fully explore the possible photochemical mechanisms for producing DMS and DMDS in dry, methane-rich, reduced environments, to address potential abiotic sources of these molecules.

Overall, our findings present an important step forward in the search for signatures of life on exoplanets. However, robustly establishing both the veracity of the present findings and their possible association with life on K2-18 b needs a dedicated community effort in multiple directions—observational, theoretical, and experimental. Observations with JWST are already demonstrating that possible hycean worlds indeed significantly expand and accelerate the search for life elsewhere. The central question now is whether we are prepared to identify the signatures of life on such worlds. The opportunity is at our doorstep.

In this work, we thoroughly test several data reduction methods and assumptions—such as the effect of the nonlinearity correction, the choice of binning, detrending, and limb darkening—to produce a robust transmission spectrum. We also apply two independently developed data reduction pipelines to check for consistency, using JExoRES and JexoPipe , as described in Sections 2.1 and 2.2, respectively. Figure 3 demonstrates the stability of the MIRI transmission spectrum of K2-18 b when comparing these two pipelines. Table 1 provides the parameters from the white-light-curve fitting with both pipelines, yielding results consistent within 1σ.

We explored how differences in nonlinearity correction can affect the transmission spectrum, as pursued in recent works (D. Grant et al. 2023; E. M. R. Kempton et al. 2023; A. Dyrek et al. 2024). As the baseline case, we use the standard nonlinearity corrections, RSCD correction, and dark subtraction. Note that the RSCD correction step does not alter the raw data. Instead, it applies a mask to the first few groups—removing these from the ramp fitting. On the other hand, the dark reference file includes both the dark signal and the RSCD effect, which are corrected by applying the dark subtraction step. For our custom approach, we empirically derive our own RSCD and nonlinearity correction. Using the average of the last 500 (out-of-transit) integrations of our observation, accounting for outliers, we first fit the ramp of each pixel using

where t is the time since the start of an integration. Here, p₁ corresponds to an overall offset/bias, p₂ is the linear part of the signal that we use to represent the true flux (J. E. Morrison et al. 2023), p₃ is the quadratic nonlinear component representing the change in response as the charge accumulates, p₄ is the strength of the RSCD signal, and p₅ is the timescale of the RSCD signal. For this fit, we exclude the first and last group, as done in the main analysis described in Sections 2.1 and 2.2. Next, we subtract the RSCD correction, including the constant term, i.e.,

from the data before applying our custom nonlinearity correction. We derive the parameters for the nonlinearity correction using a fourth-degree polynomial h, defined as

where a_i, being pixel-dependent, is used to populate the custom nonlinearity reference file. Because we have already subtracted the RSCD signal, we do not apply the dark subtraction step in this scenario.

Following this custom nonlinearity correction approach, we perform the rest of the analysis as outlined in Section 2.1. Figure 7 shows the transmission spectrum of K2-18 b using the standard and custom nonlinearity correction. Overall, we find that the spectrum is largely insensitive to the choice of nonlinearity correction. We note that other data sets may be more sensitive to the nonlinearity correction if the flux level is higher (for the present observation, it is below 51% saturation) or for targets with a larger transit depth than K2-18 b.

Previous work has shown that transmission spectra obtained with MIRI LRS can be excessively noisy when using small spectroscopic bins (T. J. Bell et al. 2024). Therefore, recent works have adopted wide bins with typical widths of 0.15–0.5 μm (e.g., T. J. Bell et al. 2024; D. Powell et al. 2024; E. Schlawin et al. 2024; L. Welbanks et al. 2024). In this work, we test the sensitivity of our transmission spectrum of K2-18 b to different bin sizes, as shown in Figure 9. We adopt four different bin widths: a maximum of 4 pixels or 0.2 μm (whichever contains the most pixels), a maximum of 5 pixels or 0.2 μm, 0.4 μm, and 0.8 μm. We limit the bin widths to at least 4 or 5 pixels to average over potential systematic effects. Note that we only use whole pixel rows, to avoid introducing correlations between different spectral channels. We find that our transmission spectrum is generally consistent between the different bin widths, apart from a small region below around 5.6 μm. Specifically, in this region, the spectrum is not consistent between the 4 and 5 pixel bin widths, i.e., the green and gray data points in Figure 9. For this reason, we do not include the data below 5.6 μm in our atmospheric retrievals. Nominally, we select a maximum of 5 pixels or 0.2 μm as our canonical binning for both JExoRES and JexoPipe . We also conduct atmospheric retrievals with varying bin widths using our canonical model. We find consistent results between the retrievals with detection significances for DMS+DMDS between 3.2σ and 3.3σ across the different bin widths, as shown in Table 2.

MIRI time-series observations are affected by detector settling, in which there is an exponential-like ramp in the light curves at the beginning of the observation (J. Bouwman et al. 2023; T. P. Greene et al. 2023). This effect can be seen in the white-light curve of our observation, as shown in Figure 1. A common approach to treating this effect is to mask some number of integrations at the start of the observation and/or to account for the detector settling when fitting the light curves (e.g., J. Bouwman et al. 2023; D. Grant et al. 2023; T. P. Greene et al. 2023; E. M. R. Kempton et al. 2023; T. J. Bell et al. 2024; A. Dyrek et al. 2024)—for example, by including an exponential trend. For this reason, we investigate the impacts of different treatments of this effect on our MIRI transmission spectrum of K2-18 b. To model the light curves, we use the following systematic model:

where τ is the time since the start of the observation, F_transit is the batman transit model, and F_out, α, β, γ, and are trend parameters. We consider three different trend prescriptions: an exponential and a linear trend (β = 0, as described in Section 2.1), an exponential and a quadratic trend, and a quadratic trend alone (γ = 0). Note that for the white-light curve, we obtain a good fit to the data using the first option. Furthermore, we nominally remove 250 integrations (17 minutes) at the start of the observation, to remove the worst part of the detector settling effect. However, for the quadratic trend, we mask 500 integrations (34 minutes) at the start, to remove even more of the settling effect. Figure 10 shows the MIRI transmission spectrum of K2-18 b obtained using different trends, showing general agreement between these scenarios. For our canonical case with an exponential and a linear trend, we also consider a wider range of the number of masked integrations, between 125 and 750 integrations, and find that the spectrum is consistent across these cases (as shown in Figure 10 for 250 and 500 masked integrations).

In our canonical data reduction, for both JExoRES and JexoPipe , we fix the quadratic limb-darkening coefficients to the values estimated from the corresponding white-light curves, as given in Table 1. This is motivated by the small and approximately wavelength-independent limb darkening in the mid-infrared, between 5 and 12 μm. To test this assumption, we generate a transmission spectrum of K2-18 b using wavelength-dependent model limb-darkening coefficients obtained via ExoCTK (M. Bourque et al. 2021). For this, we use an ATLAS9 model (F. Castelli & R. L. Kurucz 2003)—predicting u₁ = 0.028 ± 0.008 and u₂ = 0.107 ± 0.011 when integrating over the MIRI LRS bandpass (5–10 μm)—in agreement with our white-light-curve estimates. Figure 11 shows that the transmission spectrum of K2-18 b does not significantly differ between our two approaches, highlighting the robustness of the spectrum with respect to the effect of limb darkening.

We find evidence for time-correlated noise when analyzing the light curves. By inspecting the autocorrelation functions as a function of wavelength, we find that the time-correlated noise has a timescale of a few minutes and that the amplitude of this noise increases toward shorter wavelengths. To account for this, we fit a model of the noise to the residuals of the light curves after subtracting the fitted transit models and while masking the first hour of the observation. Based on the exponential-like shape of the autocorrelation functions, we employ the following wavelength-dependent GP kernel:

where σ_λ is the standard deviation of a spectral channel with central wavelength λ, ρ(λ) is a smooth function describing the amplitude of the time-correlated noise, and η is the timescale of this noise (assumed to be wavelength-independent). We parameterize ρ(λ) as a cubic polynomial, defined by four points that are equally spaced in wavelength. This gives us five parameters in total to describe the kernel. To estimate these parameters, we add up the log-likelihoods from each spectral channel, as computed via celerite (D. Foreman-Mackey et al. 2017), and use MultiNest (F. Feroz et al. 2009) to obtain the posterior distributions. Figure 12 shows the estimated ρ(λ) together with the average autocorrelation function.

Next, we use the median ρ(λ) and η to refit the light curves with the above kernel, again using celerite to evaluate the likelihood. For this, we do not use an error inflation parameter, as that is already taken into account by σ_λ. We find that the center points of the spectrum (that is, the median values of the transit depth posterior distributions) remain very similar to the case without using the kernel (Equation (A5)), while the transit depth uncertainties increase by around 0%–30%, reaching a maximum at the shortest wavelengths. Figure 13 shows the transit depth uncertainties with and without our GP model. Given these findings, we use the GP case as our canonical JExoRES spectrum, which has more conservative error bars.

As an additional robustness check, we also account for time-correlated noise using a different method in the JexoPipe pipeline, applying the “time-averaging” method (J. N. Winn et al. 2007; P. Cubillos et al. 2017). Here, we obtain Allan deviation plots for the residuals from each spectral light curve and find the average β factor for bins between 5 and 20 minutes duration, where β is the ratio between the observed fractional noise and the expected value for uncorrelated noise. Including only β values ≥1, we find a median β of 1.14 across all the spectral bins used. To account for correlated noise, we then inflate the errors on the light-curve data points by these factors and repeat the light-curve fits. This gives the final JexoPipe spectrum used for retrieval analysis with conservative error bars. Figure 13 shows the final transit depth uncertainties.

Table 3. Retrieved Atmospheric Parameters and Corresponding Prior Probability Distributions for the Canonical Retrieval with the JExoRES Spectrum

Parameter	Bayesian Prior	Description	DMS + DMDS	DMS Only	DMDS Only
	(−12, −0.3)	Mixing ratio of DMDS		⋯
	(-12, −0.3)	Mixing ratio of DMS			⋯
	(−12, −0.3)	Mixing ratio of CH4
	(−12, −0.3)	Mixing ratio of CO2
T0/K	(100, 500)	Temperature at 0.1 μbar
T1mbar/K	⋯	Temperature at 1 mbar
	(0.02, 2.00)	P–T profile curvature
	(0.02, 2.00)	P–T profile curvature
	(−6, 1)	P–T profile region limit
	(−6, 1)	P–T profile region limit
	(−2, 1)	P–T profile region limit
	(−6, 0)	Reference pressure at RP
	(−4, 10)	Rayleigh enhancement factor
γ	(−20, 2)	Scattering slope
	(−6, 1)	Cloud-top pressure
ϕ	(0, 1)	Cloud/haze coverage fraction
δMIRI/ppm	(−100, 100)	MIRI data set offset

Note. The quoted parameter values denote the retrieved median and 1σ intervals. The T_1mbar values are inferred from the retrieved P–T profile parameter constraints.

Download table as:  ASCIITypeset image