HST hot-Jupiter transmission spectral survey: detection of potassium in WASP-31b along with a cloud deck and Rayleigh scattering Purchased

Abstract

1 INTRODUCTION

Hot Jupiters are gas giant planets found in tight orbits around their parent stars, with a large number now known to transit their host star and are amenable to atmospheric investigations. This wide class of exoplanets span a very large range in equilibrium temperatures, with most planets ranging from 1000 to 3000 K depending on the albedo and efficiency of heat transfer from the day to night side. At these temperatures, theoretical models of hot Jupiters have predominately predicted atmospheres which are dominated by water opacity in the near-IR, with the optical containing alkali features and scattering by hydrogen molecules (e.g. Seager & Sasselov 2000; Sudarsky, Burrows & Hubeny 2003; Burrows et al. 2010; Fortney et al. 2010).

For highly irradiated hot-Jupiter atmospheres, the expectation of alkali metals leads to very low predicted albedos. This expectation has generally been backed up by observations, both from the detection of Na (Charbonneau et al. 2002; Redfield et al. 2008; Sing et al. 2008, 2012; Jensen et al. 2011; Huitson et al. 2012; Zhou & Bayliss 2012; Nikolov et al. 2014b) and K (Sing et al. 2011b) atoms as well as the generally low albedos found in many hot Jupiters such as HD 209458b (Rowe et al. 2008).

The atmospheric temperatures of the hot Jupiters are close to the condensation temperature of several abundant components, including silicates and iron, making clouds and hazes a natural outcome of chemistry in much the same way as abundant H₂O, CO, and CH₄ molecules in gaseous form. The possible presence of such condensation clouds was considered early on (Seager & Sasselov 2000; Hubbard et al. 2001). According to models, condensates would weaken spectral features, or mask some of them, depending on the height of the cloud deck (Marley et al. 1999; Sudarsky et al. 2003; Fortney 2005). There is a growing body of evidence for substantial clouds and hazes in hot Jupiters, including the transmission spectra of HD 189733b WASP-12b, HAT-P-32b, and WASP-6b (Sing et al. 2011a, 2013; Gibson et al. 2013; Jordan et al. 2013; Pont et al. 2013; Nikolov et al. 2014a).

Here, we present results for WASP-31b from a large Hubble Space Telescope (HST) programme (HST GO-12473; P.I. Sing), consisting of an eight planet optical atmospheric survey of transiting hot Jupiters. The overall programme goals are to detect atmospheric features across a wide range of hot-Jupiter atmospheres enabling comparative exoplanetology, detect stratosphere causing agents like TiO (Hubeny, Burrows & Sudarsky 2003; Fortney et al. 2008), and detect alkali atoms as well as hazes and clouds. In the results reported so far for this survey, Huitson et al. (2013) detected H₂O in WASP-19b and found TiO was unlikely. Water was also found in HAT-P-1b (Wakeford et al. 2013) with the broad-band spectrum also showing Na i and a significant optical absorber but no K i absorption (Nikolov et al. 2014b). Finally, TiO was ruled out in the very hot planet WASP-12b and the atmosphere was found to be dominated by aerosols, with corundum a possible candidate that is particularly relevant given its high condensation temperature applicable to this very hot Jupiter (Sing et al. 2013).

WASP-31b was discovered by Anderson et al. (2011) as part of the Wide Angle Search for Planets survey. With a low mass (0.46 M_Jup) and large inflated radius (1.55 R_Jup) the planet is of particular interest to atmospheric studies as it has a very low density (ρ_pl = 0.129 ρ_Jup) and correspondingly low surface gravity of 4.56 m s⁻¹ making the planet very favourable to transmission spectroscopy. The exoplanet orbits at 0.047 au around a F6V 6300 K low metallicity ([Fe/H] = −0.19) low activity 3–5 Gyr star, and has a nominal (zero albedo, full recirculation) equilibrium temperature of 1570 K (Anderson et al. 2011). In this paper, we present new HST transit observations with the Space Telescope Imaging Spectrograph (STIS) and Wide Field Camera 3 (WFC3) in spacial scanning mode, combining them with Spitzer photometry to construct a high signal-to-noise (S/N) near-UV to infrared transmission spectrum of WASP-31b, capable of detecting and scrutinizing atmospheric constituents. We describe our observations in Section 2, present the analysis of the transit light curves in Section 3, discuss the results in Section 4 and conclude in Section 5.

2 OBSERVATIONS

2.1 HST STIS spectroscopy

The overall observational strategy for WASP-31b is nearly identical for each of the eight targets in the large HST programme, including the WASP-19b, HAT-P-1b, and WASP-12b which have been published in Huitson et al. (2013), Wakeford et al. (2013), Nikolov et al. (2014b) and Sing et al. (2013).

We observed two transits of WASP-31b with the HST STIS G430L grating during 2012 June 13 and 26, as well as one transit with the STIS G750L during 2012 July 10. The G430L and G750L data sets each contain 43 spectra, which span five spacecraft orbits. The G430L grating covers the wavelength range from 2900 to 5700 Å, with a resolution R of λ/Δλ = 530–1040 (∼2 pixels; 5.5 Å). The G750L grating covers the wavelength range from 5240 to 10 270 Å, with a R = 530–1040 (∼2 pixels; 9.8 Å). Both the STIS data were taken with a wide 52 arcsec × 2 arcsec slit to minimize slit light losses. The visits of HST were scheduled such that the third and fourth spacecraft orbits contain the transit, providing good coverage between second and third contact, as well as an out-of-transit baseline time series before and after the transit. Exposure times of 279 s were used in conjunction with a 128-pixel wide sub-array, which reduces the readout time between exposures to 21 s, providing a 93 per cent overall duty cycle.

The data set was pipeline-reduced with the latest version of calstis, and cleaned for cosmic ray detections with a customized procedure based on median-combining difference images following Nikolov et al. (2014b). The G750L data set was defringed using contemporaneous fringe flats. The mid-time of each exposure was converted into BJD_TBD for use in the transit light curves (Eastman, Siverd & Gaudi 2010).

The spectral aperture extraction was done with iraf using a 13-pixel-wide aperture with no background subtraction, which minimizes the out-of-transit standard deviation of the white-light curves. The extracted spectra were then Doppler-corrected to a common rest frame through cross-correlation, which helped remove sub-pixel wavelength shifts in the dispersion direction. The STIS spectra were then used to create both a white-light photometric time series (see Fig. 1), and custom wavelength bands covering the spectra, integrating the appropriate wavelength flux from each exposure for different bandpasses.

Figure 1.

HST STIS normalized white-light curves for the three WASP-31b transits taken with the STIS G430L, STIS G750L, and WFC3. Top panels: raw light curves normalized to the mean out-of-transit raw flux. Middle panels: detrended light curves along with the best-fitting limb-darkened transit model superimposed with continuous lines. Lower panels: observed minus modelled light-curve residuals, compared to a null (dashed lines) and 1σ level (dotted lines).

Open in new tabDownload slide

Like the other targets in our HST programme, the STIS light curves of WASP-31 show instrument-related systematic effects primarily due to the well-known thermal breathing of HST, which warms and cools the telescope during the 96 min day/night orbital cycle, causing the focus to vary. The first exposure of each spacecraft orbit has also consistently been found to show significantly lower fluxes than the remaining exposures. Similar to our previous STIS analysis (Huitson et al. 2013; Sing et al. 2013; Nikolov et al. 2014b), we intended to discard the entire first orbit to minimize the effect of the breathing trend, but found for two of the three HST visits, a few of the exposures taken towards the end of the first orbit could be used in the analysis and fit with the same systematics model as used for the subsequent orbits. In addition, we set the exposure time of the first image of each spacecraft orbit to be 1 s in duration, such that it could be discarded without significant loss in observing time. We find that the second exposures taken for all five spacecraft orbits (each 279 s in duration) do not in general show the first-exposure systematic trend, with an exception of the third orbit of the G750L light curve. Upon analysing several targets from the large HST programme, the 1-s strategy to avoid the faint exposure systematic appears to have mixed results. For both WASP-19b, WASP-12b, and here for WASP-31b the second exposures of each spacecraft orbit do not generally show the first-exposure trend as intended, though the exposures of HAT-P-1b were still affected. Though the 1-s strategy may not always be effective, especially for brighter targets, it still appears to be a recommended practice for STIS programmes requiring high photometric accuracy and long exposures as it can often save useful observing time as was usually the case here.

2.2 HST WFC3 spectroscopy

In addition to the STIS data, observations of WASP-31b were also conducted in the infrared with WFC3 on the HST. Observations began on 2012 May 13 at 12:53 using the IR G141 grism in forward spatial scan mode over five HST orbits. The spacial scanning is done by slewing the telescope in the cross-dispersion direction during integration in a similar manner for each exposure, which increases the duty cycle and greatly increases the number of counts obtained per exposure. We used a 256 × 256 pixel subarray with the SPARS25 NSAMP = 8 readout sequence, which gives exposure times of 134.35 s. Our observations were conducted with a spatial scan rate of ∼0.15 pixels s⁻¹, where 1 pixel = 0.13 arcsec, with the final spatially scanned spectrum spanning ∼20 pixels. The G141 grism images contain only the first-order spectrum, which spans 143 pixels with a dispersion of 4.63 nm pixel⁻¹ in the direction of the x-axis. The resultant spectra have average count rates of about 3.8 × 10⁴ electrons pixel⁻¹.

We used the ‘ima’ outputs from the CALWFC3 pipeline. For each exposure, calwf3 conducts the following processes: reference pixel subtraction, zero-read and dark current subtraction, and a non-linearity correction; the resultant images are in units of electrons per second. For the spectral extraction, we trimmed a 152 × 80 box around each spectral image and the spectra were extracted using custom idl routines, similar to IRAF's APALL procedure, with an aperture of ±12 pixels determined by minimizing the standard deviation of the fitted white-light curve. The aperture is traced around a computed centring profile, which was found to be consistent in the y-axis with an error of 0.06 pixels. Background subtraction was applied using a clean region of the untrimmed image, with background accounting for 0.4 per cent of the total flux measured in each wavelength column. Subsequent data analysis was conducted with the first orbit removed (16 exposures), as it suffers from significant thermal breathing systematic effects similar to the STIS data.

For wavelength calibration, direct images were taken in the F139M narrow-band filter at the beginning of the observations. We assumed that all pixels in the same column have the same effective wavelength, as the spatial scan varied in the x-axis direction by less than one pixel, resulting in a spectral range from 1.1 to 1.7 μm. This wavelength range is later restricted for the spectroscopic light curve fits as the strongly sloped edges of the grism response result in much lower S/N light curves.

2.3 Spitzer IRAC photometry

We analyse two transit observations obtained using the Infrared Array Camera (IRAC) instrument (Programme 90092 with P.I. Desert) on the Spitzer space telescope in the 3.6 μm and 4.5 μm channels in subarray mode (32 × 32 pixel, or 39 centred on the planets host). The 3.6 μm observation was performed on ut 2013 March 9 (between 06:59 and 11:37) and the 4.5 observation was performed on ut 2013 March 19 (between 12:19 and 16:58), each transit containing 8320 images (see Fig. 2).

Figure 2.

Spitzer photometry of WASP-31b for the (left) 3.6 μm and (right) 4.5 μm channels. Plotted is the raw flux along with the best fitting transit model including systematic effects (top), the flux corrected for detector trends and binned by 8 min (middle), and residuals (bottom).

Open in new tabDownload slide

Photometry was extracted from the basic calibrated data files. These images are stacked in FITS data cubes, containing 64 exposures taken in a sequence which were produced by the IRAC pipeline (version S19.1.0) after dark subtraction, flat-fielding, linearization, and flux calibration. Both observations have effective integration times of 1.92 s per image.

We converted the images to photon counts (i.e. electrons) by multiplying the images by the gain and individual exposure times (FITS HEADER key words SAMPTIME and GAIN) and dividing by the flux conversion factor (FLUXCONV). Timing of each image was computed using the UTC-based Barycentric Julian Date (BJD_UTC) from the FITS header keyword BMJD_OBS, transforming these time stamps into BJD based on the BJD terrestrial time standard (Eastman et al. 2010). Both channels show a strong ramp feature at the beginning of the data, and we elected to trim the first 22 min of data to allow the detector to stabilize.

We performed outlier filtering for hot (energetic) or cold (low-count values) pixels in the data by examining the time series of each pixel. This task was performed in two passes, first flagging all pixels with intensity 8σ or more away from the median value computed for each frame from the 10 surrounding images. The values of these flagged pixels were replaced with the local median value. In the second pass, we flagged and replaced outliers above the 4σ level, following the same procedure. The total fraction of corrected pixels was 0.26 per cent for the 3.6 μm and 0.07 per cent for the 4.5 μm.

We estimated and subtracted the background flux from each image of the time series. To do this, we performed an iterative 3σ outlier clipping for each image to remove the pixels with values associated with the stellar point spread function (PSF), background stars or hot pixels, created a histogram from the remaining pixels, and fitted a Gaussian to determine the sky background.

We measured the position of the star on the detector in each image incorporating the flux-weighted centroiding method¹ using the background subtracted pixels from each image for a circular region with radius 3 pixels centred on the approximate position of the star. While we could perform PSF centroiding using alternative methods such as fitting a two-dimensional Gaussian function to the stellar image, previous experiences with warm Spitzer photometry showed that the flux-weighted centroiding method is either equivalent or superior (Beerer et al. 2011; Knutson et al. 2012; Lewis et al. 2013). The variation of the x and y positions of the stellar PSF on the detector were 0.17 and 0.18 pixels for the 3.6 μm channel and 0.76 and 0.24 pixels for the 4.5 μm channel.

We extract photometric measurements from our data following two methods. First, aperture photometry was performed with idl routine aper using circular apertures ranging in radius from 1.5 to 3.5 pixels in increments of 0.1. The best result was selected by measuring the flux scatter of the out-of-transit portion of the light curves for both channels after filtering the data for 5σ outliers with a width of 20 data points.

2.4 Photometric monitoring for stellar activity

We obtained nightly photometry of WASP-31 to monitor and characterize the stellar activity over the past three observing seasons with the Tennessee State University Celestron 14-inch (C14) Automated Imaging Telescope (AIT) located at Fairborn Observatory in Arizona (see e.g. Henry 1999; Eaton, Henry & Fekel 2003). The AIT uses an SBIG STL-1001E CCD camera and exposes through a Cousins R filter. Each nightly observation consisted of 4–10 consecutive exposures on WASP-31 and several comparison stars in the same field of view. The individual nightly frames were co-added and reduced to differential magnitudes in the sense WASP-31 minus the mean brightness of six constant comparison stars. Each nightly observation has been corrected for bias, flat-fielding, pier-side offset, and for differential atmospheric extinction.

A total of 227 successful nightly observations (excluding a few isolated transit observations) were collected between 2011 October 31 and 2014 January 28. The individual observations are plotted in the top panel of Fig. 3 and summarized in Table 1. The standard deviations of single observations, with respect to their corresponding seasonal means, range between 0.0031 and 0.0038 mag for the three observing seasons (Table 1; column 4). This is near the limit of our measurement precision with the C14, as determined from the constant comparison stars in the field. The seasonal means given in column 5 agree to within a standard deviation of only 0.001 mag. Since we do not standardize the ensemble differential magnitudes, this is probably consistent with the absence of year-to-year variability.

Figure 3.

Top: the composite nightly Cousins R-band photometric data set of WASP-31 from the 2011–2012, 2012–2013, and 2013–2014 observing seasons, acquired with the C14 AIT at Fairborn Observatory. The data within each observing season scatter from night to night with a standard deviation between 0.003 and 0.004 mag. This is the approximate limit of our measurement precision. The scatter of the yearly means is only 0.001 mag. Middle: the frequency spectrum of the composite data set, normalized to the first seasonal mean, suggests low-amplitude variability at about 8.6 d. Bottom: the normalized data set phased to the 8.6 d period, which we interpret as the stellar rotation period. A least-squares sine fit to the 8.6-d rotation period gives a peak-to-peak amplitude of 0.003 mag.

Open in new tabDownload slide

We performed a periodogram analysis of the composite data set in the top panel of Fig. 3 after first normalizing the data by translating the second and third observing seasons to have the same mean as the first season. To first order, this removes low-amplitude, long-term variability in WASP-31 or the comparison stars as well as any uncorrected year-to-year systematic errors. A frequency spectrum of the normalized data set, computed as the reduction in variance of the data versus trial frequency, is plotted in the middle panel of Fig. 3. Trial frequencies ranged between 0.005 and 0.95 c/d, corresponding to a period range between 1 and 200 d. The result suggests low-amplitude brightness variability with a period of 8.6 ± 0.1 d. This is consistent with the stellar rotation period of 7.9 ± 0.7 d predicted by Anderson et al. (2011) from their observed vsin i and computed stellar radius R_*. Therefore, we take 8.6 d as the stellar rotation period made visible by rotational modulation in the visibility of weak surface features.

The three-year normalized data set is plotted in the bottom panel of Fig. 3, phased to the rotation period of 8.6 d. A least-squares sine fit gives a peak-to-peak amplitude of 0.0031 ± 0.0006 mag. Thus, while we find probable low-amplitude rotational modulation in stellar brightness, we find no evidence for significant variability above the noise at or near the planetary orbital period (corresponding to a frequency of 0.2936 c/d in Fig. 1). A least-squares sine fit of the normalized data to the orbital period gives a peak-to-peak amplitude of only 0.0014 ± 0.0006 mag. Since the 8.6-d stellar rotation period and the 3.4-d planetary orbital period are well separated, the radial velocity variations are clearly due to planetary reflex motion with very little jitter caused by surface activity (e.g. Queloz et al. 2001; Paulson et al. 2004).

In addition, we also measure a low stellar activity for WASP-31 as determined by the Ca ii H & K lines, where we find an activity index of |${\rm log}(R^{\prime }_{{\rm HK}})=-5.4$| using archival High Accuracy Radial Velocity Planet Searcher (HARPS) spectra.

From equation 5 of Sing et al. (2011a), activity levels of 0.0031 mag limits the effects of the photometric variability on the measured planet-to-star radius contrast to 0.000 18 R_p/R_*. With low levels of stellar activity, we conclude their effects on measuring WASP-31b's transmission spectrum are minimal. A low-activity star in this case is particularly advantageous for transmission spectroscopy, as it minimizes or eliminates any potential uncertainty in the interpretation of the transmission spectrum across the spectral regime and from different data sets and epochs.

3 ANALYSIS

3.0.1 STIS white-light curve fits

Our overall analysis methods are similar to the HST analysis of Sing et al. (2011a, 2013), Huitson et al. (2013) and Nikolov et al. (2014b), which we also describe briefly here. The light curves were modelled with the analytical transit models of Mandel & Agol (2002). For the white-light curves, the central transit time, inclination, stellar density, planet-to-star radius contrast, stellar baseline flux, and instrument systematic trends were fit simultaneously, with flat priors assumed. The period was initially fixed to a literature value, before being updated, with our final fits adopting the value from Section 3.1. Both G430L transits were fit simultaneously with a common inclination, stellar density, and planet-to-star radius contrast. The results from the G430L and G750L white-light curve fits were then used in conjunction with literature results to refine the orbital ephemeris and overall planetary system properties (see Section 3.1).

To account for the effects of limb-darkening on the transit light curve, we adopted the four parameter non-linear limb-darkening law, calculating the coefficients with ATLAS stellar models (T_eff = 6250, log g = 4.5, [Fe/H] = −0.2) following Sing (2010). Varying the stellar model within the known parameter range has a negligible effect on the output limb-darkening coefficients and fit transit parameters.

As in our past STIS studies, we applied orbit-to-orbit flux corrections by fitting for a fourth-order polynomial to the photometric time series, phased on the HST orbital period. The systematic trends were fit simultaneously with the transit parameters in the fit. Higher order polynomial fits were not statistically justified, based upon the Bayesian information criteria (BIC; Schwarz 1978). The baseline flux level of each visit was let free to vary in time linearly, described by two fit parameters. In addition, for the G750L we found it justified by the BIC to also linearly fit for two further systematic trends which correlated with the x and y detector positions of the spectra, as determined from a linear spectral trace in iraf.

The errors on each data point were initially set to the pipeline values, which is dominated by photon noise but also includes readout noise. The best-fitting parameters were determined simultaneously with a Levenberg–Marquardt least-squares algorithm (Markwardt 2009) using the unbinned data. After the initial fits, the uncertainties for each data point were rescaled based on the standard deviation of the residuals and any measured systematic errors correlated in time (‘red noise’), thus taking into account any underestimated errors in the data points.

The red noise was measured by checking whether the binned residuals followed a N^−1/2 relation, when binning in time by N points. In the presence of red noise, the variance can be modelled to follow a |$\sigma ^2 = \sigma _{\rm w}^2/N+\sigma _{\rm r}^2$| relation, where σ_w is the uncorrelated white noise component, and σ_r characterizes the red noise (Pont, Zucker & Queloz 2006). We find that the pipeline per-exposure errors are accurate at small wavelength bin sizes, which are dominated by photon noise, but are in general somewhat underestimated at larger bin sizes. For the G750L white-light curve we find our light curve precision achieves 1.3 × the photon noise value with the scatter of the residuals measured at 201 ppm, while the G430L achieves 1.6 times the photon noise value with the residual scatter 276 ppm. We do not find any evidence for remaining systematic errors, as both G430L STIS white-light curves as well as the G750L curve give σ_r = 0. A few deviant points were cut at the 3σ level in the residuals.

The uncertainties on the fitted parameters were calculated using the covariance matrix from the Levenberg–Marquardt algorithm, which assumes that the probability space around the best-fitting solution is well described by a multivariate Gaussian distribution. Previous analyses of other HST transit observations (Berta et al. 2012; Line et al. 2013; Sing et al. 2013; Knutson et al. 2014b; Nikolov et al. 2014b) have found this to be a reasonable approximation. We then compared these uncertainties to the results of an MCMC analysis (Eastman, Gaudi & Agol 2013), which does not assume any functional form for this probability distribution. In each case, we found equivalent results between the MCMC and the Levenberg–Marquardt algorithm for both the fitted parameters and their uncertainties. Inspection of the 2D probability distributions from both methods indicate that there are no significant correlations between the systematic trend parameters and the planet-to-star radius contrast.

3.0.2 WFC3 white-light curve fits

WFC3 is now being used extensively for transiting exoplanet work (Berta et al. 2012; Gibson et al. 2012; Deming et al. 2013; Huitson et al. 2013; Mandell et al. 2013; Swain et al. 2013; Wakeford et al. 2013; Knutson et al. 2014a; Kreidberg et al. 2014). The data quality appears mainly affected by detector persistence, which causes a ramp-effect at high count levels and HST thermal-breathing variations. Thus far, parametrized methods and various common-mode techniques have been employed to handle the systematics, which largely appear repeatable orbit-to-orbit.

Our WFC3 white-light curve transit shows small intra-orbit and orbit-to-orbit trends (see Fig. 1). We fit for the systematics errors similarly as the STIS data, adopting a fourth-order polynomial of the HST orbital phase, a baseline flux allowed to vary linearly in time, and a first-order shift along the y-axis detector position as measured through cross-correlation. We found similar transit depths with the divide-OOT method (Berta et al. 2012) but found larger overall noise properties, with the BIC favouring the parametrized method (ΔBIC = 27, see Fig. 4).

Figure 4.

HST WFC3 detrended white-light curves along with the best-fitting limb-darkened transit models for the cases of detrending with a polynomial of the HST orbital phase (black circles) and the divide-OOT method (red squares). The vertical extent of the plotted symbols reflects the flux uncertainty.

Open in new tabDownload slide

We find our white-light curve precision achieves about 1.47 times the photon noise value with the scatter of the residuals measured at 157 ppm. There is no evidence for remaining systematic errors, as we estimate a value of σ_r = 0 from the binning technique. With these precision levels, the near-IR transit depth is measured to about a 53 ppm accuracy.

3.0.3 Spitzer light curve fits

The fits using both the time-variable and constant aperture methods examined were examined, including the root-mean-square residual after fitting the light curves from each channel as well as the white and red noise components measured using the wavelet method (Carter & Winn 2009). We found the 3.6 μm channel favoured the time-variable aperture while the 4.5 μm channel favoured a constant aperture (see Fig. 2), with both methods giving consistent transit depths for each channel. To compare with the HST transmission spectra, we set the transit ephemeris and system parameters to the values found in Section 3.1 and re-fit the light curves. We find residual red noise in both channels, which is comparable to the binned white noise on the transit-duration time-scale. We estimate values for σ_r of 212 and 346 ppm for channels 1 and 2, respectively, using either the binning or wavelet techniques, which we account for in the final quoted error bars.

3.1 System parameters and transit ephemeris

Figure 5.

O−C diagram for the transit times of WASP-31 using the transits listed in Table 2 and a linear ephemeris. The 1σ error envelope on the ephemeris is plotted as the dotted lines.

Open in new tabDownload slide

We also refined the transit parameters for the inclination, i, and stellar density, ρ_*, using the white-light curve fits of the G750L and G430L, in conjunction with the two Spitzer light curves and the values of Anderson et al. (2011). Given that all our measurements were highly consistent (within the errors) with both each other and with those of Anderson et al. (2011), we computed the weighted mean value of both the inclination, semimajor axis to stellar radius ratio, and the stellar density finding values of i = 84.67 ± 0| $_{.}^{\circ}$|11, a/R_* = 8.19 ± 0.10, and ρ_* = 0.895 ± 0.033 (cgs), respectively.

3.2 Transmission spectrum fits

We extracted various wavelength bins for the STIS G430L, G750L, and WFC3 G141 spectra, to create a broad-band transmission spectrum and search for the expected narrow-band Na and K features (see Figs 6 –9). In these transit fits, we fixed the transit ephemeris, inclination, and stellar density to the results from Section 3.1. The parameter R_P/R_* was allowed to fit freely, as well as the baseline flux and model parameters describing the instrument systematic trends. The limb-darkening parameters were fixed to the model values, with the four non-linear coefficients for each bin individually calculated from the stellar model, taking into account the instrument response (see Table 3). Given our overall lack of ingress or egress data, we were not able to perform a sensitive test of the stellar limb-darkening as was done on WASP-12A (Sing et al. 2013).

Figure 6.

HST STIS G430L spectral bin transit light curves for visit 15. Left: light curves for the HST STIS G430L spectral bins with the common-mode trends removed, overplotted with the best-fitting systematics model. The points have been offset in relative flux, and systematics model plotted with connecting lines for clarity purposes. The light curves are ordered in wavelength, with the shortest wavelength bin shown at the top and longest wavelength bin at the bottom. Middle: light curves fully corrected for systematic errors, with the best-fitting joint transit model over-plotted. Right: residuals plotted with 1σ error bars along with the standard deviation (dotted lines).

Open in new tabDownload slide

Figure 7.

The same as Fig. 6, but for the G430L spectral bins of visit 16.

Open in new tabDownload slide

Figure 8.

The same as Fig. 6, but for the G750L spectral bins.

Open in new tabDownload slide

Figure 9.

The same as Fig. 6, but for the WFC3 G141 spectral bins.

Open in new tabDownload slide

As done for WASP-12b (Sing et al. 2013), we handled the modelling of systematic errors in the spectral bins by two methods. We allowed each light curve to be fit independently with a parametrized model and also separately fit the spectral bins by removing common-mode trends from each spectral bin, before fitting for residual trends with a parametrized model containing fewer free parameters. For both the STIS and WFC3, removing the common trends significantly reduces the amplitude of the observed systematics, as these trends are observed to be similar wavelength-to-wavelength across a detector. The two methods give very similar results, though we choose the common-mode method to quote our subsequent results as it produces slightly higher overall S/N with more dof. The ΔBIC values are on average lower for the common-mode method by 10.5 for the G430L light curves, 13.1 for G750L light curves, and 25.7 for the WFC3 light curves. With the common-mode method, we reach average precisions of 88 per cent the photon noise limit with the G430L, 92 per cent for the G750L, and 72 per cent with the WFC3 with no significant red-noise detected in any of our spectral bins. The transmission spectral results for both HST and Spitzer are given in Table 3 and shown in Fig. 10.

Figure 10.

Combined HST STIS, HST WFC3 transmission spectrum of WASP-31b.

Open in new tabDownload slide

4 DISCUSSION

WASP-31b has a zero-albedo equilibrium temperature of T_eq = 1570 K, and is highly inflated with a radius of R_p = 1.55R_Jup, mass of M_p = 0.48M_Jup, and a surface gravity of g = 4.56 m s⁻² (Anderson et al. 2011). At these temperatures, cloud-free solar-abundance models predict significant optical absorption by Na and K atoms, covering most of the optical spectrum with pressure-broadened line wings, with the near-IR dominated by H₂O (Burrows et al. 2010; Fortney et al. 2010).

4.1 Interpreting the transmission spectrum

The broad-band transmission spectrum of WASP-31b seen in Fig. 10 has three distinct features: a significant blueward slope shortwards of ∼0.6 μm, a largely flat spectrum longwards of about 0.6 μm which extends through the red-optical out to infrared wavelengths, and a strong but narrow K line core. We compare these spectral features to both analytic models and those from the modelling suites of Burrows et al. (2010) and Fortney et al. (2010).

The models from Fortney et al. (2008, 2010) include isothermal models as well as those with a self-consistent treatment of radiative transfer and chemical equilibrium of neutral and ionic species. Chemical mixing ratios and opacities assume solar metallicity and local thermochemical equilibrium accounting for condensation and thermal ionization but no photoionization (Lodders 1999; Lodders & Fegley 2002; Freedman, Marley & Lodders 2008).

The models from Burrows et al. (2010) and Howe & Burrows (2012) use a 1D day-side TP profile with stellar irradiation, in radiative and thermochemical equilibrium. Chemical mixing ratios and opacities assume solar metallicity and local chemical equilibrium accounting for condensation and thermal ionization but no photoionization, using the opacity data base from Sharp & Burrows (2007) and the equilibrium chemical abundances from Burrows & Sharp (1999) and Burrows et al. (2001).

4.1.1 A cloud deck

The large pressure-broadened alkali features predicted by models with a clear atmosphere are not seen in the data (see Fig. 11) and cloud-free solar-abundance models can be ruled out. A 1500 K isothermal model from Fortney et al. (2008), which contains large pressure-broadened alkali line wings, gives a poor fit with a χ² value of 51.0 for 20 dof when fit to the optical HST STIS data between the wavelengths of 0.52 and 0.9 μm (excluding the potassium core). A flat line gives a better fit to the data, at 5σ confidence, with a χ² value of 25.9 for 20 dof

Figure 11.

Plotted is the broad-band transmission spectral data along with atmospheric models including a 1500 K cloud-free solar-abundance model (red) and a Rayleigh scattering and cloud-deck model (purple).

Open in new tabDownload slide

The H₂O bandhead features predicted by cloud-free models in the WFC3 spectra are also not detected, though the data are clearly of sufficient quality to have detected such a feature (see Fig. 11). In particular, a 1500 K isothermal model predicts a broad ΔR_p/R_* = 0.0058 amplitude H₂O feature between 1.25 and 1.425 μm, while the WFC3 spectra does not significantly deviate from a flat spectrum beyond the precision level of the data points which average 0.0011 R_p/R_* per spectral bin. Moreover, the 1500 K isothermal model dominated by molecular H₂O in the near-IR gives χ² value of 112 for 24 d.o.f when fit to the WFC3 data, compared to a χ² value of 34 for 24 dof for a flat line. Thus, a water dominated atmosphere can be ruled out at 8.8σ confidence.

There is not a significant difference in the overall measured radius between the four transits of the WFC3, STIS G750L, and Spitzer data sets. A single flat line fitting these the 47 data points for wavelengths longer than ∼0.52 μm (excluding the Na and K lines) results in a χ² value of 59.8 for 46 dof and R_p/R_* = 0.12519 ± 0.00015.

Approximations to hazy or cloudy atmospheres were generated varying the Fortney et al. and Burrows et al. isothermal models with either an artificially added Rayleigh scattering component, simulating a scattering haze, or a cloud deck simulated by a grey flat line at different altitudes. In these models, the artificial hazes and clouds cover up some but not all of the expected atomic and molecular gaseous species.

While a flat line is a good fit to the flat portions of the transmission spectrum, a low-amplitude water feature does slightly improve the χ² quality of the fit. A cloudy Fortney et al. model gives a χ² value of 31.2 for 23 dof (fitting for the altitude of the cloud deck and overall model altitude) when compared to the WFC3 data, compared to a flat line χ² of 34.2 for 24 dof In addition, when fit to the STIS, WFC3 and Spitzer data for wavelengths >0.52 μm (excluding the Na and K lines) the cloudy Fortney et al. model gives our best fit to those data with a χ² of 56.61 for 44 and a R_P/R_* = 0.12497 ± 0.00014.

4.1.2 Rayleigh scattering and a cloud deck

We elected to fit a two-component model which included both a power-law cross-section at short wavelengths, and a flat spectrum at longer wavelengths, with the transition wavelength, λ_T, between the two freely fit (see Fig. 11). Excluding only the Na and K line cores from the full transmission spectrum, this model gave a good fit to the 59 data points (χ² of 65.4 for 56 dof), with λ_T = 0.510 ± 0.031 μm and αT = −9284 ± 3240 K. The uncertainty in λ_T increases the error in measuring αT by about 40 per cent. The result is compatible with Rayleigh scattering, α = −4, which results in estimated temperatures of 2300 ± 800 K. While the derived temperature may appear too high to be associated with most condensates, we note that given the large error it is still within the 1σ uncertainty of the zero-albedo equilibrium temperature (T_eq = 1570 K) where condensates from Fe and Mg are found. This two-component model, simulating a short-wavelength Rayleigh haze and a long-wavelength cloud deck, results in a much better χ² fit compared to a simple flat line over the same wavelength range (χ² of 87.17 for 57 dof). Thus, the presence of a near-UV Rayleigh slope is significant at about 4σ confidence.

4.1.3 Potassium

In addition to measuring the broad-band transmission spectrum, we also searched the STIS data for narrow-band signatures, notably the line cores of the Na i and K l D lines and the H_α line. No significant absorption was found for either Na i or H_α. We quote a 34 Å wide bin across the Na line in our transmission spectrum for comparison purposes (see Table 3), choosing this width based on the HAT-P-1b Na detection (Nikolov et al. 2014b).

While no significant Na i absorption was detected, there is a strong K line evident in the data (see Fig. 11). Searching both on and around the wavelength range of the K i doublet (with line centres at 7664 and 7699 Å), we measured the transit depth in various bin sizes ranging from one resolution element (∼10 Å) up to several hundred Å, comparing the transit depths to the surrounding continuum regions and that of the white-light curve. We find a radius of R_p/R_* = 0.1334 ± 0.0020 in a 78 Å band (7645 and 7720 Å) which encompasses the K i doublet gives an optimal detection, which is larger than the white-light radius at 4.2σ confidence (ΔR_p/R_* = 0.0086 ± 0.0020). A confirmation of the K i detection is provided by the fact that each line of the potassium doublet is also individually detected, with a band from 7645 to 7683 Å giving radius of R_p/R_* = 0.1342 ± 0.0027 and a band from 7683 to 7720 Å giving a radius of R_p/R_* = 0.1325 ± 0.0032. No significant signals were detected outside the 7645 to 7720 Å region confining the potassium feature to a 78 Å region, as wider bins dilute the measured potassium absorption depth. In addition, no other similar features rising to significant levels above the broad-band continuum are seen at any other wavelengths in the G750L data.

4.2 Limits on the pressure level of the cloud deck

The pressure level of the broad-band spectrum can potentially be constrained by the shape of the potassium line and the presence (or lack) of pressure-broadened wings. The maximum expected pressures probed would correspond to the near-UV slope being due to gaseous H₂ Rayleigh scattering, which following the fits from Section 4.1.2 and Lecavelier des Etangs et al. (2008) would put the cloud deck radius of R_P/R_* = 0.124 97 at pressures of 40 mbar. However, this scenario is unlikely due to the lack of pressure-broadened K wings, and the fact that both Na and H₂O would have to both have sub-solar abundances depleted by about a factor of 100 to hide both features.

In addition to the cloud-free solar-abundance models which can be ruled out at high confidence (see Section 4.1.1), we further test for the presence of collisionally broadened wings by using the Rayleigh and cloud analytic model from Section 4.1.2, fixing both T and λ_T while adding a potassium line and adjusting P_ref via the Rayleigh slope opacity. This also tests the hypothesis that the near-UV slope is due to gaseous H₂ Rayleigh scattering (as opposed to aerosols), but a low-altitude cloud is still present. With a solar abundance of potassium and assuming H₂ Rayleigh scattering, the model gives a fit with a χ² of 71.3 with 60 dof The model reproduces the height of the K core measurement and still contains significant collisional broadening as the presumed low-altitude cloud does not fully cover the broadened alkali wings. Assuming instead the Rayleigh scattering is from an aerosol (residing at higher altitudes and lower pressures) and maintaining the amplitude of the potassium line gives a χ² of 65.2 for 60 dof, with the model lacking significant collisionally broadened wings. In conjunction with the results of Section 4.1.1, pressure-broadened potassium wings can be fully excluded from cloud-free models, and further at about the 98 per cent confidence level for the case of H₂ Rayleigh scattering with a low-altitude cloud. We note that this result can be expected from statistical theory, as the width of the potassium line calculated assuming pressures of zero in equation (8) (i.e. only natural broadening) match very well with the measured width of the potassium line.

We estimate the maximum pressure level of the cloud deck by raising the Rayleigh scattering opacity and K abundance by the minimum amounts needed to maintain the K line amplitude and produce a sufficiently narrow width. A cloud deck at pressures below about 10 mbar reduce the broadening such that the model gives a Δχ² of 1 compared the best-fitting model, which corresponds to a minimum aerosol opacity which is about 4 times that of gaseous H₂.

As the 1.4 μm H₂O feature is strong, the pressure level of the cloud deck estimated by the lack of alkali broadened wings (10 mbar) would still not be sufficient to cover H₂O at the mixing ratios predicted when assuming chemical equilibrium and solar abundances. Based on the Fortney et al. and Burrows et al. models, we estimate that a cloud deck at pressures of about 1 mbar would be needed to cover all but the very peak of the 1.4 μm H₂O feature at mixing ratios assuming solar-abundances.

4.3 The Na/K abundance ratio

To measure the relative abundances, we fit an analytic model which included Rayleigh scattering, a cloud deck as well as Na and K (see Fig. 12). We also included the small H₂O feature from the cloudy Fortney et al. model as found in Section 4.1.1 as it provided the best fits, though very similar results were obtained by removing the opacity of H₂O and fitting the entire infrared wavelength region with a flat line. Fitting for the alkali abundances, T, λ_T and the baseline radius while setting P_ref below 10 mbar gives our best-fitting overall model to all 61 data points, with a χ² of 63.0 for 56 dof We find a Na to K abundance ratio² of ln[Na/K] = −3.3 ± 2.8 which is sub-solar compared to the solar value being ln[Na/K]_⊙ = +2.8 (Asplund et al. 2009). The uncertainty in the abundance ratio is dominated by the uncertainty in the Na line core measurement, which could benefit greatly from higher resolution data clearly separating any Na signatures from the surrounding clouds.

Figure 12.

Plotted is the broad-band transmission spectral data along with atmospheric models. Two solar-composition models containing a scattering haze are shown from the modelling suits of Burrows et al. (green) and Fortney et al. (blue). Our best-fitting model is also plotted (purple) containing a Rayleigh scattering haze, a grey cloud deck at low pressures, non-pressure broadened Na and K features, and an obscured H₂O feature. The band-averaged model points are indicated with open circles.

Open in new tabDownload slide

With the K line included in the fit, the fit temperature increases somewhat to T = 2470 ± 590 K. At the high altitudes probed in the K line core, the presence of a thermosphere is possible and may be responsible for the high measured temperatures, as has been found in HD 189733b and HD 209458b using the Na line (Vidal-Madjar et al. 2011; Huitson et al. 2012). However, the errors on the retrieved temperatures are still high which limits the current interpretation.

A sub-solar value for the Na to K abundance ratio is surprising, and to the best of our knowledge has not yet been seen in any planetary atmosphere. A variety of physical processes could in principle alter the abundances, though a Na/K ratio <1 is in any case still difficult to explain. While the alkali metals are easily ionized, which could dramatically alter their measured neutral abundances, photoionization favours depleting K i over Na i presumably making it a difficult mechanism to explain the observed ratio. Condensation of Na could also potentially play a role in reducing its abundance, and would be consistent with the presence of a cloud deck or higher altitude scatterer. However, the expectations from current chemical models predict Na to condense at temperatures near ∼1000 K (Lodders 1999) which is lower than the expected temperatures of WASP-31b, though a large day/night temperature contrast could be present, making for a Na cold-trap on the night side while K is not likely to condense out at temperatures exceeding ∼800 K. This explanation of condensing Na but not K is very specific: the planet has to have a specific temperature regime so that Na₂S condenses but not KCL, which requires only a ∼200 K lower temperature. Finally, we speculate that the large element-to-element metallicity variations could be primordial in nature, or altered after formation from accretion processes.

While the absolute abundances in WASP-31b are unknown, as P_ref is unknown, limits on the pressure can in principle provide bounds on the K abundance. A cloud deck pressure level of 10 mbar or lower, required to hide the pressure-broadened alkali wings, corresponds to a minimum K abundance of ln[ε_K] = −14.77 ± 0.86, indicating a super-solar K abundance (solar ln[ε_K]_⊙ = −16.05). The K abundance would be a factor of ∼10 higher if the cloud deck is instead at 1 mbar, required to hide a solar-abundance of H₂O.

4.4 Cloud properties and vertical mixing

A broad-band transmission spectrum can give valuable constraints on aerosol particle sizes. Obtaining a flat spectrum out to wavelengths of several microns requires a particle size of at least ∼1 μm; small particles of, e.g. ∼0.1 μm or less would become transparent in the IR and thus would be inefficient at scattering the radiation, as necessary to produce a flat spectrum. The presence of Rayleigh scattering can also give constraints, as scattering occurs when the particle size is less than the wavelength of light probed. In the WASP-31b transmission spectrum, there are two distinct possible scattering slopes in the data, the near-UV to blue slope as discussed in Section 4.1.2, and second possible (though fairly flat) slope in the IR as the 4.5 μm channel exhibits the smallest radii in our spectrum (at ∼2σ significance).

To further estimate the possible particle sizes allowed in the data, we computed a series of transmission spectra assuming an opacity dominated by a single size aerosol calculated with Mie theory, and compared it with the transmission spectrum. As the cloud and haze composition cannot be easily identified, we explored a variety of plausible materials, including Al₂O₃, CaTiO₃, Fe₂O₃, MgSiO₃, and a Titan tholin (Sing et al. 2013). We find a single particle size does not give a good quality fit to the broad-band spectrum, and two distinct particle sizes are needed, a larger particle size for the cloud deck and a smaller size particle for the higher altitude haze layer. For the condensate materials explored, we find cloud particle sizes between 0.52 and 1 μm are needed to fit the flat transmission spectrum between 0.52 to 3.6 μm and the small radius at 4.5 μm. Fitting for the data <0.52 μm, we find haze particle sizes of 0.02 to 0.1 μm can reproduce the near-UV to blue optical scattering slope, with condensates containing lower values for the imaginary part of the refraction index (like MgSiO₃ or Fe-poor Mg₂SiO₄ silicates) giving larger particle sizes.

Our results can also place important constraints on the vertical mixing rate in the atmosphere of WASP-31b. As discussed previously, a flat spectrum requires the presence of cloud or haze particles down to the 10-mbar level (based on the absence of the potassium line wings) or potentially even up to the 1-mbar level (based on the weakness of the water spectral feature). Fig. 13 shows the expected particle settling velocity for conditions relevant to WASP-31b as a function of particle radius and pressure, based on Stokes–Cunningham theory. The regime in the bottom portion of the plot corresponds to Stokes flow, where settling velocity is approximately independent of pressure. The regime at the top corresponds to that where the mean-free path becomes significant compared to the particle size, such that the continuum approximation breaks down and gas-kinetic effects are important; in this regime, the settling velocity for a particle of a given size scales approximately with the inverse of pressure. The figure shows that the settling velocity of 1-μm particles is 0.01 m s⁻¹ at 10 mbar and 0.1 m s⁻¹ at 1 mbar.

Figure 13.

Downward settling velocity (colour scale, m s⁻¹) for spherical atmospheric particles as a function of particle radius and pressure, assuming an H₂ atmosphere with a temperature of 1570 K and gravity of 4.56 m s⁻², appropriate to WASP-31b (see equations 3– 7 of Parmentier, Showman & Lian 2013). Stokes flow holds at high pressure and leads to settling velocities that depend on particle size but not pressure. Mean-free-path effects become important at low pressure, where settling velocity scales as approximately one over pressure. Here, we have assumed a particle density of 3000 kg m⁻³ (appropriate to silicate) and the gas viscosity appropriate to hydrogen from Rosner (2000); see Parmentier et al. (2013) for more details.

Open in new tabDownload slide

Given these estimates of settling velocity, it is interesting to estimate the effective atmospheric eddy diffusivities that would be necessary to loft such 1-μm particles to 1–10 mbar pressures, as necessary to cause a flat transmission spectrum. Parmentier et al. (2013) considered an idealized 1D diffusion model with downward transport due to particle settling being balanced by upward transport due to mixing by atmospheric motions (see equation 10 in Parmentier et al. 2013). An order-of-magnitude consideration of this equation suggests that, crudely, to keep particles of a given size lofted to a given altitude, one requires an eddy diffusivity of order K_zz ∼ Hw_settle, where H is the atmospheric scaleheight and w_settle is the particle settling velocity (e.g. as shown in Fig. 13). Given the ∼1000-km scaleheight expected on WASP-31b and the settling velocities quoted above, we estimate that the necessary atmospheric mixing corresponds to K_zz ∼ 10⁴ m² s⁻¹ to loft 1-μm particles to 10 mbar and K_zz ∼ 10⁵ m² s⁻¹ to loft them to 1 mbar. These estimates are consistent with the more detailed analysis presented in Parmentier et al. (2013).

Determining whether the atmospheric circulation on WASP-31b can generate sufficient mixing to induce K_zz ∼ 10⁴–10⁵ m² s⁻¹ is not trivial and will require detailed 3D numerical simulations for this planet, which is a task for the future. Nevertheless, Parmentier et al. (2013) made such estimates for HD 209458b, which has a similar effective temperature to WASP-31b. Because the incident starlight drives the vertical mixing, we might expect that planets with similar incident stellar irradiation, such as WASP-31b and HD 209458b, should exhibit similar profiles of the vertical mixing rate. Parmentier et al. (2013) showed that, interestingly, K_zz is not well represented by the product of the scaleheight and the rms vertical velocity of the atmospheric motions (as has sometimes been assumed in the literature). They fit diffusion models to their 3D simulations of the atmospheric circulation to obtain the vertical profile of the effective diffusivity: |$K_{\rm zz} \approx (5\times 10^4/\sqrt{P_{\rm bar}}) {\rm \,m^2\,s}^{-1}$|⁠, where P_bar is the pressure in bar. This equation suggests that the effective eddy diffusivity should be ∼5 × 10⁵ m² s⁻¹ at 10 mbar and ∼2 × 10⁶ m² s⁻¹ at 1 mbar. These values exceed the estimates we made above for the diffusivity required to mix ∼1-μm particles to pressures of 1–10 mbar. Thus, it appears physically plausible that the circulation will indeed be sufficiently strong to mix particles upward to pressures of 1–10 mbar, as necessary to explain our flattened spectrum.

5 CONCLUSION

We have presented HST and Spitzer observations totalling six transits of WASP-31b, which cover the full-optical to infrared spectral regions. The resulting exoplanet transmission spectrum is unique, showing several atmospheric features not yet seen in conjunction including potassium line cores and no detectable sodium, a flat optical to near-IR slope indicative of a cloud deck, and Rayleigh scattering at short wavelengths. The potassium feature is strong (detected at 4.2σ confidence) though lacks significant pressure broadened wings, which we find limits the pressures probed to 10 mbar or less. Weak sodium and strong potassium lines indicate the peculiar possibility of a sub-solar [Na/K] abundance ratio below one, which is difficult to explain though may point towards Na condensation or large element-to-element metallically variations from formation/accretion processes.

The transmission spectrum is flat longwards of 0.52 μm, indicative of a high altitude broad-band cloud deck which covers the expected optical Na and near-IR H₂O features. However, despite containing a substantial cloud(s) the broad-band transmission spectrum is not flat, and a 4.2σ significant slope is observed at wavelengths shortwards of 0.52 μm, consistent with Rayleigh scattering by small sub-micron size haze particles.

WASP-31b joins a growing group of hot-Jupiter exoplanets with evidence for clouds and hazes. Clouds and hazes are now seen to be important atmospheric features in hot Jupiters across a wide range of effective temperatures, including the hottest planets like WASP-12b and much cooler planets like HD 189733b. Moreover, at temperatures near 1570 K, WASP-31b has an equilibrium temperature close to the mean for known hot Jupiters, implying atmospheric clouds and haze signatures are not confined to the extremes of hot-Jupiter parameter space.

We highlight the importance of wide wavelength coverage for exoplanet transmission spectra, as short-wavelength data can be highly sensitive to the cloud properties and the blue-optical data were necessary in this case to detect aerosol scattering features.

This work is based on observations with the NASA/ESA HST, obtained at the Space Telescope Science Institute (STScI) operated by AURA, Inc. 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. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no. 336792. We thank the anonymous referee for their comments. DS, FP, and NN acknowledge support from STFC consolidated grant ST/J0016/1. Support for this work was provided by NASA through grants under the HST-GO-12473 programme from the STScI. PW and HW acknowledge support from the UK Science and Technology Facilities Council (STFC). ALE and AVM acknowledge support from the French Agence Nationale de la Recherche (ANR), under programme ANR-12-BS05-0012 ‘Exo-Atmos’.

REFERENCES