arXiv:1403.4602v2 [astro-ph.EP] 27 Aug 2014
Draft version August 29, 2014
Preprint typeset using L
A
T
E
X style emulateapj v. 08/22/09
HUBBLE SPACE TELESCOPE
NEAR-IR TRANSMISSION SPECTROSCOPY OF THE SUPER-EARTH HD
97658B
Heather A. Knutson
1,2
, Diana Dragomir
3,4
, Laura Kreidberg
5
, Eliza M.-R. Kempton
6
, P. R. McCullough
7
,
Jonathan J. Fortney
8
, Jacob L. Bean
5
, Michael Gillon
9
, Derek Homeier
10
, Andrew W. Howard
11
Draft version August 29, 2014
ABSTRACT
Recent results from the
Kepler
mission indicate that super-Earths (planets with masses between
1
−
10 times that of the Earth) are the most common kind of planet arou
nd nearby Sun-like stars.
These planets have no direct solar system analogue, and are curre
ntly one of the least well-understood
classes of extrasolar planets. Many super-Earths have average
densities that are consistent with a
broad range of bulk compositions, including both water-dominated w
orlds and rocky planets covered
by a thick hydrogen and helium atmosphere. Measurements of the t
ransmission spectra of these
planets offer the opportunity to resolve this degeneracy by direct
ly constraining the scale heights and
corresponding mean molecular weights of their atmospheres. We pr
esent
Hubble Space Telescope
near-
infrared spectroscopy of two transits of the newly discovered tr
ansiting super-Earth HD 97658b. We
use the Wide Field Camera 3’s scanning mode to measure the wavelengt
h-dependent transit depth
in thirty individual bandpasses. Our averaged differential transmis
sion spectrum has a median 1
σ
uncertainty of 23 ppm in individual bins, making this the most precise o
bservation of an exoplanetary
transmission spectrum obtained with WFC3 to date. Our data are inc
onsistent with a cloud-free
solar metallicity atmosphere at the 10
σ
level. They are consistent at the 0
.
4
σ
level with a flat line
model, as well as effectively flat models corresponding to a metal-rich
atmosphere or a solar metallicity
atmosphere with a cloud or haze layer located at pressures of 10 mb
ar or higher.
Subject headings:
binaries: eclipsing — planetary systems — techniques: spectroscop
y
1.
INTRODUCTION
The
Kepler
mission has resulted in the discovery of
more than three thousand 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 inter-
mediate 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 composi-
tions including primarily rocky with a thin atmosphere
(true “super-Earths”), a rocky or icy core surrounded
1
Division of Geological and Planetary Sciences, California
In-
stitute of Technology, Pasadena, CA 91125, USA
2
hknutson@caltech.edu
3
Las Cumbres Observatory Global Telescope Network, Goleta,
CA 93117, USA
4
Department of Physics, Broida Hall, UC Santa Barbara, CA,
USA
5
Department of Astronomy and Astrophysics, University of
Chicago, Chicago, IL 60637, USA
6
Department of Physics, Grinnell College, Grinnell, IA 5011
2
USA
7
Space Telescope Science Institute, Baltimore, MD 21218, US
A
8
Department of Astronomy and Astrophysics, University of Ca
l-
ifornia, Santa Cruz, CA 95064 USA
9
Institut d’Astrophysique et de G ́eophysique, Universite ́
e de
Li ́ege, Li ́ege 1, Belgium
10
Centre de Recherche Astrophysique de Lyon, 69364 Lyon,
France
11
Institute for Astronomy, University of Hawaii at Manoa, Hon
-
olulu, HI, USA
by a thick hydrogen atmosphere (“mini-Neptunes”), or
water-dominated with a thick steam atmosphere (“wa-
ter 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 dis-
tant orbits (Hansen & Murray 2012; Chiang & Laughlin
2013; Raymond & Cossou 2014).
Detailed studies of super-Earth compositions offer im-
portant clues on their origins: presumably water-rich
planets must have formed beyond the ice line, while in-
situ formation models predict primarily rocky compo-
sitions with relatively water-poor hydrogen-dominated
atmospheres (e.g., Raymond et al. 2008). By combin-
ing radius measurements from
Kepler
with mass esti-
mates obtained using either the radial velocity or transit
timing techniques, it is possible to constrain the aver-
age densities and corresponding bulk compositions of the
super-Earths in the Kepler sample (e.g., Lithwick et al.
2012; Hadden & Lithwick 2013; 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 towards
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 mod-
2
els of planet formation and atmospheric mass loss for
close-in planets, this still leaves a wide range of plausi-
ble 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 com-
position. Planets with cloud-free, hydrogen dominated
atmospheres will have relatively large scale heights and
correspondingly strong absorption features during tran-
sit, 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 dif-
ficult 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 Can-
cri 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 sys-
tems 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 ́esert et al.
2011; Berta et al. 2012; Kreidberg et al. 2014). Al-
though 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 Cam-
era 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 con-
straints 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 us-
ing 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 den-
sity of 3
.
4
±
0
.
9 g cm
−
1
(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 temper-
ature between 700
−
1000 K depending on the effi-
ciency of energy transport to the planet’s night side.
If this planet has the same atmospheric composition
as GJ 1214b, its modestly higher atmospheric temper-
ature might prevent the formation of the cloud layer de-
tected 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 pres-
ence or absence of the water absorption band located
at 1.4
μ
m; this feature has now been robustly de-
tected in the atmospheres of several hot Jupiters (e.g.,
-0.15
-0.10
-0.05
-0.00
0.05
0.10
Time from Predicted Transit Center (d)
0.996
0.997
0.998
0.999
1.000
1.001
Relative Flux
Fig. 1.—
Raw white light photometry for the UT 2013 Dec 19
visit (top) and UT 2014 Jan 7 transit (bottom). Different scan
directions for the December visit are indicated as light (fo
rward
scan) and dark (reverse scan) blue filled circles. Scan direc
tions
for the January visit are plotted as yellow (forward scan) an
d red
(reverse scan) filled circles. The two light curves have been
offset
by a relative flux of 0.0015 for clarity. The first spacecraft o
rbit
has been trimmed from each visit, leaving four orbits per tra
nsit.
We have also normalized the light curves from each scan direc
tion
to one by dividing each by the median flux value; we do this in
order to remove a small offset in the measured fluxes from the tw
o
scan directions.
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
§
2 and the implications for the prop-
erties of this planet’s atmosphere in
§
3.
2.
OBSERVATIONS
We observed one transit of HD 97658b on UT 2013
Dec 19 and another on UT 2014 Jan 7 (GO 130501, PI
Knutson) using the G141 grism on the
HST
WFC3 in-
strument, 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 ap-
proximately seven hours per visit; our target was visi-
ble 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 uti-
lized 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 Fig. 1. We do not include these short ex-
3
0.19
0.20
0.21
0.22
0.23
0.24
0.25
-0.15
-0.10
-0.05
-0.00
0.05
0.10
Time from Predicted Transit Center (d)
0.22
0.23
0.24
0.25
0.26
0.27
Relative Sky Background (% )
Fig. 2.—
Estimated sky background as a percentage of the total
flux in the white-light aperture (i.e., summed over all wavel
engths)
for stacked images created using the method described in sec
tion
2.1. Sky backgrounds for the December visit are shown in the
top panel and background for the January visit are shown in th
e
bottom panel.
posures in our analysis, as they were taken in imaging
rather than spectroscopic mode.
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 1
.
4
′′
s
−
1
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
−
1
, 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 sum-
marize here. We begin with the raw sample up the ramp
images from the ima.fits files, which are bias and dark
1.2
1.4
1.6
Wavelength (microns)
2x10
6
4x10
6
6x10
6
8x10
6
Total Counts (electrons)
Fig. 3.—
A representative spectrum from the UT 2013 Dec 19
visit; this spectrum was created from the two dimensional im
age
by summing in the
y
(cross-dispersion) direction.
subtracted, and apply our own flat-fielding and wave-
length 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 cre-
ate 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 con-
tributions 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 wave-
length element (see Fig. 2), and subtract this estimated
sky background from each image.
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 ex-
tract 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 ex-
tended wings of the point spread function (see Fig. 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 trans-
mission 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
T DB
time standard following the
methods of Eastman et al. (2010).
Our next step is to create a spectral template by av-
eraging the ten spectra immediately before and after the
transit. We then fit this template spectrum to each indi-
vidual spectrum in our time series, allowing the relative
position and amplitude of the template spectrum to vary
4
-0.15
-0.10
-0.05
-0.00
0.05
0.10
Time from Predicted Transit Center (d)
0.9970
0.9975
0.9980
0.9985
0.9990
0.9995
1.0000
Relative Flux
Fig. 4.—
Normalized white light photometry with best-fit detec-
tor effects removed for the UT 2013 Dec 19 visit (top) and UT 201
4
Jan 7 transit (bottom). Different scan directions for the Dec
ember
visit are indicated as light (forward scan) and dark (revers
e scan)
blue filled circles. Scan directions for the January visit ar
e plotted
as yellow (forward scan) and red (reverse scan) filled circle
s. The
first spacecraft orbit has been trimmed from each visit, leav
ing four
orbits per transit, and the two light curves have been offset b
y a
relative flux of 0.0015 for clarity.
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 Fig. 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 direc-
tions, 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 inte-
gration time will be slightly longer than in the case of a
reverse scan. The forward and reverse spectra in our im-
ages 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 ben-
efit 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 in-
dividual 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.
-400
-200
0
200
400
-0.15
-0.10
-0.05
-0.00
0.05
0.10
Time from Predicted Transit Center (d)
-400
-200
0
200
Relative Flux (ppm)
Fig. 5.—
White light residuals after the best-fit detector and
transit light curves are removed for the UT 2013 Dec 19 visit (
top
panel) and UT 2014 Jan 7 transit (bottom panel). Different sca
n
directions for the December visit are indicated as light (fo
rward
scan) and dark (reverse scan) blue filled circles. Scan direc
tions
for the January visit are plotted as yellow (forward scan) an
d red
(reverse scan) filled circles.
2.1.1.
White Light Transit Fits
We take the two raw transit light curves plotted in
Fig. 1 and fit them each with a transit light curve, a lin-
ear function of time, and an exponential ramp in orbital
phase:
F
(
t
) =
c
1
(1 +
c
2
t
+
c
3
e
−
p/c
4
)
F
transit
(
t
)
(1)
where
c
1
−
c
4
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 usu-
ally 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 re-
set. 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 in-
strumental systematics is slightly different for the white
light curves calculated from the forward and reverse scan
directions (see Fig. 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 or-
bital 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 un-
certainty on our reported transit time for the first visit,
which spanned the transit ingress, was particularly sen-
sitive 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 ad-
dressed 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-
5
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 uncer-
tainties contributed by these two parameters. We al-
lowed 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 in-
clude 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 cal-
culate four-parameter nonlinear limb-darkening parame-
ters (Claret 2000) for our transit light curves, with the
relative intensity at each position given as:
I
(
μ
)
I
(1)
= 1
−
4
∑
k
=1
a
k
(1
−
μ
k/
2
)
(2)
where
I
(1) is the specific intensity at the center of the
stellar disk,
a
k
is the
k
th limb-darkening coefficient, and
μ
=
cos
(
θ
) where
θ
is the angle between the 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 temper-
ature 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 con-
straints on the stellar density from the transit light curve.
We also tried models with effective temperatures rang-
ing 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 inclu-
sion 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 band-
pass. 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
6
steps. In our initial fits
with fixed
i
and
a/R
⋆
, we set the uncertainties on indi-
TABLE 1
Transit Parameters from Joint
HST
and
Spitzer
Fits
Parameter
Value
Global Values
i
(
◦
)
89
.
85
±
0
.
48
a
a/R
⋆
26
.
24
±
1
.
21
P
(days)
b
9
.
489264
±
0
.
000064
T
0
(BJD
T DB
)
b
2456665
.
46415
±
0
.
00078
UT 2013 Aug 10 Spitzer Transit
c
R
p
/R
⋆
0
.
02778
±
0
.
00073
T
c
(BJD
T DB
)
2456523
.
12537
±
0
.
00049
UT 2013 Dec 19 HST Transit
R
p
/R
⋆
0
.
03012
±
0
.
00087
T
c
(BJD
T DB
)
2456646
.
48556
±
0
.
00071
UT 2014 Jan 7 HST Transit
R
p
/R
⋆
0
.
03111
±
0
.
00080
T
c
(BJD
T DB
)
2456665
.
4584
±
0
.
0056
a
We limited the value of the inclination to be less than or equa
l
to 90
◦
in our fits. We find a 1
σ
lower limit on the inclination of
89.65
◦
, and a 2
σ
limit of 88.79
◦
.
b
We calculate our updated ephemeris using the published tran
sit
times from Dragomir et al. (2013) and Van Grootel et al. (2014
)
and our two new transit measurements.
T
0
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.
c
We find that for the
Spitzer
data the uncertainties from the co-
variance matrix are larger than those from the prayer bead an
alysis,
and report the covariance errors here.
vidual 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 tech-
nique results in uncertainties that are 2
−
3 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 co-
variance methods implicitly assume random Gaussian-
distributed noise (Carter & Winn 2009). We plot the
2D probability distributions from the MCMC analysis
and confirm that there are no significant correlations be-
tween any of the fit parameters in this case. For the
case where we fit the
Spitzer
and
HST
transits simul-
taneously while allowing
i
and
a/R
⋆
to vary as free pa-
rameters, 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 pre-
ferred 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 normal-
ized transit light curves after dividing out our best fit for
6
300
400
500
600
700
Transit Center [BJD-2456000]
-20
-15
-10
-5
0
5
O-C (min)
Fig. 6.—
Observed minus calculated transit times using our up-
dated ephemeris. Previously published
MOST
(Dragomir et al.
2013) and
Spitzer
transit times (Van Grootel et al. 2014) are
shown as open circles, with our new
HST
observations are shown
as filled circles. The second
HST
transit has a larger uncertainty
on the transit time because it does not include any points dur
ing
ingress or egress, which provide the strongest constraints
on the
measured center of transit time.
detector effects in Fig. 4. The white light curve resid-
uals from UT 2013 Dec 19 and UT 2014 Jan 7 have a
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 resid-
uals from these fits in Fig. 5, and the offset between the
observed and calculated transit times from our updated
orbital ephemeris in Fig. 6.
2.1.2.
Differential Transit Fits
We estimate the planet’s wavelength-dependent tran-
sit 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 differ-
ential 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 calcu-
lated 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 al-
lowed 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 di-
rections. 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
TABLE 2
Binned Four Parameter Nonlinear
Limb-Darkening Coefficients
a
Wavelength
b
c
1
c
2
c
3
c
4
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.434
c
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.491
c
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.5
d
0.8713 -1.4634 1.4955 -0.5581
a
The spectra extracted from the UT 2013 Dec 19 and UT
2014 Jan 7 visits have slightly different wavelength solu-
tions, and we therefore created custom limb darkening co-
efficients for each visit. A full table of both the binned
and unbinned coefficients for both visits is available upon
request.
b
These 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 ap-
propriately for this resolution (4x higher than shown here)
.
c
There 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 conti
n-
uous across all wavelengths considered here.
d
Limb-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).
effective temperature by
±
50 K changes the resulting
transit depths by 1 ppm, which we include in our analy-
sis as a separate systematic error term.
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 permuta-
tion 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
7
1.1
1.2
1.3
1.4
1.5
1.6
Wavelength (micron)
-100
-50
0
50
100
150
Differential Transit Depth (ppm)
Fig. 7.—
Wavelength-dependent relative transit depths measured
for the UT 2013 Dec 19 (blue filled circles) and UT 2014 Jan 7
(red filled circles) visits. The white-light transit depths
have been
subtracted from each visit.
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 chan-
nel as our final error estimate, which allows us to account
for time-correlated noise in our data while avoiding pos-
sible under-estimation of uncertainties due to the limited
number of measurements available for the residual per-
mutation analysis.
Next we bin the resulting wavelength-dependent ra-
dius values by a factor of four to create the transmis-
sion spectrum shown in Fig. 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 chan-
nels 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 Gaus-
sian 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 av-
erage 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.
We take the error-weighted average of the two indi-
vidual transmission spectra to create the combined spec-
trum shown in Fig. 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
TABLE 3
Wavelength Dependent Transit Depths from
Template Fitting Method
Wavelength Depth White Noise Error
a
Total Error
a
μ
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
a
White noise measurement errors are estimated using the cova
ri-
ance matrix, which implicitly assumes that individual meas
ure-
ment errors are Gaussian distributed and uncorrelated. The
total
measurement errors are calculated by comparing the covaria
nce
errors to a residual permutation method that better account
s for
time-correlated noise and taking the larger of the two in eac
h
wavelength channel. Uncertainty in the limb-darkening mod
els
from the stellar effective temperature contributes an addit
ional
error of 1 ppm in each bin, which we do not include here.
the two error estimates as previously described. As an
additional consistency check, we also carried out a ver-
sion of our analysis in which we allowed the light curves
derived from the forward and reverse scans to have dif-
ferent planet-star radius ratios. We found that we ob-
tained consistent radius estimates from both scan direc-
tions, indicating that the different behaviors visible in
the white-light curves are effectively removed in the dif-
ferential time series by our template fitting technique.
2.2.
White Light Residual Fitting Method
We also obtained an independent estimate of
the transmission spectrum following the approach of
Kreidberg et al. (2014) and Stevenson et al. (2014). Be-
ginning with the raw images, this analysis used an en-
tirely different set of codes than those described in the
previous sections. In this version of the analysis we
treated each up-the-ramp sample in the ima.fits images
as an independent subexposure. For each subexposure,
we interpolated all rows to a common wavelength scale to
account for the changing dispersion solution with spatial
position on the detector. We estimated the background
by making conservative masks around the stellar spec-
tra and taking the median of the unmasked pixels. We
subtracted the background and optimally extracted the
spectra. To combine the data from individual subexpo-
8
TABLE 4
Wavelength Dependent Transit Depths from White Light
Residual Method
Wavelength Depth White Error
a
Total Error
a
Reduced
χ
2
μ
m
ppm
ppm
1.145
909
20
23
0.81
1.163
940
19
24
0.77
1.182
896
19
23
0.93
1.200
928
19
22
0.77
1.218
910
19
26
0.82
1.237
895
18
24
0.82
1.255
922
18
20
0.76
1.274
936
19
23
0.88
1.292
944
18
30
0.83
1.311
970
18
29
1.04
1.329
933
18
21
0.87
1.348
922
18
26
0.80
1.366
929
18
26
0.81
1.384
922
18
28
0.83
1.403
915
18
18
0.89
1.421
983
18
18
0.84
1.440
978
18
19
0.85
1.458
960
18
21
0.75
1.477
936
19
26
0.68
1.495
924
19
26
0.89
1.513
962
19
25
0.91
1.532
941
19
29
0.89
1.550
942
19
32
0.97
1.569
960
20
42
0.95
1.587
970
20
26
0.87
1.606
1001
20
29
0.87
a
White noise measurement errors are estimated using a Markov
Chain
Monte Carlo analysis, which implicitly assumes that indivi
dual measure-
ment errors are Gaussian distributed and uncorrelated. Tot
al measure-
ment errors are calculated by comparing the MCMC errors to a r
esidual
permutation method that better accounts for time-correlat
ed noise and
taking the larger of the two in each wavelength bin.
1.1
1.2
1.3
1.4
1.5
1.6
Wavelength (micron)
850
900
950
1000
Normalized Transit Depth (ppm)
Fig. 8.—
Wavelength-dependent transit depths averaged over
the two visits, where the depths are defined as the square of th
e
planet-star radius ratio
R
p
/R
⋆
in each band. Depths derived us-
ing the spectral template fitting technique (Deming et al. 20
13;
Knutson et al. 2014) are shown as filled circles, and depths fr
om
the white light residual fitting technique (Kreidberg et al.
2014)
are shown as open circles. No offset has been applied to either
data set, demonstrating that the average transit depths are
also in
good agreement.
sures, we summed the spectra by column. The final step
in the reduction process is to correct for drift of the spec-
tra in the dispersion direction over the course of a visit.
We used the first exposure from the first visit as a tem-
plate and shifted all subsequent spectra to the template
wavelength scale. The spectra shifted by a total of 0.3
pixels over the five orbits contained in our observations,
1.1
1.2
1.3
1.4
1.5
1.6
Wavelength (micron)
15
20
25
30
35
Transit Depth Uncertainty (ppm)
Fig. 9.—
Comparison of uncertainties on the reported tran-
sit depths under the assumption of either white, Gaussian no
ise
(open circles) or allowing for time-correlated noise using
the resid-
ual permutation method for estimating uncertainties (fille
d cir-
cles). Errors derived using the spectral template fitting te
chnique
(Deming et al. 2013; Knutson et al. 2014) are shown as blue cir
-
cles, and depths from the white light residual fitting techni
que
(Kreidberg et al. 2014) are shown as red circles.
which is larger than the approximately 0.01 pixel drift
observed in previous scanning mode observations of GJ
1214b (Kreidberg et al. 2014). This increased drift may
be related to the longer scan length and faster scan rate
utilized for these observations as compared to GJ 1214b.
We binned the spectra in four-pixel-wide channels,
yielding 26 spectrophotometric light curves between 1.15
and 1.61
μ
m. The light curves show orbit-long ramp-
like systematics that are characteristic of WFC3 data.
We correct for these systematics using the divide-white
technique, which assumes that the observed effects have
the same shape across all wavelengths. We fit each spec-
troscopic light curve with a transit model multiplied by a
scaled vector of systematics from the best-fit white light
curve and a linear function of time. We fix
i
and
a/R
⋆
to
the values reported in Table 1. The fit to each channel
has six free parameters: one transit depth, four scaling
factors for the systematics (one for each visit and scan
direction), and the linear slope. As before, we calculate
the four-parameter nonlinear limb-darkening coefficients
using a
PHOENIX
stellar atmosphere model where we take
the flux-weighted average of the theoretical stellar inten-
sity profile within each photometric bandpass. We set
the uncertainties on individual points equal to the sum
of the the photon noise and read noise in quadrature. We
list the best-fit values and corresponding errors in Table
4.
We report uncertainties on the transit depths corre-
sponding to 1
σ
confidence intervals from either a Markov
chain Monte Carlo (MCMC) fit, which implicitly assumes
white (Gaussian and uncorrelated) noise, or a residual
permutation analysis that better accounts for any time-
correlated noise present in the data. We take the larger
of the two errors in each wavelength bin as our final un-
certainties and provide the MCMC only errors separately
in Table 4 for comparison. Our residual permutation er-
rors are on average 30% larger than those obtained with
MCMC. This suggests that there is some time-correlated
noise in the light curves, which is most likely the result of
imperfect corrections for visit- and orbit-long systematic
trends in the data.
9
3.
DISCUSSION AND CONCLUSIONS
1.1
1.2
1.3
1.4
1.5
1.6
Wavelength (micron)
850
900
950
1000
Normalized Transit Depth (ppm)
Fig. 10.—
Wavelength-dependent transit depths averaged over
the two visits (black filled circles). Four different atmosph
ere mod-
els are shown for comparison: a solar metallicity model (red
solid
line), a 50
×
solar metallicity model (red dashed line), a pure wa-
ter model (solid blue line), and a solar metallicity model wi
th an
opaque cloud deck at 1 mbar (solid green line). The average de
pth
for each model has been normalized to match the average measu
red
transit depth in this plot.
We used the the WFC3 instrument on the
Hubble
Space Telescope
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 dom-
inated by systematic noise, as illustrated in Fig. 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 pe-
riod with an uncertainty that is a factor of thirty smaller
than in Van Grootel et al. (2014). Our average radius ra-
tio 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 ac-
tivity as this star has a log(
R
′
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 ob-
servation 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 Fig. 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 Fig. 9), making these observations
the most precise measurement of a planetary transmis-
sion 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 up-
ward slope visible towards redder wavelengths. These
features could be instrumental in nature, as the columns
that make up the wavelength bins closest to the discon-
tinuity are located on a boundary between quadrants on
the WFC3 array. The statistical significance of these fea-
tures is marginal, and we therefore do not consider them
in our subsequent comparison to atmosphere models for
this planet.
In Fig. 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 pro-
files of these models are calculated assuming full redis-
tribution of energy to the planet’s night side and an
albedo which varies according to composition. We con-
sider both cloud-free models with solar and 50
×
solar
metallicities, as well as a pure water model and a so-
lar 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):
Signif icance
=
χ
2
−
ν
√
2
ν
(3)
where
χ
2
is calculated by comparing our averaged trans-
mission 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 met-
ric 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 metal-
licity 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 so-
lar metallicity model with clouds at 1 mbar is also con-
sistent 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 atmo-
sphere 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 corre-
spondingly small scale height, or the presence of a high
cloud deck that obscures the expected water absorption
feature in a hydrogen-dominated atmosphere. We con-
strain 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 hy-
drogen and water, and find that in this scenario the at-
mosphere has to be at least 20% water by number in
order to be consistent with our data at the 3
σ
level. We
list the
χ
2
values and level of disagreement for all models
in Table 5.
The conclusion that HD 97658b’s transmission spec-
trum 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 plan-
ets, a more precise measurement of HD 97658b’s trans-
mission spectrum will eventually allow us to distinguish
between high clouds and a cloud-free, metal-rich atmo-