Water vapour in the atmosphere of the habitable-zone eight Earth-mass planet K2-18 b

Angelos Tsiaras, Ingo P. Waldmann, Giovanna Tinetti, Jonathan Tennyson, Sergey N. Yurchenko

References

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreements 758892, ExoAI; 776403/ExoplANETS A) and under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ ERC grant agreement numbers 617119 (ExoLights) and 267219 (ExoMol). We furthermore acknowledge funding by the Science and Technology Funding Council (STFC) grants: ST/K502406/1 and ST/P000282/1. The data used here were obtained by the Hubble Space Telescope as part of the 13665 and 14682 GO proposals (PI: Björn Benneke).

Author Contributions

A.T. performed the data analysis and developed the HST analysis software Iraclis; I.P.W developed the atmospheric retrieval software TauREx; G.T. contributed to the interpretation of the results; J.T. and S.N.Y. coordinate the ExoMol project. All authors discussed the results and commented on the manuscript.

Financial Interests

The authors declare no competing financial interests.

Materials & Correspondence

Correspondence and requests for materials should be addressed to A.T. (atsiaras@star.ucl.ac.uk) or I.P.W. (ingo@star.ucl.ac.uk).

Methods

Observations

Nine transits of K2-18 b were observed as part of the HST proposals 13665 and 14682 (PI: Björn Benneke) and the data are available through the MAST Archive. More specifically, the relevant HST visits are: visit 29 (06/12/2015), visit 35 (14/03/2016), and visit 30 (19/05/1016) from proposal 13665; visit 3 (02/12/2016), visit 1 (04/01/2017), visit 2 (06/02/2017), visit 4 (13/04/2017), visit 5 (30/11/2017), and visit 6 (13/05/2018) from proposal 14682. Out of these nine visits, we decided not to include the last one, as it suffered from pointing instabilities.

Each transit was observed during five HST orbits, with the G141 infrared grism of the WFC3 camera (1.1 - 1.7 $\mu$ m), in the spatially scanning mode. During an exposure using the spatial scanning mode the instrument slews along the cross-dispersion direction, allowing for longer exposure times and increased signal-to-noise ratio (S/N), without the risk of saturation . Both forward (increasing row number) and reverse (decreasing row number) scanning were used for these observations.

The detector settings were: SUBTYPE=SQ256SUB, SAMP_SEQ=SPARS10, NSAMP=16, APERTURE=GRISM256, and the scanning speed was 1.4 ′′s-1. The final images had a total exposure time of 103.128586 seconds, a maximum signal level of $1.9\times 10^{4}$ electrons per pixel, and a total scanning length of approximately 120 pixels. Finally, for calibration reasons, a 0.833445 s non-dispersed (direct) image of the target was taken at the beginning of each visit, using the F130N filter and the following settings: SUBTYPE=SQ256SUB, SAMP_SEQ=RAPID, NSAMP=4, APERTURE=IRSUB256.

Extracting the planetary spectrum

We carried out the analysis of the eight K2-18 b transits using our specialised software for the analysis of WFC3, spatially scanned spectroscopic images, which has been integrated into the Iraclis package (see Code Availability section). The reduction process included the following steps: zero-read subtraction, reference pixels correction, nonlinearity correction, dark current subtraction, gain conversion, sky background subtraction, calibration, flat-field correction, and bad pixels/cosmic rays correction.

We extracted the white (1.088 – 1.68 $\mu$ m) and the spectral (Supplementary Table 1) light curves from the reduced images, taking into account the geometric distortions caused by the tilted detector of the WFC3/IR channel. The wavelength range of the white light curve corresponds to the edges of the WFC3/G141 throughput (where the throughput drops to 30% of the maximum). In addition, we tested two wavelength grids for the spectral light curves, with a resolving power of 20 and 50. We decided to use the latter as it was able to capture the observed water feature more precisely – i.e there were enough data points within the wavelength range of the water feature to produce a statistically significant result.

We fitted the light curves using our transit model package PyLightcurve, the transit parameters shown in Supplementary Table 2, and limb-darkening coefficients (Supplementary Table 1) calculated based on the PHOENIX model, the nonlinear formula, and the stellar parameters in Supplementary Table 2.

More specifically, we fitted the white light curves with a transit model (with the planet-to-star radius ratio and the mid-transit time being the only free parameters) alongside with a model for the systematics. It is common for WFC3 exoplanets observations to be affected by two kinds of time-dependent systematics: the long-term and short-term “ramps”. The first affects each HST visit and has a linear behaviour, while the second affects each HST orbit and has an exponential behaviour. The formula we used for the systematics was the following:

At a first stage we fitted the white light curves using the formulas above and the uncertainties per pixel, as propagated through the data reduction process. However, it is common in HST/WFC3 data to have additional scatter that cannot be explained by the ramp model. For this reason, we scaled-up the uncertainties on the individual data points, in order for their median to match the standard deviation of the residuals, and repeated the fitting. From this second step of analysis, we found that the measured mid-transit times were not consistent with the expected ephemeris, which we found to be: $P=32.94007\pm 0.00003$ days and $T_{0}=2457363.2109\pm 0.0004$ BJDTDB. Supplementary Figure 1, shows the difference between the predicted and the observed transit times using the ephemeris in the literature and the one calculated in this work. We used the de-trended and time-aligned – i.e. with TTVs removed – white light curves to also refine the orbital parameters ( $a/R_{*}=81.3\pm 1.5$ and $i=89.56\pm 0.02$ deg). At a final step, we used all the new parameters (ephemeris, and orbital parameters) to perform a final fit on the white light curves (again having the planet-to-star radius ratio and the mid-transit time being the only free parameters).

Supplementary Figure 2 shows the raw white light curves, the detrended white light curves and the fitting residuals as well as a number of diagnostics, while Supplementary Table 3 presents the fitting results. From these, we can see that:

the final planet-to-star radius ratio is consistent among the eight different transits, demonstrating the stability of both the instrument and the analysis process,

on average, the white light curve residuals show an autocorrelation of 0.17, which is a low number relatively to the currently published observations of transiting exoplanets with HST (up to 0.7), indicating a good fit,

uncorrected systematics are still present in the residuals which, on average, show a scatter two times larger that the expected photon noise, and

this extra noise component is taken into account by the increased uncertainties, as the reduced $\chi^{2}$ is, on average, 1.16.

Furthermore, we fitted the spectral light curves with a transit model (with the planet-to-star radius ratio being the only free parameter) alongside with a model for the systematics that included the white light curve (divide-white method), and a wavelength-dependent, visit-long slope:

Supplementary Figures 3 to 19 show the raw spectral light curves, the detrended spectral light curves and the fitting residuals as well as a number of diagnostics, while Supplementary Table 4 presents all the fitting results and average diagnostics per spectral channel. From these, we can see that:

the spectral light curves residuals show, on average, standard deviations much closer to the photon noise and lower values of autocorrelation, proving the advantage of using the white light curve as a model compared to the ramp model, and

any extra noise component is taken into account by the scaled-up uncertainties, as the reduced $\chi^{2}$ is for all channels, on average, close to unity.

Finally, the eight spectra of K2-18 b (Supplementary Figure 20) were combined, using a weighted average, to produce the final spectrum (Table 1).

Stellar contamination

K2-18 is a moderately active M2.5V star, with a variability of 1.7% in the B band and 1.38% in the R band. Hence, to make sure that the observed water feature is not the effect of stellar contamination we fitted the observed spectrum with a model that assumes a flat planetary spectrum and contribution only from the star (M2V star as described in Rackham et al. 2018). The model that best describes our data has a spot coverage of 26% and a faculae coverage of 73%. We plot this spectrum versus the observed one in Supplementary Figure 21. In addition, we plot the spectrum produced by the spot and faculae combination reported by Rackham et al. 2018 and correspond to a 1% I-band variability, for reference. However, as Supplementary Figure 21 shows, the best-fit model cannot describe the observed water feature. From these we conclude that there is no combination of stellar properties that could introduce the observed water feature.

Atmospheric retrieval

We fitted the final planetary spectrum using our Bayesian atmospheric retrieval framework TauREx, which fully maps the correlations between the fitted atmospheric parameters through nested sampling.

We identified three solutions: a) a cloud-free atmosphere containing only H2O and H2/He, b) a cloud-free atmosphere containing H2O, H2/He and N2, and c) a cloudy atmosphere containing only H2O and H2/He. The best-fit spectra and the posterior plots are shown in Figure 3. In all cases, a statistically significant atmosphere around K2-18 b was retrieved with an ADI of 5.0, 4.7 and 4.0, respectively. The ADI is the positively defined logarithmic Bayes Factor, where the null hypothesis is a model that contains no active trace gases, Rayleigh scattering or collision induced absorption – i.e. a flat spectrum. An ADI of 5 corresponds to approximately a 3.6 $\sigma$ detection of an atmosphere. The values are too similar to distinguish between the three scenarios.

Data Availability

The data analysed in this work are available through the NASA MAST HST archive (https://archive.stsci.edu) programs 13665 and 14682. The molecular line lists used are available from the ExoMol webpage (www.exomol.com). The final and intermediate results (reduced data, extracted light curves, light curve fitting results and atmospheric fitting results) are available through the UCL-Exoplanets webpage (https://www.ucl.ac.uk/exoplanets).

Code Availability

All the software used to produced the presented results are publicly available through the UCL-Exoplanets GitHub page (https://github.com/ucl-exoplanets). More specifically, the codes used were:

TauREx (https://github.com/ucl-exoplanets/TauREx_public),

Iraclis (https://github.com/ucl-exoplanets/Iraclis), and

PyLightcurve (https://github.com/ucl-exoplanets/pylightcurve).