A DETECTION OF WATER IN THE TRANSMISSION SPECTRUM OF THE HOT JUPITER WASP-12b AND IMPLICATIONS FOR ITS ATMOSPHERIC COMPOSITION

We obtained time series spectroscopy during six transits of WASP-12b between UT 2014 January 1 and March 4 using the Wide Field Camera 3 (WFC3) IR detector on the Hubble Space Telescope (HST) as part of HST GO Program 13467. Three of the transit observations used the G102 grism, which provides spectral coverage from 0.82 to 1.12 μm, and the other three used the G141 grism, which spans the range 1.12–1.65 μm. The G141 grism has been widely used for exoplanet transit spectroscopy, but these observations are the first to use the G102 grism for this purpose. The additional wavelength coverage from the G102 grism provides access to features from more molecules and absorption bands, giving us greater leverage in constraining the atmospheric composition of the planet. Each transit observation (called a visit) consisted of five consecutive, 96-minute HST orbits. WASP-12 was visible during approximately 45 minutes per orbit and occulted by the Earth for the remainder of the time. We took a direct image of the target with the F126N narrow-band filter at the beginning of each orbit for wavelength calibration. Example staring mode and spatial scan images are shown in Figure 1.

Zoom In Zoom Out Reset image size

The spectroscopic data were obtained in spatial scan mode with the 256 × 256 subarray, using the SPARS10, NSAMP = 16 readout pattern, which has an exposure time of 103.1 s. The scan rates were 0.04 and 0.05 arcsec s⁻¹ for the G102 and G141 grisms, respectively. The spectra extend roughly 40 pixels in the spatial direction, with peak per-pixel counts below 25,000 electrons for both grisms. We alternated between forward and reverse scanning along the detector to decrease instrumental overhead time. This setup yielded 19 exposures per orbit and a duty cycle of 74%.

In addition to the spatial scan data, we also obtained 10 staring mode spectra in each grism during the first orbit of the first visit. WASP-12 is a triple star system: WASP-12 A hosts the planet, and WASP-12 BC is an M-dwarf binary separated from WASP-12 A by about 1'' (Bergfors et al. 2013; Bechter et al. 2014). We use these staring mode data to resolve the spectrum of WASP-12 A from WASP-12 BC, enabling us to correct for dilution to the planet's transit light curve due to the binary. The detector orientation was set to 1787 in order to spatially separate the spectra. At this orientation, the spectra are separated by 7.9 pixels (compared to the FWHM of 1.1 pixel at 1.4 μm). We describe the dilution correction in detail in Section 3.3.

We also obtained photometric monitoring of the WASP-12 system to search for stellar activity. We acquired 428 out-of-transit R-band images over the 2011–2012, 2012–2013, and 2013–2014 observing seasons with Tennessee State University's Celestron 14-inch automated imaging telescope. See Table 1 for a summary of the photometric monitoring campaign. We tested for variability from star spots by calculating differential magnitudes for WASP-12 relative to 17 comparison stars (shown in Figure 2). The standard deviation of the differential magnitudes is 0.3%, which is comparable to the photon noise for the data. The brightness of the star increases by approximately 0.001 mag per year. There are no significant periodicities between 1 and 200 days. To calculate the impact of star spots on our transmission spectrum, we took 0.3% as an upper limit for the variability of WASP-12. Variability of this amplitude could be produced by star spots 300 K cooler than the star's effective temperature (6300 K; Hebb et al. 2009), covering 3% of the photosphere. Based on the formalism outlined in Berta et al. (2011) and Désert et al. (2011), we calculated that the maximum variation in transit depth due to spots with these properties is of order 10⁻⁵, which is below the precision of our transit depth measurements.

Zoom In Zoom Out Reset image size

We reduced the HST/WFC3 data using the custom pipeline described in Kreidberg et al. (2014b). Example extracted spectra are shown in Figure 3. We bin the spectra into 24- and 15-pixel wide channels (for G102 and G141, respectively) to obtain a total of 13 spectrophotometric light curves at resolution R ≡ λ/Δλ = 15–25. We also sum the spectra over the full wavelength range to create broadband ("white") light curves for each grism. The uncertainty in the flux per exposure is determined by adding in quadrature the photon noise, the read noise (22 electrons per differenced image), and the error in the estimate of the background, which we determine from the median absolute deviation of the flux values for the background pixels.

Zoom In Zoom Out Reset image size

We also extract spectra from the staring mode data to determine the wavelength-dependent flux ratio of WASP-12 BC to WASP-12 A. The spectra are spatially separated on the detector, so to extract them we fit a double Gaussian model to each column of the final detector read. We obtain final estimates and uncertainties for the spectra of WASP-12 A and WASP-12 BC by taking the mean and standard deviation of the spectra from all 10 exposures. These extracted spectra are each assigned a wavelength solution using the centroids of the stars in the direct image.

We fit the broadband light curves with a transit model (Mandel & Agol 2002) and an analytic function to correct for instrument systematics. Transit light curve observations with HST/WFC3 have several well-documented systematic trends in flux with time, including visit-long slopes and orbit-long exponential ramps (Berta et al. 2012; Swain et al. 2013; Kreidberg et al. 2014b; Stevenson et al. 2014a). The first orbit of a visit has a larger amplitude ramp than subsequent orbits, so we follow established practice and do not use this initial orbit in the light curve fits. We also discard the first exposure from the remaining orbits, which improves the fit quality. We fit the remaining data (216 exposures for each grism) with a systematics model based on the model-ramp parameterization of Berta et al. (2012). The model M(t) has the form:

$$\begin{eqnarray}M(t) & = & T(t)\times (c\;S(t)+v\;{t}_{{\rm{v}}})\\ & & \times (1-\mathrm{exp}(-a\;{t}_{\mathrm{orb}}-b-D(t))).\end{eqnarray} \tag{ 1 }$$

In our fits, T(t) is the transit model, or relative stellar flux as a function of time t (in BJD_TDB). The free parameters for the transit model are the planet-to-star radius ratio k, a linear limb darkening parameter u, the ratio of semimajor axis to stellar radius a/R_s, the orbital inclination i, and the time of mid-transit T₀. The light curves have poor coverage of the planet's ingress, so we put priors on a/R_s and i to enable the measurement of T₀. We use Gaussian priors with mean and standard deviation 2.91 ± 0.02 and 8056 ± 003 (for a/R_s and i, respectively), based on estimates of those parameters from Stevenson et al. (2014c). We use an orbital period P = 1.091424 days and assume a circular orbit (Campo et al. 2011). The data from the two grisms are fit separately, but for each grism the three transits are fit simultaneously. We fit unique values of k and T₀ to each transit, but tie the values for u, a/R_s, and i over all the transits.

We fit the instrument systematics with a constant normalization term c, a scaling factor S(t), a visit-long linear slope v, and an orbit-long exponential ramp with rate constant a, amplitude b, and delay D(t). The timescale t_v corresponds to time relative to the expected transit midpoint for each visit, and t_orb is time since the first exposure in an orbit (both in BJD_TDB). The scaling factor S(t) is equal to 1 for exposures with forward spatial scanning and s for reverse scanning; this accounts for a small offset in normalization between the scan directions caused by the upstream-downstream effect of the detector readout (McCullough & MacKenty 2012). The function D(t) is equal to d for times t during the first fitted orbit and 0 elsewhere. A negative value for d implies that the ramp amplitude is larger in the first orbit than in subsequent orbits. The parameters u, a, b, and d were constrained to the same value for all the transits, whereas k, s, c, and v were allowed to vary between transits.

There are a total of 21 free parameters in the fits to each grism's broadband light curve. We estimated the parameters and their uncertainties with a Markov chain Monte Carlo (MCMC) fit to the data, using the emcee package for Python (Foreman-Mackey et al. 2013). For the best-fit light curves, we calculated the Durbin-Watson statistic to test for time-correlated noise. The values were 1.77 and 2.17 for the G102 and G141 white light curves, indicating that the residuals are uncorrelated at the 1% significance level. The reduced chi-squared values (χ²_ν) for the best-fit light curves are 1.45 and 1.34, and the residuals are 164 and 115 ppm (37 and 25% larger than the predicted photon+read noise) for the G102 and G141 data, respectively.

To obtain conservative estimates of the errors for the white light curve transit parameters, we redid the MCMC analyses with the per-exposure uncertainties scaled by a constant factor (1.20 and 1.16) chosen to yield ${\chi }_{\nu }^{2}=1$. We report the transit times from the white light curve fits in Table 2. These values extend the baseline of precise transit times by two years from the Stevenson et al. (2014c) measurements and will aid in testing the possibility of perihelion precession for this system (first studied by Campo et al. 2011). The white light curves are consistent with the priors on a/R_s and i but do not yield improved values for those parameters due to the poor phase coverage of the transit ingress.

We fit the spectroscopic light curves with a similar model as we use for the broadband data. The only differences are that we hold a/R_s and i fixed to the prior mean values and fix the mid-transit times to the white light curve best fit values. We also solve for a constant rescaling parameter for the photometric uncertainties in each spectroscopic channel to ensure that the reduced χ² for the light curve fits is unity. We chose to fix a/R_s, i, and T_c because they impact the mean spectroscopic transit depth only, not the relative depths. Changes in the mean transit depth do not significantly affect our retrieval results, as the planet-to-star radius ratio is included as a free parameter in the atmospheric retrieval.

We fit each of the spectroscopic channel light curves independently. We achieve nearly photon-limited precision: the median rescaling factor for the photometric uncertainties is 1.1. The light curves do not exhibit statistically significant time-correlated noise, based on a Durbin-Watson test at the 1% significance level. An example spectroscopic light curve fit for the 1.320–1.389 μm channel from the G141 grism is shown in Figure 4. We show a pairs plot of the fit parameters for the same channel in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We tested an analytic model for the systematics that included a quadratic term for the visit-long trend. Previous analyses of WFC3 data have suggested that a quadratic model produces better light curve fits (Stevenson et al. 2014a). However, for this data set we find that the quadratic model is disfavored according to the Bayesian Information Criterion (BIC), with typical Δ BIC values of 10 compared to the linear model. In any case, the main consequence of adding a quadratic term is to shift the transmission spectrum up or down. This affects our estimate of the planetary radius, but the effect is small relative to the error introduced by uncertainties in the stellar radius.

We also explored modeling the instrument systematics with the divide-white technique, which has been applied successfully to several other data sets (Knutson et al. 2014; Kreidberg et al. 2014b; Stevenson et al. 2014c). This method assumes the systematics are independent of wavelength. For the WASP-12 data, however, this assumption is not appropriate. WFC3 instrument systematics are known to depend on detector illumination (Berta et al. 2012; Swain et al. 2013), and in our data the mean pixel fluence varies by 30% between the spectroscopic channels (see Figure 3).

To illustrate the dependence of the systematics on the illumination level, we show in Figure 6 the systematics decorrelation parameters from the analytic model as a function of light curve normalization c. The value c represents the baseline flux level in each spectroscopic light curve. This value is directly proportional to the mean detector illumination in the channel. We note several qualitative trends in the decorrelation parameters: with increasing illumination, the visit-long slope decreases, the delay term increases toward zero, and the ramp rate increases.

Zoom In Zoom Out Reset image size

The upshot of these trends is that there are residual systematics for the divide-white light curve fits that are correlated with detector illumination. We therefore only report the results for the analytic model.

3.2.1. Limb Darkening Models

The transit light curves are affected by dilution from WASP-12 A's two companion stars, WASP-12 BC, and from the planet's nightside emission. Following Stevenson et al. (2014c), we calculated a corrected transit depth δ'for each spectroscopic channel using the formula

$$\begin{eqnarray}&&\delta ^{\prime} =\delta (1+{\alpha }_{\mathrm{Comp}}+{\alpha }_{p})\end{eqnarray} \tag{ 2 }$$

where δ is the measured transit depth, α_Comp is the ratio of flux from WASP-12 BC to WASP-12 A, and α_p is the ratio of the planet's nightside emission to the flux of WASP-12 A. We calculated α_p using the same model as Stevenson et al. (2014c). We determined α_Comp empirically using the spectra extracted from the staring mode observations. Each visit is assigned a unique α_Comp value to account for differences in the orientation of the spectra on the detector. We plot the dilution factor as a function of wavelength in Figure 7 (calculated by interpolating the companion spectrum to the same wavelength scale as the WASP-12 A spectrum). The measured dilution is roughly 10% larger than the values used by Stevenson et al. (2014c); this difference could result from systematic uncertainty in the aperture photometry used to scale the dilution (Stevenson et al. 2014c).

Zoom In Zoom Out Reset image size

To obtain final transit depths, we take the weighted average of the corrected transit depths from each visit. Uncertainties in the dilution factor corrections are propagated through to the final transit depth uncertainties. We report the corrected transit depth measurements in Table 3, and we show the transmission spectrum in Figure 8. The corrected transit depths are consistent between the visits, indicating that the estimated uncertainties are appropriate, and confirming that stellar activity does not significantly influence the measured spectrum.

Zoom In Zoom Out Reset image size

Table 3.  Transit Fit Parameters

Bandpassa	(Rp/R⋆)2b	Limb-	rms
(μm)	(%)	darkening u	(ppm)
0.838–0.896	1.4441 ± 0.0069	0.27 ± 0.02	349
0.896–0.954	1.4422 ± 0.0055	0.28 ± 0.01	300
0.954–1.012	1.4402 ± 0.0052	0.27 ± 0.01	296
1.012–1.070	1.4428 ± 0.0051	0.26 ± 0.01	301
1.070–1.129	1.4391 ± 0.0053	0.25 ± 0.01	299
1.112–1.182	1.4386 ± 0.0047	0.26 ± 0.01	267
1.182–1.251	1.4365 ± 0.0045	0.26 ± 0.01	267
1.251–1.320	1.4327 ± 0.0041	0.21 ± 0.01	253
1.320–1.389	1.4582 ± 0.0040	0.22 ± 0.01	242
1.389–1.458	1.4600 ± 0.0043	0.18 ± 0.01	255
1.458–1.527	1.4530 ± 0.0045	0.20 ± 0.01	270
1.527–1.597	1.4475 ± 0.0058	0.16 ± 0.02	332
1.597–1.666	1.4332 ± 0.0055	0.16 ± 0.02	321

Notes.

^aThe measurements between 0.838 and 1.129 μm are from the G102 grism; those from 1.112 to 1.666 μm are from the G141 grism. ^bTransit depths corrected for dilution from companion stars and planet nightside emission. The light curve fits had a/R_s and i fixed to 291 and 8056, respectively.

Download table as:  ASCIITypeset image

The FREE parameterization retrieves molecular abundances, cloud and haze properties, and the altitude-independent scale height temperature.

Our model explored the dominant molecular opacities expected for a solar composition gas at the temperatures and pressures probed by the observations. These are H₂O, TiO, Na, and K over the wavelength range of the WFC3 spectrum (Burrows & Sharp 1999; Fortney et al. 2008). The remaining gas is assumed to be a solar composition mixture of H₂ and He. We discuss results for including additional molecular species in Section 5.1.2.

In addition to the molecular opacities, we included opacity from clouds and haze. Previous analyses of WASP-12b's transmission spectrum have suggested that hazes are present in the atmosphere, but the haze composition is poorly constrained (Sing et al. 2013; Stevenson et al. 2014c). We therefore used a flexible model that includes a power law haze and an opaque gray cloud deck. We modeled clouds as a gray opacity source that masks transmission through the atmosphere below a fixed pressure level P_c. The haze opacity was parameterized by σ(λ) = σ₀ (λ/λ₀)^γ, where the scattering amplitude σ₀ and slope γ are free parameters (as in Lecavelier Des Etangs et al. 2008). The scattering is presumed to happen throughout the entire atmosphere as if it were a well mixed gas. The scattering amplitude is a scaling to the H₂ Rayleigh scattering cross section times the solar H₂ abundance at 0.43 μm. More sophisticated models for clouds and haze would require additional free parameters that are not justified by the precision of our data. Our simple parameterization is sufficient to capture degeneracies between clouds/haze and the water abundance, which is the primary goal of this investigation.

We also retrieved a "scale height temperature" parameter T_s. We assumed the atmosphere is isothermal, with temperature equal to T_s at all pressures. We tested fitting for a more complex temperature–pressure profile but found it did not affect our results because the transmission spectrum is only weakly sensitive to the atmosphere's thermal structure (see Barstow et al. 2013, 2014).

The final free parameter in the retrieval was a scale factor for the planet's 10 bar radius. We assumed a baseline planet radius of 1.79 R_J and a stellar radius of 1.57 R_⊙ from Hebb et al. (2009). The planet radius was scaled by a factor R_scale to account for uncertainty in the pressure level in the atmosphere at a given radius. To first order, the effect of this scaling is to shift the model transit depths by a constant factor. A second order effect is that scaling the radius changes the amplitude of spectral features (see Equation (1) in Lecavelier Des Etangs et al. 2008).

In sum, the retrieval had 8 free parameters: the abundances of H₂O, Na+K (fixed at the solar abundance ratio), and TiO, as well as clouds and haze, the scale height temperature, and the planet radius scale factor. We refer to this 8-parameter fit as the FULL model. The best fits for the FULL model and nested models within it are shown in Figure 10. The best-fit models are those that produced the lowest χ² values in the MCMC. We list the retrieved water abundances, temperatures, and best fit χ² values for these models in Table 4. The distributions of retrieved parameters from the FULL model are shown in Figure 11. We note that all of the nested models are nearly indistinguishable near the water absorption feature at 1.4 μm. The best constrained molecular abundance is that of water: we retrieved a water volume mixing ratio (VMR) of 1.6 × 10⁻⁴ − 2.0 × 10⁻² at 1σ, which is consistent with expectations for a solar composition gas. The other molecular abundances are not as well constrained. We retrieved a 3σ upper limit on the VMR of TiO equal to 2 × 10⁻⁴, and the abundance of Na+K is unbounded. The cloud and haze properties are also poorly constrained. We discuss the implications of these measurements for the atmosphere composition in detail in Section 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.1.1. Detection Significances

5.1.2. Retrievals with Additional Molecules/Data

In addition to retrievals with the FREE parameterization (described in Section 5.1), we also developed a reparameterized model to retrieve atmospheric properties that are consistent with chemical equilibrium. Rather than varying the mixing ratios of individual gases, we made C/O and metallicity free parameters and computed thermochemical equilibrium abundances along the temperature/pressure (T/P) profile for major species. We retrieved the T/P profile rather than a scale height temperature in order to explore any additional degeneracies due to the profile shape. The T/P profile is modeled with a 5-parameter analytic function (Line et al. 2013b) which produces physically realistic thermal profiles consistent with radiative equilibrium (e.g., Guillot 2010; Heng et al. 2012; Robinson & Catling 2012; Parmentier & Guillot 2014). We calculated chemical profiles for H₂, He, H₂O, CH₄, CO, CO₂, NH₃, H₂S, PH₃, C₂H₂, HCN, Na, K, FeH, TiO, VO, and N₂ to cover the full range of gases that could contribute significant opacity in the near-IR. The equilibrium abundances are computed using the NASA Chemical Equilibrium with Applications code (McBride & Gordon 1996; Line et al. 2011; Moses et al. 2011). The chemical profiles were fed into the radiative transfer model to calculate model transmission spectra for comparison with our data. The model also included opacity from clouds and haze (using the same formalism described in Section 5.1), and the reference radius R_p, for a total of 11 free parameters. While we reduced the number of free absorber parameters to just C/O and metallicity, the overall number of free parameters increased due to the more flexible T/P profile.

We ran retrievals using three different priors for C/O. The fiducial case had an uninformative prior constraint on C/O. We also ran retrievals that constrained the atmospheric composition to be either oxygen-rich or carbon-rich. For the O-rich scenario, the prior probability was set to zero for C/O values greater than unity. Correspondingly, the C-rich scenario had zero prior probability for C/O < 1. We show the best fit models for all three cases in Figure 10. Figure 12 shows marginalized distributions of C/O and metallicity for each scenario, as well as constraints on the temperature–pressure profiles. Table 5 gives the χ² values of the best fits in each scenario. It also lists the 68% credible intervals for C/O, metallicity, and the temperature and water abundance at 1 mbar pressure.

Zoom In Zoom Out Reset image size

The fiducial and oxygen-rich scenarios give nearly identical constraints on C/O and the atmospheric metallicity. Even though the fiducial model has an uninformative prior on C/O, 100% of the retrieved C/O values are less than unity for this case (suggesting the result is data-driven rather than prior-driven). Both scenarios yield constraints on the molecular abundances and thermal profile that agree well with results from the FULL model. The median retrieved water abundances are 6.3 × 10⁻⁴ and 1.0 × 10⁻³ for the fiducial and O-rich scenarios, respectively. The median 1 mbar temperatures are 1410 and 1450 K. We retrieved metallicities in the range 0.3–30× solar (at 1σ). The retrieved cloud and haze properties for these scenarios are consistent with the FULL model at 1σ.

By contrast, the carbon-rich scenario produced significantly different atmospheric properties. To reproduce the water absorption feature in the spectrum, the temperature was driven to much lower values (the median is 320 K at 1 mbar). These low temperatures allow for higher water abundances (see Figure 13), though they are unlikely for the terminator region given how highly irradiated the planet is. The median water abundance, 1.4 × 10⁻⁴, is comparable to that for the O-rich and fiducial scenarios. On the other hand, the median methane abundance increases to 2.4 × 10⁻³, versus 1.3 × 10⁻¹¹ for the fiducial case. These differences have two main effects on the model transmission spectra. One is that absorption features have smaller amplitude because lower temperatures decrease the atmospheric scale height. The second effect is that methane absorption is present in the spectrum. It is especially noticeable in the window between water features at >1.6 μm wavelengths (see the best-fit C-rich model in Figure 10). The consequence of these changes is that C-rich models do not fit the measured spectrum as well as O-rich models. We discuss the strength of the evidence for one scenario over the other in quantitative detail in Section 6.

Zoom In Zoom Out Reset image size

As a test, we also considered a scenario with clouds and haze removed. We computed the Bayes factor for this nested model and found a detection significance of 3.7σ for the cloud and haze parameters. The detection significance is higher than for the FREE model parameterization because the assumption of chemical equilibrium breaks the degeneracy between the cloud/haze and the other molecular abundances (Na+K, TiO); i.e., the other molecular abundances can not increase arbitrarily to make up for the absence of clouds/haze.

We also compared the results from CHIMERA with output from an independent retrievel code to test the robustness of our measurements against a different modeling approach (optimal estimation versus MCMC). We fit the high-resolution WFC3 spectrum with the NEMESIS code (Irwin et al. 2008). The model included H₂O as the only molecular opacity source. We modeled the atmosphere's thermal structure as a scalar multiple of the dayside temperature–pressure profile from Stevenson et al. (2014b). NEMESIS does not currently incorporate cloud-top pressure as a free parameter, so we fixed the altitude of an opaque gray cloud deck over a grid of pressures ranging from 10 to 10⁻³ mbar and retrieved the atmospheric properties for each case (for further discussion of retrieving cloud properties with NEMESIS, see Barstow et al. 2013). The best fit model had a cloud deck at 1 mbar and a retrieved water abundance of 4.0 × 10⁻⁴–1.7 × 10⁻³ at 1σ. Models with higher altitude cloud decks (0.1–0.001 mbar) achieved nearly as good a fit, and had 1σ upper bounds on the water abundance of 3 × 10⁻². These results are in excellent agreement with the results from CHIMERA.

Broadly speaking, a carbon-rich atmosphere is expected to have lower water abundance than an oxygen-rich atmosphere, because most of the oxygen atoms are bound in CO in chemical equilibrium at high temperatures (Madhusudhan et al. 2011b). Our water abundance measurement is a qualitatively better match to predictions for an oxygen-rich composition. In Figure 13, we show the retrieved H₂O abundance for the FULL retrieval and nested models in comparison with thermochemical equilibrium predictions for oxygen-rich (C/O = 0.55 = solar) and carbon-rich (C/O = 1) atmospheres. The predicted abundance depends on pressure and temperature, so we plot the span of predictions over pressures from 0.1 to 10 mbar as a function of temperature.

To quantitatively compare our measurement to the models, we must marginalize over the uncertainty in the pressure and temperature probed by the observations. To do this, we compute an equivalent pressure level p_eq corresponding to an optical depth τ = 0.56 at 1.4 μm at each step in the MCMC chain. This quantity is representative of the typical pressure level at which photons are absorbed by water molecules (Lecavelier Des Etangs et al. 2008). Figure 14 shows the distribution of equivalent pressures from the MCMC chain. Note that p_eq is not a retrieved quantity, but rather derived from the opacities at each MCMC step. We show the distribution of equivalent pressures obtained from this method in Figure 14. The pressures have a 1σ range of 0.05–5 mbar, with a peak at 0.5 mbar.

Zoom In Zoom Out Reset image size

For a specified atmospheric composition, we can then calculate the predicted equilibrium H₂O abundance for a temperature and pressure equal to T_s and p_eq. This procedure yields a "calculated" water abundance to compare with the "observed" water abundance at each step in the MCMC chain. We show the distribution of observed minus calculated ("O − C") values for the oxygen-rich and carbon-rich models in Figure 13 (right panel).

To test how well each model agrees with the retrieved H₂O abundance, we assess whether the O–C distribution is consistent with zero. For the oxygen-rich model, zero is contained in the 1 σ credible interval centered on the median, indicating that this model is a good match. By contrast, the O–C distribution for the carbon-rich model is in tension with zero at approximately 2σ (the 95% credible interval is 0.1–6.4). The median O–C value is 3.7, implying that the typical retrieved water abundances are nearly four orders of magnitude larger than predicted for a carbon-rich composition.

We emphasize that these results are strongly dependent on the temperature and pressure probed by the observations. Assuming a different local thermal profile can result in order-of-magnitude differences in the predicted abundances (cf. Figure 2, Madhusudhan 2012). It is therefore essential to account for the temperature and pressure (and their uncertainties) when estimating C/O.

We obtain more stringent constraints on the atmospheric C/O from the C-C retrieval than from the FREE parameterization. By every metric we enumerate below, C/O values greater than one are ruled out at high confidence. First, the entire posterior probability distribution for C/O from the fiducial model is less than one. We calculate that there are roughly 3000 independent samples in the MCMC chain for that model, which is a sufficient number to rule out C/O > 1 at greater than 3σ confidence. Second, the χ² values for the best fit models provide additional evidence in favor of the O-rich scenario. The best fit O-rich model has χ² = 2.75 versus χ² = 29.2 for the best fit C-rich case (for 2 degrees of freedom). This difference in χ² implies the O-rich model is 10⁶ times more likely than the best fit C-rich model, assuming Gaussian statistics. Third, the Bayes factor (the ratio of integrated posterior probabilities) for the C-rich scenario to the O-rich scenario is 14, which constitutes strong evidence in favor of the O-rich model (Jeffreys 1998). In addition, the fact that the temperature range drops unphysically low—below the condensation temperature of water—is further evidence that the carbon-rich model is not appropriate for these data.

Taken together, these results rule out a carbon-rich composition for the atmosphere at high confidence. This is a more definitive constraint than what we obtained from studying the water abundance alone because the model is sensitive to both the presence of water and the absence of absorption features from other molecules. For example, the fact that no methane features are detected in the spectrum strengthens the case for an O-rich composition beyond the constraints from the presence of water alone. The C-C model also assumes more prior knowledge by imposing chemical equilibrium in the calculation of the model spectra.