Water Vapor and Clouds on the Habitable-zone Sub-Neptune Exoplanet K2-18b

The recent discovery of the transiting 8.63 ± 1.35 M_⊕ exoplanet K2-18b in the habitable zone of a bright, nearby M3-dwarf provides us with an opportunity to carry out a spectroscopic study of the atmosphere of a habitable-zone planet outside our solar system (Montet et al. 2015; Benneke et al. 2017; Cloutier et al. 2017, 2019; Sarkis et al. 2018). K2-18b is an intriguing planet because its equilibrium temperature (255 ± 4 K at an albedo of A = 0.3) is potentially very close to that of the Earth (257 K). The planet's predicted temperature provides the right conditions for water vapor to condense to the liquid phase in its atmosphere. While K2-18b has a much shorter orbital period (33 days) than the Earth, its host star is also cooler and smaller than the Sun (3457 K, 0.45 R_⊙) resulting in only 2.53% of the Sun's luminosity. As a result, the total insolation received by K2-18b (${1368}_{-107}^{+114}$ W m⁻²) is very close to that received by the Earth (1361 W m⁻²), putting K2-18b firmly in the habitable zone.

Compared to habitable-zone planets around Sun-like stars, habitable-zone planets around M dwarfs offer two key advantages for atmosphere studies (Nutzman & Charbonneau 2008). The small diameter of the star results in larger transit signatures as the amplitude of transit and the atmospheric signal scale inversely with the square of the stellar radius. Furthermore, the short orbital periods for habitable-zone planets around mid M dwarfs (30–70 days) enables the observation of repeated transits within a relatively short time frame. This means that for K2-18b, we were able to observe nine transits with Wide Field Camera 3 on the Hubble Space Telescope (HST/WFC3) within a period of 3 yr.

Kepler showed that as much as 80% of M dwarfs host small planets (<3R_⊕) in the habitable zone (Dressing & Charbonneau 2013, 2015; Kopparapu 2013; Farr et al. 2015; Silburt et al. 2015), suggesting that planets like K2-18b may be common. However, most of the M dwarf planetary systems detected in the original Kepler survey are extremely faint, making spectroscopic characterization of these planets prohibitively inefficient. Fortunately, K2-18b orbits a relatively bright (K = 8.89) host star, permitting detailed characterization of its atmosphere.

Beyond its potentially temperate climate, K2-18b also occupies an interesting niche in mass–radius space. Very little is currently known about the bulk and atmospheric compositions of planets with masses between those of Earth and Neptune. These planets have no analogs in the solar system, and aside from the recent atmospheric detection for GJ 3470b (Benneke et al. 2019), most atmospheric studies resulted in nondetections due to the prevalence of high-altitude clouds (Knutson et al. 2014b; Kreidberg et al. 2014; Crossfield et al. 2017). Refined population studies of the radius distribution of sub-Neptune-sized planets have revealed a significant drop in the planet occurrence rates between near 1.5–2.0 R_⊕ (Fulton et al. 2017; Fulton & Petigura 2018). Photoevaporation (Lopez & Fortney 2013; Owen & Wu 2013; Lopez & Rice 2018) and core-powered mass loss (Schlichting 2018) are two theories that can reproduce this observed drop in the occurrence rate. According to photoevaporation models, the most highly irradiated super-Earths are expected to have primarily rocky compositions and relatively small radii, while less irradiated super-Earths are able to retain a modest (few percent in mass) primordial hydrogen-rich atmosphere that inflates their observed radii to values greater than 2.0 R_⊕. With a relatively low incident flux and a measured radius of 2.61 R_⊕, K2-18b would then be expected to host an extended hydrogen-rich atmosphere, making it a favorable target for atmospheric characterization studies using the transmission spectroscopy technique.

In this work, we present the detection of water vapor and clouds in the atmosphere of the habitable-zone exoplanet K2-18b. In Section 2 we describe our observations obtained with HST, Spitzer, and K2, as well as the techniques used to reduce the data and produce spectrophotometric lightcurves. In Section 3 we discuss the data analysis and present the best-fit white-light-curve parameters and transmission spectrum. Importantly, we also provide a refined stellar characterization. Our atmospheric modeling analysis is described in Section 4. The main results are presented in Section 5, and the possibility of liquid water clouds is discussed in Section 6.

Our team observed the transiting habitable-zone exoplanet K2-18b with Wide Field Camera 3 on the HST/WFC3 as part of two large spectral surveys of low-mass exoplanets (GO 13665 and GO 14682; PI Benneke). Building up sufficient signal-to-noise was possible for this habitable-zone planet because the shorter orbital period (33 days) enabled us to observe nine independent transits within 3 yr. The SNR of the transmission signal was further boosted by the small stellar radius, which amplifies the signal of the planet and atmosphere during transit. We complement the HST/WFC3 observations with two new Spitzer transit observations taken at 3.6 μm (Program 12081, PI Benneke) as well as previously published Spitzer (4.5 μm) and Kepler/K2 transit observations from Benneke et al. (2017; see our Table 1).

2.1. HST/WFC3 transits

2.2. Spitzer/IRAC Transits

2.3. Kepler/K2 Transits

3.1. Updated Stellar Parameters

For the correct analysis and interpretation of the K2-18b observations, it is absolutely crucial to update the stellar and planetary bulk parameters compared to the discovery and confirmation papers (Montet et al. 2015; Benneke et al. 2017). This is particularly important for K2-18b because GAIA DR2 substantially updated the parallax distance from the previously estimated 34 to 38.025 ± 0.079 pc (see also Cloutier et al. 2019). In its final consequence, this results in a substantial update to the planet radius from the previous value of R_p = 2.29 R_⊕ to R_p = 2.61 R_⊕. Not correcting the stellar radius led Tsiaras et al. (2019) to use a planetary density that was 27% too high and a surface gravity that was 18% too high. The overestimated planetary density may have led to the continued classification of K2-18b as a super-Earth, even though the radius of R_p = 2.61 R_⊕ puts K2-18b directly in the center of the intriguing sub-Neptune population (Fulton & Petigura 2018). Besides, the overestimated surface gravity unavoidably resulted in a systematic offset in the atmospheric scale height for each model in the atmospheric retrieval.

To revise the stellar parameters, we first calculate the distance modulus from the Gaia DR2 stellar parallax to be μ = 2.900 ± 0.005, where we add 30 μas systematic offset in the measured parallax (Cloutier et al. 2019). By propagating the uncertainties in K2-18's K-band magnitude (K = 8.899 ± 0.019), we then find an absolute K-band magnitude of M_K = 5.999 ± 0.020. Using this absolute K-band magnitude, we infer the stellar radius of K2-18b from the empirical M star correlation (Mann et al. 2015) and the stellar mass from the empirical M dwarf mass–luminosity relation by Benedict et al. (2016; see our Table 2). We choose this approach for the stellar radius because it does not rely on the empirical mass–radius relationship, which we agreed with Dr. Cloutier to be the preferred approach. Either way the stellar radii presented here and the one from Cloutier et al. (2019) are consistent at 1σ. However, both estimates present a substantial update to the old pre-GAIA stellar radius estimate, which led Tsiaras et al. (2019) to use an overestimated planetary density and surface gravity for their analyses. Finally, we use the effective stellar temperature from Benneke et al. (2017) to compute the stellar luminosity as well as the insolation of the planet (Table 2).

Table 2.  Parameters of the Host Star K2-18, the Planet K2-18b, and Its Atmosphere

Parameter	Unit	Value	Comment
Host star:
M⋆	M⊙	0.4951 ± 0.0043	Cloutier et al. (2019)
R⋆	R⊙	0.4445 ± 0.0148	this paper
${T}_{\star ,\mathrm{eff}}$	M⊙	3457 ± 39	Benneke et al. (2017)
L⋆	L⊙	0.0253 ± 0.0021	this paper
Distance	pc	38.025 ± 0.079	Cloutier et al. (2019)
Orbit:
RP/R⋆	⋯	${0.05375}_{-0.00012}^{+0.00012}$	this paper
T0	BJDUTC	${2457725.551189}_{-0.000068}^{+0.000067}$	this paper
P	days	${32.940045}_{-0.000010}^{+0.000010}$	this paper
ap	au	${0.15910}_{-0.00047}^{+0.00046}$	this paper
b	⋯	${0.641}_{-0.015}^{+0.015}$	this paper
Planet:
Rp	R⊕	2.610 ± 0.087	this paper
Mp	M⊕	8.63 ± 1.35	Cloutier et al. (2019)
${\rho }_{p}$	${\rm{g}}\,{\mathrm{cm}}^{-3}$	${2.67}_{-0.47}^{+0.52}$	this paper
gp	m s−2	${12.43}_{-2.07}^{+2.17}$	at 10 mbar, this paper
Sp	S⊕	${1.005}_{-0.079}^{+0.084}$	this paper, ${1368}_{-107}^{+114}\,{\rm{W}}\,{{\rm{m}}}^{-2}$
Teq	K	254.9 ± 3.9	this paper, for AB = 0.3
Atmosphere:
${\mathrm{logH}}_{2}{\rm{O}}$	⋯	$-{2.08}_{-1.39}^{+1.03}$	0.0337%–8.88%
$\mathrm{log}({{\rm{H}}}_{2}{\rm{O}}/{{\rm{H}}}_{2}{{\rm{O}}}_{\mathrm{solar}})$	⋯	${1.18}_{-1.39}^{+1.03}$	0.615–162 × solar
$\mathrm{log}{p}_{\mathrm{cloud} \mbox{-} \mathrm{top}}$	mbar	${1.37}_{-0.48}^{+0.77}$	7.74–139 mbar
μatm	⋯	${2.42}_{-0.12}^{+1.27}$

Download table as:  ASCIITypeset image

3.2. Transit White-light-curve Fitting

We carry out a global analysis of our HST/WFC3, K2, and Spitzer transit light curves within the ExoTEP analysis framework (Benneke et al. 2017, 2019). ExoTEP jointly fits the transit and systematics models of all transits along with photometric noise parameters using a Markov Chain Monte Carlo (MCMC) method. The main astrophysical outputs of the white-light-curve analysis are the global transit parameters (a/R_⋆, b), the transit ephemeris (T₀, P), and the transit depths in the K2, HST/WFC3, and Spitzer/IRAC bandpasses (Table 2). Prior to the global MCMC fit, we analyze each transit light curve individually and then initialize the corresponding systematics parameters in the global fit at their best-fit values to ease convergence. We find that the planet-to-star radius ratio estimates in the same band are consistent within their statistical uncertainties across all epochs. The white-light-curve fits to the HST/WFC3 and Spitzer transits are depicted in Figures 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.1. HST/WFC3 Instrument Model

We correct for systematic trends in the uncorrected WFC3 transit light curves by simultaneously fitting an analytical model-ramp function along with the astrophysical transit model. Following previous studies (e.g., Berta et al. 2012; Deming et al. 2013; Kreidberg et al. 2014), we account for the possible presence of both visit-long slopes and orbit-long exponential ramps using the parametric instrumental systematics model:

$$\begin{eqnarray}&&{S}_{\mathrm{WFC}3}(t)=(c\,S(t)+v\,{t}_{v})\times (1-\exp (-a\,{t}_{\mathrm{orb}}-b)).\end{eqnarray} \tag{ 1 }$$

Here, c is a normalization constant, v is the visit-long linear slope, a and b describe the rate and amplitude of the orbit-long exponential slope, S(t) is set to 1 for forward scans and s for backward scans to account for the systematic flux offset between the scan directions, and t_v and t_orb are the time in hours since the start of the visit and the start of the observations within the current orbit. Following standard procedure, we discard the first HST orbit of each visit, as it systematically exhibits a stronger ramp than the three subsequent orbits. We also remove the first forward and backward scan exposures of each orbit.

3.2.2. Spitzer/IRAC Instrument Model

We account for the presence of systematic variations in the Spitzer time series due to subpixel inhomogeneities in the detector response using the pixel-level decorrelation (PLD) model (Deming et al. 2015; Benneke et al. 2017). For each Spitzer visit, the systematics model

$$\begin{eqnarray}&&{S}_{{Spitzer}}({t}_{i})=1+A\,{e}^{-{t}_{i}/\tau }+m\,t+\displaystyle \frac{{\displaystyle \sum }_{k=1}^{9}{w}_{k}{D}_{k}({t}_{i})}{{\displaystyle \sum }_{k=1}^{9}{D}_{k}({t}_{i})}\end{eqnarray} \tag{ 2 }$$

includes both a linear-exponential ramp in time ($A\,{e}^{-{t}_{i}/\tau }+{mt}$) and the PLD term. In this analytical model, D_k(t_i) are the registered counts in each of the 3 × 3 pixels covering the central region of the point-spread function. The coefficients w_k are time-independent PLD weights. The 10 parameters in this model are fitted along with the transit model for all Spitzer data sets.

3.2.3. Kepler/K2 Instrument Model

We use the detrended K2 light curves of K2-18b from Benneke et al. (2017) and additionally fit a linear trend with time to allow for a residual slope in each visit superimposed with the transit signal.

3.2.4. Transit Model and MCMC Analysis

Our astrophysical transit light-curve model f(t_i) is computed using the Batman module (Kreidberg 2015). In the joint white-light-curve fit, four distinct transit depths are fitted for the four bandpasses of K2 and HST/WFC3, Spitzer/IRAC Channel 1, and Spitzer/IRAC Channel 2. The orbital parameters are assumed to be consistent for all 12 transit light curves. We also fit limb-darkening coefficients using a quadratic law because the eight HST/WFC3 combined provide full time coverage of the transit light curve including ingress and egress (Figure 1). The cadence of the observations is accounted for by integrating the model over time within each exposure. Finally, the total log-likelihood for a single set of parameters is calculated using

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}{ \mathcal L } & = & \displaystyle \sum _{V=1}^{N}-{n}_{V}\mathrm{ln}{\sigma }_{V}-\displaystyle \frac{{n}_{V}}{2}\mathrm{ln}2\pi \\ & & -\,\displaystyle \sum _{i=1}^{{n}_{V}}\displaystyle \frac{{\left[{d}_{V}({t}_{i})-{S}_{V}({t}_{i})\times {f}_{V}({t}_{i})\right]}^{2}}{{\sigma }_{V}^{2}},\end{array}\end{eqnarray} \tag{ 3 }$$

where N is the number of visits, n_V is the number of data points d_V (t_i) in visit V, and σ_V is a photometric noise parameter associated with each visit, simultaneously fitted to account for the possibility of variations in the scatter between independent visits. Each visit has a different systematics model S_V (t_i) with additional free parameters specific to the instrument that performed the observations, as described in the previous sections. The log-likelihood is used for finding best-fit parameters and to obtain the joint posterior distribution of all astrophysical and instrumental model parameters using the emcee package (Foreman-Mackey et al. 2013). All parameters have flat priors and initial values for the joint MCMC fit are set to the best-fitting parameters from the fit to the individual transit light curves. The median and 1σ values of the orbital parameters and transit depths in the K2 and Spitzer bandpasses are quoted in Tables 2 and 3, respectively. Within emcee, we use four walkers per free parameter and deem the chains converged once their lengths exceed at least 50 times the autocorrelation time computed using the new "Farad" method defined in the emcee 3.0.0 documentation. The latter is a more conservative approach compared to autocorrelation time defined in Goodman & Weare (2010).

Table 3.  Optical/IR Transmission Spectrum of K2-18b

Instrument	Wavelength	Depth	+1σ	−1σ
	(μm)	(ppm)	(ppm)	(ppm)
K2	0.43–0.89	2872.9	51.9	53.1
HST/WFC3	1.12–1.15	2916.2	27.9	26.8
	1.15–1.18	2900.0	24.7	25.7
	1.18–1.21	2872.5	24.5	24.6
	1.21–1.24	2909.7	23.0	23.0
	1.24–1.27	2909.2	24.5	24.0
	1.27–1.30	2898.1	22.4	23.1
	1.30–1.33	2896.9	22.8	22.7
	1.33–1.36	2969.2	23.2	22.8
	1.36–1.39	2954.5	22.8	22.3
	1.39–1.42	2953.4	23.4	24.4
	1.42–1.45	2970.9	23.0	22.2
	1.45–1.48	2907.1	22.8	22.8
	1.48–1.51	2934.7	23.3	22.5
	1.51–1.54	2890.9	22.2	22.9
	1.54–1.57	2868.6	24.7	25.0
	1.57–1.60	2895.7	21.9	22.1
	1.60–1.63	2878.2	23.5	24.5
Spitzer/IRAC Ch1	3.15–3.94	2849.3	91.0	89.3
Spitzer/IRAC Ch2	3.96–5.02	2882.2	93.6	91.6

Download table as:  ASCIITypeset image

3.3. HST/WFC3 Transit Spectroscopy

Following standard practice we use the results of the white-light-curve fitting to pre-correct the systematics in each spectroscopic light curve. This is possible because the ramp-like systematic variations in the HST/WFC3 data are to first-order independent of wavelength. We test two methods for the correction: either dividing each of the spectroscopic time series by its corresponding best-fit systematics model from the white-light-curve analysis, or dividing it by the ratio of the white light curve to its best-fitting transit model. We find no significant difference between the two approaches for pre-correcting the light curves in terms of the derived parameters.

In each spectroscopic channel, we then perform a joint MCMC fit to the eight pre-corrected spectroscopic light curves from the eight HST/WFC3 transits. We keep the parameters a/R_⋆, b, T₀, and P fixed to the best-fitting values from the white-light-curve fit and fit a simplified systematics model that detrends linearly against the x position of the trapezoidal spectral trace on the detector. We obtain consistent transmission spectra when fitting either a single linear limb-darkening coefficient at each wavelength or two coefficients for a quadratic limb-darkening profile. A consistent transmission spectrum is also obtained when using pre-computed limb-darkening coefficients from stellar atmosphere models (see Lothringer et al. 2018 and Benneke et al. 2019 for details). The resulting WFC3 transit depths for the latter are quoted in Table 3 and depicted in Figure 3 together with the K2 and Spitzer transit depths. A typical spectral light-curve fit for one wavelength bin is depicted in Figure 1.

Zoom In Zoom Out Reset image size

3.4. Stellar Activity

K2-18 is a moderately active M dwarf with a measured rotation period of approximately 39 days and corresponding photometric variability of 8 mmag in the Kepler bandpass (Cloutier et al. 2017), 9 mmag in B band, and 7 mmag in R band (Sarkis et al. 2018). We see no evidence for spot occultations in any of our transit light curves, and when comparing the measured transit depth across multiple visits in the same bandpass we obtain consistent values, implying that the observed transit depths are relatively insensitive to variable spot coverage on the star.

Still, to further rule out a stellar origin, we attempt to fit the transit spectrum with a stellar contamination model. Following Rackham et al. (2018) we consider the hypothesis that the observed transit depth variations are due to inhomogeneities on the stellar surface. The observed wavelength-dependent transit depths ${D}_{\lambda ,\mathrm{obs}}$ are modeled as:

$$\begin{eqnarray}&&{D}_{\lambda ,\mathrm{obs}}=\displaystyle \frac{{D}_{\lambda }}{1-{f}_{\mathrm{spot}}\left(1-\tfrac{{F}_{\lambda ,\mathrm{spot}}}{{F}_{\lambda ,\mathrm{phot}}}\right)-{f}_{\mathrm{fac}}\left(1-\tfrac{{F}_{\lambda ,\mathrm{fac}}}{{F}_{\lambda ,\mathrm{phot}}}\right)},\end{eqnarray} \tag{ 4 }$$

where D_λ is the geometric radius ratio ${\left({R}_{p}/{R}_{\star }\right)}^{2}$, assumed here to be constant with wavelength (=D), f_spot and f_fac are the spot and faculae covering fractions, and the F_λ refer to bandpass-integrated PHOENIX model spectra of the photosphere, spots, and faculae (Hauschildt et al. 1999). Using K2-18b stellar parameters (Table 2), we explore a wide range of possible stellar contamination spectra for M2–M3 stars with spot and faculae covering fractions consistent with 1% I-band variability. Based on Rackham et al. (2018), we then explore scenarios with giant spots and solar-like spots as well as with and without faculae, resulting in spot fractions of up to 20%. However, none of those models can explain the amplitude of the observed transit spectrum. The amplitude of resulting transit depth variations due to the stellar inhomogeneities are up to 3–8 parts-per-million (ppm), about an order of magnitude smaller than the observed 80–90 ppm transit depth variations in the observed WFC3 spectrum (Figure 3). Finally, we also explore the most extreme scenarios where the spot and faculae covering fractions can be as high as 100%, but even those stellar inhomogeneity models fail to explain the amplitude of the observed transit depth variation. They deliver an absolute maximum of 20 ppm at 1.4 μm, which still only corresponds to less than a quarter of the transit depth variation in the observations. We conclude that stellar inhomogeneities and activity cannot explain the measured transmission spectrum.

We compute quantitative constraints on the atmosphere of K2-18b using the SCARLET atmospheric retrieval framework (Benneke & Seager 2012, 2013; Knutson et al. 2014a; Kreidberg et al. 2014; Benneke 2015; Benneke et al. 2019). To be as independent of model assumption as possible, we employ the "free retrieval" mode, which parameterizes the mole fractions of the molecular gases, the pressure of the cloud deck, and the atmospheric temperature as free fitting parameters. SCARLET then determines their posterior constraints by combining the atmospheric forward model with a Bayesian MCMC analysis.

To evaluate the likelihood for a particular set of parameters, the atmospheric forward model first computes a model atmosphere in hydrostatic equilibrium, then determines the opacities of molecules at each layer, and finally computes the transmission spectrum. Beyond H₂/He, our model allows for H₂O, CH₄, CO, CO₂, NH₃, HCN, and N₂ with a log-uniform prior for mixing ratios between 10⁻¹⁰ and 1. We find that only H₂O is required by the data and that including the other molecules has virtually no impact on the best fit to the data. Following Benneke & Seager (2012, 2013), we also include a cloud deck at a freely parameterized cloud top pressure with a log-uniform prior between 0.1 mbar and 10 bar. The cloud deck is assumed to be opaque to grazing light beams below the cloud top pressure as would occur for large droplets. We also explored a more complex three-parameter Mie-scattering cloud description as introduced in Benneke et al. (2019); however, we find no significant improvement in the fit to the observed transmission spectrum compared to gray clouds. Our atmospheric temperature is parameterized using a single free parameter for the mid-atmosphere probed by the observations because low-resolution transmission spectra are largely insensitive to the exact vertical temperature structure (Benneke & Seager 2012). We also considered a five-parameter analytic model (Parmentier & Guillot 2014; Benneke et al. 2019), but no gradients in the temperature–pressure profile can be constrained.

High-resolution synthetic transmission spectra are computed using line-by-line radiative transfer, which are integrated over the appropriate instrument response functions to obtain synthetic observations to be compared to the observations. Sufficient wavelength resolution in the synthetic spectra is ensured by repeatedly verifying that the likelihood for a given model is not significantly affected by the finite wavelength resolution (Δχ² < 0.001). Reference models are computed at λ/Δλ = 250,000. Our final retrieval run computes all atmospheric models at λ/Δλ = 30,000.

We perform all retrieval analyses with 100 walkers using log-uniform priors on all parameters. Following best practice for emcee (Foreman-Mackey et al. 2013), our formal convergence test is based on the integrated autocorrelation time (Goodman & Weare 2010), but we run the chains well beyond formal convergence to obtain smooth posterior distribution contours. We double-check that the number of iterations of our 100 walker chain exceeds at least 50 times the autocorrelation time computed based on Goodman & Weare (2010) as well as the autocorrelation time computed computed based on the more conservative "Farad" method in the emcee 3.0.0 documentation.

Our transmission spectrum of K2-18b reveals a slightly attenuated but statistically significant water absorption feature at 1.4 μm in our HST/WFC3 data (Figure 3). The water absorption is detected in multiple neighboring spectroscopic channels covering the 1.4 μm water band and protrudes over an otherwise relatively flat visible to near-IR transmission spectrum. Quantitatively, retrieval models that include molecular absorption by water are favored by the Bayesian factor at 459:1 (3.93σ). We evaluate the Bayes factor by comparing the Bayes evidence for models with and without water opacity (Benneke & Seager 2013). We also evaluate the AIC of our model in comparison to the simple flat line model with R_P/R_⋆ as the only free parameter and find a strong preference of ΔAIC = 18.1 in favor of the model with water absorption.

Our retrieval modeling also shows that the data are best matched by a hydrogen-dominated atmosphere with water vapor absorbing above clouds in the mid-atmosphere (Figure 4). Water abundances between 0.033% and 8.9% are most consistent with the data at 1σ, corresponding to 0.6 times and 162 times the value for a solar abundance atmosphere, respectively. The clouds become optically thick below the 7.7–139 mbar level. The water abundance and cloud top pressure are strongly correlated because models with similar water vapor column density above the cloud deck result in similar transmission spectra (see also Benneke & Seager 2013). For the lowest water abundances, pressure levels up to a few hundred millibars are probed by our observations because K2-18b's radius is relatively small and the low temperature results in a relatively small atmospheric scale height. Both effects reduce the path lengths for starlight grazing through the atmosphere compared to most previously studied exoplanets.

Zoom In Zoom Out Reset image size

Cloud-free models with water mixing ratios greater than several hundred times solar are substantially disfavored by our combined data set indicating that K2-18b hosts an atmosphere rich in hydrogen and helium. This is consistent with our new measurement of the planet density (Section 3.1) which supports the presence of a thick hydrogen-rich gas envelope (Lopez & Fortney 2014). Quantitatively, the retrieval model with clouds in the mid-atmosphere is favored by at least 8.3:1 (2.6σ) in the Bayesian evidence over retrieval models without clouds. The significance rises further if extreme subsolar O/H ratios for K2-18b are not considered. We conclude that it is crucial to report the water abundance constraints from retrieval models that consider clouds. Cloud-free retrieval models would misleadingly favor extremely high water mole fractions and narrower constraints.

Finally, while no molecular species other than H₂O is directly detected in our observations of K2-18b, we can still infer useful upper limits on their mole fractions. The two dominant effects of additional molecules are the introduction of additional wavelength-dependent opacities and, at high mole fraction, the increase of the overall atmospheric mean molecular weight (Benneke & Seager 2012). The latter, in turn, dampens the amplitudes of all molecular absorption features in the transmission spectrum. The upper limit on the spectrally inactive gas N₂, for example, is 10.9% at 2σ (97.5%) because above that level the atmospheric mean molecular mass would be sufficiently increased to not allow for the amplitude of the detected 1.4 μm H₂O absorption feature. The other molecular gases (CO, CO₂, NH₃, and CH₄) are similarly constrained to 7.45%, 2.4%, 13.5%, and 0.248% at 2σ, respectively. Molecular opacities also played an important role for these constraints because CO, CO₂, NH₃, and CH₄ can result in substantial opacities within the HST/WFC3 and/or Spitzer bandpasses.

The simultaneous detection of water vapor and clouds in the mid-atmosphere of K2-18b is particularly intriguing because K2-18b's stellar insolation and atmospheric temperature are sufficiently low for water vapor to condense. Water clouds are therefore a very plausible explanation for the detected clouds. To demonstrate the plausibility of liquid droplet formation, we explore self-consistent temperature structures for atmospheric compositions consistent with the observations (Figure 5). We computed a set of fiducial temperature-pressure profiles using the implementation described in Benneke et al. (2019) and Morley et al. (2013). Both models iteratively solve the radiative–convective heat transport. The models predict the onset of liquid and icy condensation between approximatively 10 and 1000 mbar at locations where the water vapor is supersaturated. The fact that the best-fitting retrieval models are obtained for cloud top pressures in the same range further supports the scenario of liquid water clouds (Figure 4).

Zoom In Zoom Out Reset image size

The water vapor mixing ratio retrieved from our spectroscopic observations should be regarded as the mole fraction of water vapor in the gaseous atmosphere above the clouds. Because water condenses into clouds, this presents a lower limit on the overall water fraction in the atmosphere K2-18b. The water abundance in the upper atmosphere is set by a complex interplay between upward transport through vertical mixing and homogeneous/heterogeneous condensation and coagulation processes. Detailed modeling of these cloud formation processes will be presented in a follow-up study.

The detection of water vapor makes K2-18b a key target for more detailed follow-up studies with the upcoming James Webb Space Telescope (JWST). Unlike any other temperate and low-mass planet, we now know that K2-18b shows evidence for atmospheric water vapor and is amenable to characterization via transmission spectroscopy. Our result suggests K2-18b orbits far enough away from its moderately active host star to avoid high levels of haze production, while also remaining cool enough that cloud species like KCl and ZnS that may exist in warmer atmospheres do not form. This may indicate that cooler planets like K2-18b occupy a sought after niche in the low-mass planet parameter space that is accessible to spectroscopic characterization with JWST. JWST's wavelength coverage will extend from 0.55 μm to the thermal infrared, where many other molecular species like CH₄, CO, CO₂, and NH₃ can be probed directly. The higher precision and spectral resolution obtained from repeated JWST transit observations will allow us to better constrain K2-18b's atmospheric composition and cloud properties, and potentially even look for biomarkers in the gas envelope of a habitable-zone exoplanet (e.g., Seager et al. 2013).

This work is based on observations with the NASA/ESA HST, obtained at the Space Telescope Science Institute (STScI) operated by AURA, Inc. We received support for the analysis by NASA through grants under the HST-GO-13665 and HST-GO-14682 programs (PI Benneke). This work is also based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA (PIs Benneke and Werner). B.B. further acknowledges financial supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada and the Fond de Recherche Québécois-Nature et Technologie (FRQNT; Québec). D.D. acknowledges support provided by NASA through Hubble Fellowship grant HST-HF2-51372.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.