HUBBLE SPACE TELESCOPE NEAR-IR TRANSMISSION SPECTROSCOPY OF THE SUPER-EARTH HD 97658B

The Kepler mission has resulted in the discovery of more than 3000 transiting planets and planet candidates (Batalha et al. 2013; Burke et al. 2014) to date, with a majority of the sample consisting of sub-Neptune-sized planets. Analyses of this set of Kepler planet candidates indicate that planets with radii intermediate between that of Neptune and the Earth appear to be the most common kind of extrasolar planet orbiting nearby FGK stars, with a peak at radii between 2–3 × that of the Earth (Howard et al. 2012; Fressin et al. 2013; Petigura et al. 2013). Planets in this size range are typically referred to as "super-Earths," although they could potentially form with a broad range of compositions including primarily rocky with a thin atmosphere (true "super-Earths"), a rocky or icy core surrounded by a thick hydrogen atmosphere ("mini-Neptunes"), or water-dominated with a thick steam atmosphere ("water worlds"). Many of these super-Earths are found in close-in, tightly packed multiple planet systems (e.g., Lissauer et al. 2011; Fabrycky et al. 2012; Steffen et al. 2013), and there is an ongoing debate as to whether these systems formed in place or migrated in from more distant orbits (Hansen & Murray 2012; Chiang & Laughlin 2013; Raymond & Cossou 2014).

Detailed studies of super-Earth compositions offer important clues on their origins: presumably water-rich planets must have formed beyond the ice line, while in-situ formation models predict primarily rocky compositions with relatively water-poor hydrogen-dominated atmospheres (e.g., Raymond et al. 2008). By combining radius measurements from Kepler with mass estimates obtained using either the radial velocity or transit timing techniques, it is possible to constrain the average densities and corresponding bulk compositions of the super-Earths in the Kepler sample (e.g., Lithwick et al. 2012; Hadden & Lithwick 2014; Weiss & Marcy 2014; Marcy et al. 2014). These observations indicate that the super-Earths in the Kepler sample display a broad range of average densities, with a transition toward denser, primarily rocky compositions below 1.5–2 Earth radii (Weiss & Marcy 2014; Marcy et al. 2014). For the larger, lower-density super-Earths in the Kepler sample it is possible to match their measured densities with a broad range of compositions, including both water-rich and water-poor scenarios, simply by varying the amount of hydrogen in their atmospheres (Seager et al. 2007; Valencia et al. 2007, 2013; Rogers & Seager 2010a; Zeng & Sasselov 2013). Although it can be argued that some compositions are unlikely based on models of planet formation and atmospheric mass loss for close-in planets, this still leaves a wide range of plausible models (e.g., Rogers & Seager 2010b; Rogers et al. 2011; Nettelmann et al. 2011; Heng & Kopparla 2012; Lopez et al. 2012; Lopez & Fortney 2013; Fortney et al. 2013).

Measurements of the transmission spectra of super-Earths allow us to directly estimate the mean molecular weight of the planet's atmosphere (e.g., Miller-Ricci et al. 2008; Fortney et al. 2013), which in turn provides improved constraints on its interior composition. Planets with cloud-free, hydrogen dominated atmospheres will have relatively large-scale heights and correspondingly strong absorption features during transit, while planets with hydrogen-poor atmospheres will have relatively small-scale heights and weak absorption. Unfortunately, the majority of the super-Earths detected by Kepler orbit faint stars (m_K > 10), making it difficult to accurately measure their transmission spectra using existing facilities. There are currently three super-Earths known to transit relatively bright stars, including GJ 1214b (m_K = 8.8; Charbonneau et al. 2009), 55 Cancri e (m_K = 4.0; Winn et al. 2011; Demory et al. 2011), and HD 97658b (m_K = 5.7; Dragomir et al. 2013).

GJ 1214b is currently the only one of these three systems with a well-characterized transmission spectrum, which is flat and featureless across a broad range of wavelengths (e.g., Bean et al. 2010, 2011; Désert et al. 2011; Berta et al. 2012; Kreidberg et al. 2014). Although initial observations were consistent with either a cloud-free, hydrogen-poor atmosphere or a hydrogen-rich atmosphere with a high cloud deck (e.g., Bean et al. 2010), the most recent near-IR observations obtained by Kreidberg et al. (2014) using the Wide Field Camera 3 (WFC3) instrument on the Hubble Space Telescope (HST) are precise enough to rule out cloud-free models over a broad range of atmospheric metallicities (also see Benneke & Seager 2013). The presence of a high altitude cloud or haze layer means that for this planet, at least, transmission spectroscopy provides relatively weak constraints on the mean molecular weight of its atmosphere and, by extension, on its interior composition.

In this paper, we present HST WFC3 near-infrared transmission spectroscopy of the transiting super-Earth HD 97658b. This planet was first detected using the radial velocity technique (Howard et al. 2011), and later found to transit using MOST photometry (Dragomir et al. 2013). It has a mass of 7.9  ±  0.7 M_⊕, a radius of 2.3  ±  0.2 R_⊕, and an average density of 3.4  ±  0.9 g cm⁻¹ (Dragomir et al. 2013), making it modestly denser and more massive than GJ 1214b. HD 97658b orbits its K star primary with a period of 9.5 days, and has a predicted zero-albedo temperature between 700–1000 K depending on the efficiency of energy transport to the planet's night side. If this planet has the same atmospheric composition as GJ 1214b, its modestly higher atmospheric temperature might prevent the formation of the cloud layer detected in GJ 1214b's atmosphere (Morley et al. 2013). Our observations are obtained at wavelengths between 1.2–1.6 μm, and are primarily sensitive to the presence or absence of the water absorption band located at 1.4 μm; this feature has now been robustly detected in the atmospheres of several hot Jupiters (e.g., Deming et al. 2013; Wakeford et al. 2013; Mandell et al. 2013; McCullough et al. 2014), and is expected to be present in super-Earths as well over a broad range of atmosphere compositions (e.g., Benneke & Seager 2012; Kreidberg et al. 2014; Hu & Seager 2014). We discuss our observations in Section 2 and the implications for the properties of this planet's atmosphere in Section 3.

We observed one transit of HD 97658b on UT 2013 December 19 and another on UT 2014 January 7 (GO 130501, PI Knutson) using the G141 grism on the HST WFC3 instrument, which provides low-resolution spectroscopy at wavelengths between 1.1–1.7 μm. Each observation consists of five HST orbits with a total duration of approximately seven hours per visit; our target was visible for approximately half of each 96 minute HST orbit. Scheduling constraints resulted in a slightly shorter first orbit during the January visit, which consisted of 203 spectra instead of the 206 spectra obtained during the December visit. During the January visit, we also utilized a series of short exposures at the end of each orbit in order to force a buffer dump, which avoided the mid-orbit gaps in coverage visible in the light curves for the first visit in Figure 1. We do not include these short exposures in our analysis, as they were taken in imaging rather than spectroscopic mode.

Zoom In Zoom Out Reset image size

Our spectra were obtained using the 256 × 256 pixel subarray and SPARS10 mode with four samples, giving a total integration time of 14.97 s in each image. Our target is one of the brightest transiting planet host stars, and we therefore utilized the new scanning mode (e.g., McCullough & MacKenty 2012; Deming et al. 2013; Kreidberg et al. 2014; Knutson et al. 2014) with a scan rate of 14 s⁻¹ in order to achieve a higher observing efficiency while remaining well below saturation. This results in a scanned spectral image that fills most of the subarray image, with a fractional coverage similar to that of the HD 209458 observations from Deming et al. (2013). We also alternated between forward and reverse scan directions in order to further reduce overheads; this approach was previously used by Kreidberg et al. (2014) for GJ 1214b. The resulting spectra have peak counts around 40,000 electrons pixel⁻¹, which is high enough to cause a modestly steep asymptotic ramp in the measured fluxes within individual orbits (see Wilkins et al. 2014, for a discussion of the relationship between the ramp slope and total flux).

We reduced the data using two independent methods, which allow us to test whether or not the transmission spectrum we obtain from our analysis is sensitive to our choice of analysis technique. We discuss each approach separately below.

2.1. Spectral Template Fitting Method

In this approach, we utilized the spectral template method first presented in Deming et al. (2013), which we describe in detail in Knutson et al. (2014) and summarize here. We begin with the raw sample up the ramp images from the ima.fits files, which are bias and dark subtracted, and apply our own flat-fielding and wavelength calibrations based on the standard WFC3 pipeline as discussed in Wilkins et al. (2014). We then subtract successive non-destructive readout pairs in order to create a series of difference images, trim the region around the spectral scan in each image, and sum the resulting images to create a composite containing the full scan. The benefit of this approach is that it minimizes the contributions of pixels that are not directly illuminated by the star in a given readout time step (see Deming et al. 2013, for more discussion on this point). We estimate that the sky background in our final composite images is 0.18%–0.26% of the total flux measured in each wavelength element (see Figure 2), and subtract this estimated sky background from each image.

Zoom In Zoom Out Reset image size

Once we have created a stacked image from each set of readouts, we apply a filter to correct for bad pixels and cosmic ray hits as described in Knutson et al. (2014). We then sum the spectrum along the y axis in order to extract a one-dimensional spectrum, using an aperture that extends fifteen pixels above and below the edges of the spectrum (200 pixels in total) in order to include the extended wings of the point-spread function (see Figure 3 for a representative example). Because the WFC3 spectra are undersampled (Deming et al. 2013), we convolve each of our one-dimensional spectra with a Gaussian function with a width (FWHM) of four pixels in order to mitigate effects related to the shifting position of the spectrum on the detector. Although this modestly degrades the spectral resolution of our data, we later bin our transmission spectrum by an equivalent amount in order to increase the signal-to-noise ratio. We calculate the MJD mid-exposure time corresponding to each image using the information from the flt.fits image headers, and convert these times to the BJD_TDB time standard following the methods of Eastman et al. (2010).

Zoom In Zoom Out Reset image size

Our next step is to create a spectral template by averaging the ten spectra immediately before and after the transit. We then fit this template spectrum to each individual spectrum in our time series, allowing the relative position and amplitude of the template spectrum to vary as free parameters. We find no difference in our results if we create a separate spectral template for the forward and reverse scan images, and we therefore use the same template for all images. The resulting series of best-fit amplitudes are plotted in Figure 1, and are identical to the white-light curves obtained by summing the fluxes across all wavelengths. We find that there is a small flux offset between the light curves for the two scan directions, which we remove by dividing each time series by its median flux value. McCullough & MacKenty (2012) suggest that this offset is a consequence of the order in which columns are read out by the detector (the "up-stream/down-stream effect"); when the scan moves in the same direction as the readout then the effective integration time will be slightly longer than in the case of a reverse scan. The forward and reverse spectra in our images also occupy slightly different positions on the array, which might also contribute to this offset.

We subtract the best-fit spectral template from each individual spectrum in order to create a differential time series for each individual wavelength element. The benefit to this approach is that it effectively removes all common-mode detector effects from the differential light curve. We find that the scatter in our light curves for individual wavelength elements is within 5% of the photon noise limit in all cases, indicating that there is minimal color-dependence in the systematic noise sources.

2.1.1. White-light Transit Fits

We take the two raw transit light curves plotted in Figure 1 and fit them each with a transit light curve, a linear function of time, and an exponential ramp in orbital phase:

$$\begin{equation} F(t) = c_1(1+ c_2t+c_3e^{-p/c_4})F_{{\rm transit}}(t), \end{equation} \tag{ 1 }$$

where c₁ − c₄ are free parameters in the fit, t is the time from the measured transit center in days, and p is the time in days from either the first observation in each spacecraft orbit or the last observation before a mid-orbit buffer dump. We use this definition because we find that the ramp after a mid-orbit buffer dump is usually not as steep as the initial ramp at the start of the orbit, and using this definition provided a good fit to the observed behavior by reducing the amplitude of the mid-orbit ramp model. The latter definition is only used for the December transit observation, which includes two mid-orbit buffer dumps that cause the ramp effect to reset. We allow the planet–star radius ratio R_p/R_⋆ and center of transit time to vary as free parameters for each individual transit. We find that the behavior of the instrumental systematics is slightly different for the white-light curves calculated from the forward and reverse scan directions (see Figure 1). We therefore allow the linear function of time and exponential ramp functions to vary independently for each scan direction, while keeping the same R_p/R_⋆ and center of transit time.

We calculate the transit light curve F_transit following Mandel & Agol (2002) where we initially fixed the orbital inclination i and the ratio of the semi-major axis to the stellar radius a/R_⋆ equal to their best-fit values from Van Grootel et al. (2014). We found that the uncertainty on our reported transit time for the first visit, which spanned the transit ingress, was particularly sensitive to our choice of values for these parameters. Our transit timing uncertainties for this visit ranged between 30 s up to 5 minutes for the cases where i and a/R_⋆ were either fixed or allowed to vary freely in our fits. We addressed this by acquiring the normalized 4.5 μm Spitzer transit photometry from Van Grootel et al. (2014) and fitting this light curve simultaneously with our white-light HST photometry assuming a single global value for i and for a/R_⋆. This ensured that our final reported transit times correctly accounted for the added uncertainties contributed by these two parameters. We allowed the Spitzer planet–star radius ratio and center of transit time to vary independently in our fits, and found that our results for these two parameters were consistent with those reported by Van Grootel et al. (2014). In our final fits, the uncertainty on the transit time for the first visit increased from 30 s to one minute as compared to the fixed i and a/R_⋆ case and the best-fit transit time was 3.8 minutes earlier. The second visit does not include observations during ingress or egress, and has a correspondingly large uncertainty on the best-fit transit time. We find that for this visit, allowing i and a/R_⋆ to vary has a negligible effect on the reported uncertainties in the best-fit transit time.

We calculate our initial predicted transit times for the HST observations using the ephemeris from Van Grootel et al. (2014), and then derive an updated estimate for the orbital period and center of transit time based on our new observations. We then repeat our transit fits using this updated orbital period. We calculate four-parameter nonlinear limb-darkening parameters (Claret 2000) for our transit light curves, with the relative intensity at each position given as

$$\begin{equation} \frac{I(\mu)}{I(1)} = 1-\sum _{k=1}^{4}a_k(1-\mu ^{k/2}), \end{equation} \tag{ 2 }$$

where I(1) is the specific intensity at the center of the stellar disk, a_k is the kth limb-darkening coefficient, and μ = cos(θ) where θ is the angle between the line of sight and the location of the emerging flux. We calculate our limb-darkening coefficients using a PHOENIX stellar atmosphere model (Allard et al. 2012), where we take the flux-weighted average of the theoretical stellar intensity profile across the band and then fit for the limb-darkening coefficients. We use the best-fit stellar parameters from Van Grootel et al. (2014), who find an effective temperature of 5170  ±  50 K, a surface gravity of 4.58  ±  0.05, and a metallicity of −0.23  ±  0.03. These values are also consistent with a recent analysis by Mortier et al. (2013), although this study prefers a lower metallicity of −0.35  ±  0.02 and does not take into account the constraints on the stellar density from the transit light curve. We also tried models with effective temperatures ranging between 5120–5220 K and a metallicity of −0.35, but found that these had a negligible effect (2 ppm or less) on our resulting transmission spectrum. The inclusion of limb-darkening in our fits creates a small offset in the average transit depth as compared to fits without limb-darkening, and has a negligible effect on the slope of the resulting transmission spectrum across the bandpass. We find that the uncertainty in the stellar effective temperature contributes a systematic error of 1 ppm to our estimate of the differential transit depths.

We estimate the uncertainties in our fitted white-light parameters using a residual permutation ("prayer bead") method, the covariance matrix from our Levenberg–Marquart minimization, and a Markov Chain Monte Carlo (MCMC) analysis with 10⁶ steps. In our initial fits with fixed i and a/R_⋆, we set the uncertainties on individual points in our HST white-light time series equal to the standard deviation of the residuals from our best-fit solution. We find that the residual permutation technique results in uncertainties that are two to three times larger than the corresponding values from both the covariance matrix and MCMC analysis. This is expected, as the noise in our white-light curves is dominated by time-correlated instrument effects, while the MCMC and covariance methods implicitly assume random Gaussian-distributed noise (Carter & Winn 2009). We plot the two-dimensional probability distributions from the MCMC analysis and confirm that there are no significant correlations between any of the fit parameters in this case. For the case where we fit the Spitzer and HST transits simultaneously while allowing i and a/R_⋆ to vary as free parameters, we increase the per-point uncertainties on the HST transit light curves by a factor of two in order to more accurately reflect the uncertainties contributed by the time-correlated component of the noise. We selected this scaling factor by requiring that the uncertainties on R_p/R_⋆ derived from the covariance matrix be equal to the uncertainties from the prayer bead method. This step is critical for the joint fits, because otherwise the HST data have a disproportionate influence on the preferred values of i and a/R_⋆. We find that in this version of the fit the values of the inclination and a/R* show a high degree of correlation, which is a well-known property of this particular parameterization of the transit light curve shape. None of the transit parameters in either version of the fits are correlated with the parameters used to describe the time-varying instrumental signal.

We give our best-fit transit parameters for the joint HST and Spitzer fits in Table 1, and plot the normalized transit light curves after dividing out our best fit for detector effects in Figure 4. The white-light curve residuals from UT 2013 December 19 and UT 2014 January 7 have an rms of [73, 68] ppm and [63, 57] ppm, respectively, where we calculate this number separately for data from the forward and reverse scan directions. We show the residuals from these fits in Figure 5, and the offset between the observed and calculated transit times from our updated orbital ephemeris in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Transit Parameters from Joint HST and Spitzer Fits

Parameter	Value
Global values
i (°)	89.85 ± 0.48a
a/R⋆	26.24 ± 1.21
P (days)b	9.489264 ± 0.000064
T0 (BJDTDB)b	2456665.46415 ± 0.00078
UT 2013 Aug 10 Spitzer transitc
Rp/R⋆	0.02778 ± 0.00073
Tc (BJDTDB)	2456523.12537 ± 0.00049
UT 2013 Dec 19 HST transit
Rp/R⋆	0.03012 ± 0.00087
Tc (BJDTDB)	2456646.48556 ± 0.00071
UT 2014 Jan 7 HST transit
Rp/R⋆	0.03111 ± 0.00080
Tc (BJDTDB)	2456665.4584 ± 0.0056

Notes. ^aWe limited the value of the inclination to be less than or equal to 90° in our fits. We find a 1σ lower limit on the inclination of 8965, and a 2σ limit of 8879. ^bWe calculate our updated ephemeris using the published transit times from Dragomir et al. (2013) and Van Grootel et al. (2014) and our two new transit measurements. T₀ is the zero epoch transit center time from our best-fit ephemeris, and T_c is the measured transit center time from a given observation. ^cWe find that for the Spitzer data, the uncertainties from the covariance matrix are larger than those from the prayer bead analysis, and we report the covariance errors here.

Download table as:  ASCIITypeset image

2.1.2. Differential Transit Fits

We estimate the planet's wavelength-dependent transit depth by fitting the differential time series for each individual wavelength element. In this case, we fit the differential time series for all five orbits rather than the four used in the white-light fits, as we find that the differential light curves do not show any detectible systematic trends during the first orbit. The light curves for the forward and reverse scans are offset by a constant flux value, so we divide each light curve by its median value before carrying out our fit. We fit each light curve with a linear function of time and a transit function calculated as the difference between the transit light curve in that band and the white-light transit curve, where we fix the values of i and a/R_⋆ to their best-fit values from the white-light analysis. We also tried fits where we allowed an independent linear function of time for each scan direction, but found that this gave a transmission spectrum that was indistinguishable from the case where we assumed the same linear function for both scan directions. We use the best-fit transit time from the white-light fits, and fit for the planet–star radius ratio R_p/R_⋆ corresponding to each wavelength element. The best-fit transit depths reported here are simply the square of these values.

We calculate the appropriate four-parameter nonlinear limb-darkening coefficients for each wavelength element using the same PHOENIX model used for the white-light fits (see Table 2), where we convolve the model spectrum at each position on the star with the same Gaussian function used on our data before fitting for the limb-darkening coefficients. We find that varying the stellar effective temperature by ±50 K changes the resulting transit depths by 1 ppm, which we include in our analysis as a separate systematic error term.

Table 2. Binned Four Parameter Nonlinear Limb-darkening Coefficients^a

Wavelengthb	c1	c2	c3	c4
1.133	0.6247	−0.4071	0.6890	−0.3099
1.152	0.6489	−0.4619	0.7545	−0.3395
1.171	0.6710	−0.5019	0.7856	−0.3530
1.190	0.6925	−0.5644	0.8445	−0.3743
1.208	0.7275	−0.6442	0.9250	−0.4065
1.227	0.7631	−0.7246	1.0242	−0.4507
1.246	0.8020	−0.8191	1.1222	−0.4894
1.265	0.8576	−0.9453	1.2552	−0.5437
1.284	0.4703	0.1247	0.0695	−0.0947
1.303	0.5182	−0.01457	0.2336	−0.1563
1.321	0.5640	−0.1255	0.3410	−0.1967
1.340	0.6366	−0.2849	0.5109	−0.2661
1.359	0.7125	−0.4597	0.6892	−0.3356
1.378	0.8012	−0.6639	0.8937	−0.4138
1.397	0.8842	−0.8546	1.0907	−0.4911
1.416	0.8901	−0.8644	1.1018	−0.4975
1.434c	1.0237	−1.1944	1.4386	−0.6244
1.453	1.1298	−1.4436	1.6910	−0.7202
1.472	1.2482	−1.7152	1.9657	−0.8248
1.491c	0.4779	0.4024	−0.4201	0.1113
1.510	0.5016	0.3492	−0.3886	0.1055
1.529	0.5557	0.3381	−0.4372	0.1307
1.547	0.6129	0.2520	−0.3962	0.1249
1.566	0.6479	0.1523	−0.3120	0.1000
1.585	0.6325	0.1655	−0.3350	0.1121
1.604	0.7095	0.0517	−0.2727	0.0993
1.623	0.7208	0.0239	−0.2472	0.0905
1.642	0.7129	0.0421	−0.2933	0.1144
4.5d	0.8713	−1.4634	1.4955	−0.5581

Notes. ^aThe spectra extracted from the UT 2013 December 19 and UT 2014 January 7 visits have slightly different wavelength solutions, and we therefore created custom limb darkening coefficients for each visit. A full table of both the binned and unbinned coefficients for both visits is available upon request. ^bThese binned coefficients are only used in the white-light residual fitting analysis. In the spectral template fitting analysis we fit the light curves for individual wavelength elements and use limb-darkening coefficients calculated appropriately for this resolution (4× higher than shown here). ^cThere is a well-known degeneracy in the values of the four-parameter nonlinear limb-darkening coefficients; although the values of the coefficients do not vary smoothly, we find that the corresponding stellar intensity profiles are continuous across all wavelengths considered here. ^dLimb-darkening coefficients calculated for the 4.5 μm Spitzer IRAC bandpass, which we use to fit the transit light curve from Van Grootel et al. (2014).

Download table as:  ASCIITypeset image

We estimate uncertainties in the best-fit planet–star radius ratios for each wavelength element using both a residual permutation method and the covariance matrix from our Levenberg–Marquart minimization, where we set the uncertainties on individual data points in a given wavelength channel equal to the standard deviation of the residuals in that channel after the best-fit model has been subtracted. We find that the residual permutation errors for both visits are on average higher than the corresponding covariance matrix errors, suggesting the presence of time-correlated noise in our light curves. However, the residual permutation errors also display a greater variation in size from one wavelength channel to the next, likely as a result of sampling error due to the limited number of data points in each light curve. We therefore take a conservative approach and select the larger of the two error estimates in each wavelength channel as our final error estimate, which allows us to account for time-correlated noise in our data while avoiding possible under-estimation of uncertainties due to the limited number of measurements available for the residual permutation analysis.

Next, we bin the resulting wavelength-dependent radius values by a factor of four to create the transmission spectrum shown in Figure 7, where we average the estimated radii and corresponding errors in each of the twenty eight bins. Because we have smoothed our spectra in wavelength space prior to fitting for the transit depths in each individual channel, the transit depths in adjacent wavelength channels are correlated with each other and the errors on the binned transit depths do not decrease as the square root of the number of wavelength channels in the bin. We estimate the effect of this smoothing on the binned errors by creating a simulated data set with the same time sampling and number of wavelength channels as our original data. We add noise to these simulated data and fit them in exactly the same way as our original data in order to derive a binned transmission spectrum. We then repeat this analysis with the same simulated data, but including an additional step where we smooth the simulated "spectra" using the same Gaussian function that we apply to our actual spectra. When creating our binned transmission spectrum, we find that we must average the four individual error estimates from the smoothed version of the data and multiply this average by a factor of 1.22 in order to obtain the same binned uncertainties as in the unsmoothed version. We adopt this same scaling when calculating the errors on the binned transit depth estimates from our real data.

Zoom In Zoom Out Reset image size

We take the error-weighted average of the two individual transmission spectra to create the combined spectrum shown in Figure 8 and in Table 3. We present two separate estimates for the measurement uncertainties in Table 3, including one calculated using the covariance matrix errors alone and the other taking the larger of the two error estimates as previously described. As an additional consistency check, we also carried out a version of our analysis in which we allowed the light curves derived from the forward and reverse scans to have different planet–star radius ratios. We found that we obtained consistent radius estimates from both scan directions, indicating that the different behaviors visible in the white-light curves are effectively removed in the differential time series by our template fitting technique.

Zoom In Zoom Out Reset image size

Table 3. Wavelength Dependent Transit Depths from Template Fitting Method

Wavelength	Depth	White Noise Errora	Total Errora
(μm)	(ppm)	(ppm)	(ppm)
1.132	912	20	27
1.151	896	20	34
1.170	916	20	22
1.188	917	19	19
1.207	913	20	20
1.226	897	19	20
1.245	898	18	19
1.264	934	18	18
1.283	948	18	19
1.301	952	19	20
1.320	952	19	21
1.339	936	18	25
1.358	933	18	25
1.377	925	19	29
1.396	914	19	20
1.415	922	19	32
1.433	981	20	27
1.452	946	19	25
1.471	956	18	22
1.490	945	18	23
1.509	964	20	27
1.528	965	20	21
1.546	947	19	20
1.565	946	21	28
1.584	963	21	23
1.603	973	20	22
1.622	968	20	23
1.641	953	24	30

Notes. ^aWhite noise measurement errors are estimated using the covariance matrix, which implicitly assumes that individual measurement errors are Gaussian distributed and uncorrelated. The total measurement errors are calculated by comparing the covariance errors to a residual permutation method that better accounts for time-correlated noise and taking the larger of the two in each wavelength channel. Uncertainty in the limb-darkening models from the stellar effective temperature contributes an additional error of 1 ppm in each bin, which we do not include here.

Download table as:  ASCIITypeset image

2.2. White- light Residual Fitting Method

We used the WFC3 instrument on the HST to observe transits of the super-Earth HD 97658b at wavelengths between 1.1–1.7 μm. Our white-light transit fits produce consistent estimates of the transit times and planet–star radius ratios between the two visits. Our errors for these parameters are dominated by systematic noise, as illustrated in Figures 1–3. We combine our measured transit times with previously published values from MOST and Spitzer and show that they are consistent with a linear ephemeris. We also derive an improved estimate for the planet's orbital period with an uncertainty that is a factor of 30 smaller than in Van Grootel et al. (2014). Our average radius ratio is 2.5σ larger than the Spitzer 4.5 μm measurement from Van Grootel et al. (2014), and is in good agreement with the MOST visible-light value from Dragomir et al. (2013). It is unlikely that this offset is due to stellar activity as this star has a log($R^{\prime }_{\scriptsize {HK}}$)value between −4.95 and −5.00 (Howard et al. 2011) and varied by less than 0.2% in brightness over a single 1.5 day visible-light observation by MOST (Dragomir et al. 2013).

We obtain nearly identical versions of the transmission spectrum using the spectral template fitting technique (Deming et al. 2013; Knutson et al. 2014) and the white-light residual method (Kreidberg et al. 2014), as shown in Figure 8. Our median uncertainties on the differential wavelength-dependent transit depth are 23 ppm from the spectral template fitting and 26 ppm from the white-light residual method (see Figure 9), making these observations the most precise measurement of a planetary transmission spectrum that we are aware of with WFC3 to date. Both versions of our spectrum display a discontinuity around 1.42 μm, which may be related to a slight upward slope visible toward redder wavelengths. These features could be instrumental in nature, as the columns that make up the wavelength bins closest to the discontinuity are located on a boundary between quadrants on the WFC3 array. The statistical significance of these features is marginal, and we therefore do not consider them in our subsequent comparison to atmosphere models for this planet.

Zoom In Zoom Out Reset image size

In Figure 10, we compare our measured transmission spectrum from the spectral template fitting technique to several representative atmosphere models for HD 97658b, which are calculated following Kempton et al. (2012). The effective temperatures (720–730 K, depending on the model) and corresponding pressure-temperature profiles of these models are calculated assuming full redistribution of energy to the planet's night side and an albedo which varies according to composition. We consider both cloud-free models with solar and 50 × solar metallicities, as well as a pure water model and a solar metallicity model with a high cloud deck located at one mbar. We calculate the significance with which our data can rule out a given model using the equation from Gregory (2005):

$$\begin{equation} {\rm Significance} = \frac{\chi ^2-\nu }{\sqrt{2\nu }}, \end{equation} \tag{ 3 }$$

where χ² is calculated by comparing our averaged transmission spectrum to each model and ν is the number of degrees of freedom in the fit (27 in our case, as there are 28 points and we normalize the models to match the average measured transit depth of our data). This metric assumes that our measurement errors are Gaussian and uncorrelated from one wavelength bin to the next; although this is almost certainly untrue at some level, it represents a reasonable starting point for comparing different models. Following this approach we find that our measured transmission spectrum is inconsistent with the solar and 50 × solar cloud-free models at the 10σ and 9σ levels, respectively. It is equally well described by the water-dominated (0.6σ) model and the solar metallicity model with optically thick clouds at a pressure of one mbar (0.6σ), as well as a flat line at the average transit depth across the band (0.4σ). We find that a solar metallicity model with clouds at 1 mbar is also consistent with the data (1.8σ), indicating that the clouds could be located slightly deeper in the atmosphere. We note that there are any number of high metallicity atmosphere models that could provide a fit comparable to that of the pure water model; all our data appear to require is either a relatively metal-rich atmosphere with a correspondingly small-scale height, or the presence of a high cloud deck that obscures the expected water absorption feature in a hydrogen-dominated atmosphere. We constrain the maximum hydrogen content of the atmosphere in the first scenario by considering a series of cloud-free models with varying number fractions of molecular hydrogen and water, and find that in this scenario the atmosphere has to be at least 20% water by number in order to be consistent with our data at the 3σ level. We list the χ² values and level of disagreement for all models in Table 5.

Zoom In Zoom Out Reset image size

The conclusion that HD 97658b's transmission spectrum appears to be flat at the precision of our data places it in the same category as both the super-Earth GJ 1214b (Kreidberg et al. 2014) and the Neptune-mass GJ 436b (Knutson et al. 2014). As with these two planets, a more precise measurement of HD 97658b's transmission spectrum will eventually allow us to distinguish between high clouds and a cloud-free, metal-rich atmosphere. Our constraints on the atmospheric scale height in the cloud-free scenario are relatively weak compared to those obtained for GJ 1214b and GJ 436b, despite the fact that we achieve smaller errors (20 ppm versus 30 ppm) in our measurement of the differential transmission spectrum. This is primarily because HD 97658b has a smaller planet–star radius ratio than either of these systems, and the predicted amplitude of the transmission spectrum is correspondingly small. Fortunately, it also orbits a brighter star than either GJ 1214b or GJ 436b, making it possible to achieve high precision transit measurements with relatively few observations. Unfortunately, this makes ground-based observations particularly challenging as the nearest comparison star with a comparable brightness is located more than 40' away. For space telescopes such as Hubble and Spitzer, achieving the precision required to study this planet in detail will mean pushing the systematic noise floor to unprecedentedly low levels. There is every reason to believe that this level of performance should be achievable, and given the unique nature of this planet it is likely that this will be put to the test in the near future.

We are grateful to STScI Director Matt Mountain and the Director's Discretionary Time program for making it possible to obtain these observations prior to the next regular HST proposal cycle. We also thank Beth Perriello and the rest of the HST scheduling team, as they worked overtime in the days leading up to the Christmas holiday in order to make sure that we obtained the best possible data from this program. We also wish to acknowledge Valerie Van Grootel, who was kind enough to share an advance copy of her paper prior to publication. H.K. acknowledges support from NASA through grant GO-13501. L.K. received funding for this work from the National Science Foundation in the form of a Graduate Research Fellowship. J.L.B. acknowledges support from the Alfred P. Sloan Foundation and NASA through grants NNX13AJ16G, GO-13021, and GO-13467. D.H. has been supported through the European Research Council advanced grant PEPS awarded to Gilles Chabrier under the European Community's Seventh Framework Program (FP7/2007-2013 grant agreement No. 247060).