ORE Open Research Exeter
TITLE
A ground-based optical transmission spectrum of WASP-6b
AUTHORS
Jordán, A; Espinoza, N; Rabus, M; et al.
JOURNAL
Astrophysical Journal
DEPOSITED IN ORE
16 June 2016
This version available at
http://hdl.handle.net/10871/22131
COPYRIGHT AND REUSE
Open Research Exeter makes this work available in accordance with publisher policies.
A NOTE ON VERSIONS
The version presented here may differ from the published version. If citing, you are advised to consult the published version for pagination, volume/issue and date of
publication
The Astrophysical Journal, 778:184 (13pp), 2013 December 1 doi:10.1088/0004-637X/778/2/184
C
2013. The American Astronomical Society. All rights reserved. Printed in the U.S.A.
A GROUND-BASED OPTICAL TRANSMISSION SPECTRUM OF WASP-6b
Andr
´
es Jord
´
an
1
,N
´
estor Espinoza
1
, Markus Rabus
1
, Susana Eyheramendy
2
, David K. Sing
3
,
Jean-Michel D
´
esert
4,5
,G
´
asp
´
ar
´
A. Bakos
6,11,12
, Jonathan J. Fortney
7
, Mercedes L
´
opez-Morales
8
,
Pierre F. L. Maxted
9
, Amaury H. M. J. Triaud
10,13
, and Andrew Szentgyorgyi
8
1
Instituto de Astrof
´
ısica, Facultad de F
´
ısica, Pontificia Universidad Cat
´
olica de Chile, Av. Vicu
˜
na Mackenna 4860, 7820436 Macul, Santiago, Chile
2
Departmento de Estad
´
ıstica, Facultad de Matem
´
aticas, Pontificia Universidad Cat
´
olica de Chile, Av. Vicu
˜
na Mackenna 4860, 7820436 Macul, Santiago, Chile
3
School of Physics, University of Exeter, Stocker Road, Exeter EX4 4QL, UK
4
CASA, Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA
5
Division of Geological and Planetary Sciences, California Institute of Technology, MC 170-25 1200, East California Boulevard, Pasadena, CA 91125, USA
6
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA
7
Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA
8
Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
9
Astrophysics Group, Keele University, Staffordshire ST5 5BG, UK
10
Department of Physics, and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
Received 2013 September 29; accepted 2013 October 21; published 2013 November 13
ABSTRACT
We present a ground-based optical transmission spectrum of the inflated sub-Jupiter-mass planet WASP-6b. The
spectrum was measured in 20 spectral channels from 480 nm to 860 nm using a series of 91 spectra over a
complete transit event. The observations were carried out using multi-object differential spectrophotometry with
the Inamori-Magellan Areal Camera and Spectrograph on the Baade Telescope at Las Campanas Observatory. We
model systematic effects on the observed light curves using principal component analysis on the comparison stars
and allow for the presence of short and long memory correlation structure in our Monte Carlo Markov Chain analysis
of the transit light curves for WASP-6. The measured transmission spectrum presents a general trend of decreasing
apparent planetary size with wavelength and lacks evidence for broad spectral features of Na and K predicted by
clear atmosphere models. The spectrum is consistent with that expected for scattering that is more efficient in the
blue, as could be caused by hazes or condensates in the atmosphere of WASP-6b. WASP-6b therefore appears to
be yet another massive exoplanet with evidence for a mostly featureless transmission spectrum, underscoring the
importance that hazes and condensates can have in determining the transmission spectra of exoplanets.
Key words: planetary systems – planets and satellites: atmospheres – techniques: spectroscopic
Online-only material: color figures
1. INTRODUCTION
Due to their fortuitous geometry, transiting exoplanets allow
the determination of physical properties that are inaccessible or
hard to reach for non-transiting systems. One of the most excit-
ing possibilities enabled by the transiting geometry is to measure
atmospheric properties of exoplanets without the need to resolve
them from their parent star through the technique of transmis-
sion spectroscopy. In this technique, the atmospheric opacity at
the planet terminator is probed by measuring the planetary size
via transit light-curve observations at different wavelengths. The
measurable quantity is the planet-to-star radius ratio as a func-
tion of wavelength, (R
p
/R
∗
)(λ) ≡ k(λ), and is termed the trans-
mission spectrum. The measurement of a transmission spectrum
is a challenging one, with one atmospheric scale height H trans-
lating to a signal of order 2Hk ≈ 10
−4
for hot Jupiters (e.g.,
Brown 2001). The requirements on precision favor exoplanets
with large atmospheric scale heights, large values of k (e.g.,
systems transiting M dwarfs), and orbiting bright targets due to
the necessity of acquiring a large number of photons to reach
the needed precision.
The first successful measurement by transmission spec-
troscopy was the detection with the Hubble Space Telescope
11
Alfred P. Sloan Fellow.
12
Packard Fellow.
13
Fellow of the Swiss National Science Foundation.
(HST) of absorption by Na i in the hot Jupiter HD 209458b
(Charbonneau et al. 2002). The signature of Na was 2–3 times
weaker than expected from clear atmosphere models, providing
the first indications that condensates can play an important role
in determining the opacity of their atmospheres as seen in trans-
mission (e.g., Fortney 2005, and references therein). Subsequent
space-based studies have concentrated largely on the planets or-
biting the stars HD 209458 and HD 189733 due to the fact that
they are very bright stars and therefore allow the collection of
a large number of photons even with the modest aperture of
space-based telescopes. A recent study of all the transmission
spectra available for HD 189733, spanning the range from 0.32
to 24 μm, points to a spectrum dominated by Rayleigh scat-
tering over the visible and near-infrared range, with the only
detected feature being a narrow resonance line of Na (Pont
et al. 2013). For HD 209458, Deming et al. (2013) present new
WFC3 data combined with previous Space Telescope Imaging
Spectrograph data (Sing et al. 2008), resulting in a transmission
spectrum spanning the wavelength range 0.3–1.6 μm. They con-
clude that the broad features of the spectrum are dominated by
haze and/or dust opacity. In both cases the spectra are different
from those predicted by clear atmosphere models that do not
incorporate condensates.
In order to further our understanding of gas giant atmo-
spheres, it is necessary to build a larger sample of systems with
measured transmission spectra. Hundreds of transiting exoplan-
ets, mostly hot gas giants, havebeendiscovered by ground-based
1
The Astrophysical Journal, 778:184 (13pp), 2013 December 1 Jord
´
an et al.
Tab le 1
List of Comparison Stars
2MASS Identifier
2MASS-23124095-2243232
2MASS-23124836-2252099
2MASS-23124448-2253190
2MASS-23124428-2256403
2MASS-23114068-2248130
2MASS-23113937-2250334
2MASS-23114820-2256592
surveys such as HATNet (Bakos et al. 2004), WASP (Pollacco
et al. 2006), KELT (Pepper et al. 2007), XO (McCullough et al.
2005), TRES (Alonso et al. 2004), and HATSouth (Bakos et al.
2013), with magnitudes within reach of the larger collecting ar-
eas afforded by ground-based telescopes but often too faint for
HST.
14
The ground-based observations have to contend with the
atmosphere and instruments lacking the space-based stability of
HST, but despite these extra hurdles the pace of ground-based
transmission spectra studies is steadily increasing. Following
the ground-based detection of Na i in HD 189733b (Redfield
et al. 2008) and confirmation of Na i in HD 209358b (Snellen
et al. 2008), Na i has been additionally reported from the ground
in WASP-17b (Wood et al. 2011; Zhou & Bayliss 2012) and
XO-2b (Sing et al. 2012). K i has been detected in XO-2b (Sing
et al. 2011a) and the highly eccentric exoplanet HD 80606b
(Col
´
on et al. 2012). All of these studies have used high-
resolution spectroscopy or narrowband photometry to specif-
ically target resonant lines of alkali elements. Recently, a detec-
tion of Hα has been reported from the ground for HD 189733b
(Jensen et al. 2012), complementing previous space-based de-
tection of Lyα and atomic lines in the UV with HST for
HD 189733b and HD 209358b (Vidal-Madjar et al. 2003, 2004;
Lecavelier Des Etangs et al. 2010).
Differential spectrophotometry using multi-object spectro-
graphs offers an attractive means to obtain transmission spectra
given the possibility of using comparison stars to account for
the various systematic effects that affect the spectral time series
obtained. Using such spectrographs, transmission spectra in the
optical have been obtained for GJ 1214b (Bean et al. 2011;
610–1000 nm with VLT/FORS) and recently for WASP29-b
(Gibson et al. 2013; 515–720 nm with Gemini/GMOS), with
both studies finding featureless spectra. In the near-infrared
Bean et al. (2013) present a transmission spectrum in the range
1.25–2.35 μm for WASP-19b, using MMIRS on Magellan. In
this work we present an optical transmission spectrum of an-
other planet, WASP-6b, an inflated sub-Jupiter-mass (0.504 M
J
)
planet orbiting a V = 11.9 G dwarf (Gillon et al. 2009), in the
in the range 471–863 nm.
2. OBSERVATIONS
The transmission spectrum of WASP-6b was obtained per-
forming multi-object differential spectrophotometry with the
Inamori-Magellan Areal Camera and Spectrograph (IMACS;
Dressler et al. 2011) mounted on the 6.5 m Baade telescope at
Las Campanas Observatory. A series of 91 spectra of WASP-6
and a set of comparison stars were obtained during a tran-
sit of the hot Jupiter WASP-6b in 2010 October 3 with the
14
The Kepler mission (Borucki et al. 2010) has discovered thousands of
transiting exoplanet candidates, but the magnitudes of the hosts are usually
significantly fainter than the systems discovered by ground-based surveys,
making detection of their atmospheres more challenging.
3
3.5
4
4.5
5
5.5
5000 6000 7000 8000
log
10
(Counts)
Wavelength
(
˚
A
)
Comparison spectra
WASP-6 spectrum
Figure 1. Extracted spectra for WASP-6 and the seven comparison stars used
in this work for a typical exposure.
f/2 camera of IMACS, which provides an unvignetted circu-
lar field of view of radius r ≈ 12 arcmin. The large field
of view makes IMACS a very attractive instrument for multi-
object differential spectrophotometry as it allows us to search
for suitable comparison stars that have as much as possi-
ble similar magnitude and colors as the target star. The me-
dian cadence of our observations was 224 s, and the expo-
sure time was set to 140 s, except for the first eight exposures
when we were tuning the exposure level and whose exposure
times were {30, 120, 150, 150, 150, 130, 130, 130} s. The count
level of the brightest pixel in the spectrum of WASP-6 was
≈43,000 ADU, i.e., ≈65% of the saturation level. In addition
to WASP-6, we observed 10 comparison stars of comparable
magnitude, seven of which had the whole wavelength range
of interest (≈4700–8600 Å) recorded in the CCD with enough
signal-to-noise ratio. The seven comparison stars we used are
listed in Table 1. The integrated counts over the wavelength
range of interest for the spectrum of WASP-6 were typically
≈3.6 × 10
8
electrons, giving a Poisson noise limit for the
white-light light curve of ≈0.06 mmag. Each star was observed
through a 10 × 10 arcsec
2
slit in order to avoid the adverse
effects of variable slit losses. We used the 300-l+17 grating as
dispersing element, which gave us a seeing-dependent resolu-
tion Δλ that was ≈5 Å under 0.7 arcsec seeing and a dispersion
of 1.34 Å pixel
−1
. In addition to the science mask, we obtained
HeNeAr arc lamps through a mask that had slits at the same
position as the science mask but with slit widths in the spectral
direction of 0.7 arcsec. Observing such masks is necessary in
order to produce well-defined lines that are then used to define
the wavelength solution.
The extracted spectra of WASP-6 and the seven comparison
stars we used are shown for a typical exposure in Figure 1.The
conditions throughout the night were variable. The raw light
curves constructed with the integrated counts over the whole
spectral range for WASP-6 and the comparison stars are shown
as a function of time in Figure 2. Besides the variation due to
varying airmass (and the transit for WASP-6), there were periods
with strongly varying levels of transparency concentrated in
the period of time 0–2 hr after mid-transit. The seeing was in
the range ≈0.
6–0.
8. In order to maintain good sampling of the
point-spread function in the spatial direction, we defocused the
telescope slightly in the periods of best seeing. Changes in
seeing and transparency left no noticeable traces in the final
light curves.
2
The Astrophysical Journal, 778:184 (13pp), 2013 December 1 Jord
´
an et al.
−0.25
0
0.25
0.5
0.75
−2 −10123
Magnitude - mean magnitude + constant
Time from mid-transit
(
hours
)
Raw comparison lightcurves
Raw WASP-6 lightcurve
Figure 2. Raw light curves for WASP-6 and seven comparison stars used in this
work as a function of time.
3. DATA REDUCTION
3.1. Background and Sky Subtraction
After subtracting the median value of the overscan region to
every image, an initial trace of each spectrum was obtained by
calculating the centroid of each row, which are perpendicular to
the dispersion direction. Each row was then divided into three
regions: a central region, which contains the bulk of the light
of the star; a middle, on-slit region, which is dominated by sky
continuum and line emission; and an outer, out-of-slit region,
which contains a smooth background outside the slit arising
from, e.g., scattered light. The middle and outer regions have
components on each side of the spectrum. The outermost region
was used to determine a smooth background that varies slowly
along the dispersion direction. The median level was obtained in
the outer regions on either side of the slit, and then a third-order
polynomial was used to estimate the average background level
as a function of pixel in the dispersion direction. This smooth
background component was then subtracted from the central
and middle regions. Then a Moffat function plus a constant
level c
i
was fit robustly to each background-subtracted row in
those regions. The estimated c
i
(one per row) was then subtracted
from the central and middle regions in order to obtain a spectrum
where only the stellar contribution remains. It is necessary to
estimate the sky emission on a row-per-row basis as skyemission
lines have a wide, box-shaped form with sharp boundaries due
to the fact that they fully illuminate the wide slit.
3.2. Fine Tracing and Spectrum Extraction
The background- and sky-subtracted spectrum was traced by
an algorithm that cross-correlates each slice perpendicular to the
wavelength direction with a Gaussian in order to find the spectral
trace. The centers of the trace were then fitted robustly with a
fourth-order polynomial. This new tracing procedure served as
a double check for the centers obtained via the centroid method
in the background and sky subtraction part of the data reduction
process; both methods gave traces consistent with each other.
With the trace in hand, the spectrum was extracted by using
a simple extraction procedure, i.e., summing the flux on each
row ±15 pixels from the trace position at that row. We also
tried optimal extraction (Marsh 1989), but it led to additional
systematic effects when analyzing the light curves,
15
and in
any case optimal extraction is not expected to give significant
gains over simple extraction at the high signal-to-noise levels
we are working with here. We also took spectroscopic flats at
the beginning of the night with a quartz lamp and reduced the
data using both flat-fielded and non-flat-fielded spectra. The
results where consistent when using both alternate reductions,
but the flat-fielded spectra showed higher dispersion in the final
transmission spectrum. We therefore used the non-flat-fielded
spectra in the present work.
3.3. Wavelength Calibration
The extracted spectra were calibrated using NeHeAr lamps
taken at the start of the night. The wavelength solution was
obtained by the following iterative procedure: pixel centers
of lines with known wavelengths were obtained by fitting
Gaussians to them, and then all the pixel centers, along with
the known wavelengths of the lines, were fitted by a sixth-order
Chebyshev polynomial. We checked the absolute deviation of
each line from the fit and removed the most deviant one from our
sample, repeating the fit without it. This process was iterated,
removing one line at a time, until an rms of less than 2000 m s
−1
was obtained. The rms of the final wavelength solution was
≈1200 m s
−1
, using 27 lines.
The procedure explained in the preceding paragraph served
to wavelength-calibrate the first spectrum of the night closest
in time to the NeHeAr lamps. In order to measure and correct
for wavelength shifts throughout the night, the first spectrum
was cross-correlated with the subsequent ones in pixel-space in
order to find the shifts in wavelength-space. If λ
t
0
,s
(p)isthe
wavelength solution at time t
0
(the beginning of the night) for
star s as a function of the pixel p, then the wavelength solution at
time t is just λ
t,s
(p +δp
t,s
), where δp
t,s
is the shift in pixel-space
found by cross-correlating the spectrum of star s taken at time
t
0
with the one taken at time t. Finally, each spectrum was fitted
with a b-spline in order to interpolate each of the spectra into a
common wavelength grid with pixel size 0.75 Å.
4. MODELING FRAMEWORK
The observed signal of WASP-6 is perturbed with respect to
its intrinsic shape, which we assume ideally to be a constant
flux, F. This constant flux is multiplied by the transit signal,
f (t; θ ), which we describe parametrically using the formalism
of Mandel & Agol (2002). In what follows θ represents the
vector of transit parameters. The largest departure from this
idealized model in our observations will be given by systematic
effects arising from atmospheric and instrumental effects, which
are assumed to act multiplicativelyon our signals. We will model
the logarithm of the observed flux, L(t), as
L(t; θ ) = S(t)+log
10
f (t; θ )+log
10
F + (t), (1)
where S(t) represents the (multiplicative) perturbation to the
star’s flux, which we will refer to in what follows as the
perturbation signal, and (t) is a stochastic signal that represents
the noise in our measurements (under the term noise we will also
include potential variations of the star that are not accounted for
in the estimate of the deterministic S(t) and that can be modeled
by a stochastic signal).
15
Optimal extraction assumes that the profile along the wavelength direction
is smooth enough to be approximated by a low-order polynomial. However,
this assumption is not always valid. In particular, we found that fringing in the
reddest part of the spectra induces fluctuations in the extracted flux with
wavelength due to the inadequacy of the smoothness assumption.
3
The Astrophysical Journal, 778:184 (13pp), 2013 December 1 Jord
´
an et al.
4.1. Modeling the Perturbation Signal
4.1.1. Estimation of Systematic Effects via Principal
Component Analysis of the Comparison Stars
Each star in the field is affected by a different perturbation
signal. However, these perturbation signals have in common
that they arise from the same physical and instrumental sources.
In terms of information, this is something we want to take ad-
vantage of. We model this by assuming that a given pertur-
bation signal is in fact a linear combination of a set of signals
s
i
(t), which represent the different instrumental and atmospheric
effects affecting all of our light curves, i.e.,
S
k
(t) =
K
i=1
α
k,i
s
i
(t). (2)
Note that this model for the perturbation signal so far includes
the popular linear and polynomial trends (e.g., s
i
(t) = t
i
).
According to this model, the logarithm of the flux of each
of N stars without a transiting planet in our field can be
modeled as
L
k
(t; α) = S
k
(t; α)+log
10
F
k
+
k
(t), (3)
where α denotes the set of parameters {α
k,i
}
K
i=1
. In the case in
which we have a set of comparison stars, we can see each of them
as an independent (noisy) measurement of a linear combination
of the signals s
i
(t) in Equation (2). A way of obtaining those
signals is by assuming that the s
i
(t) are uncorrelated random
variables, in which case these signals are easily estimated
by performing a principal component analysis (PCA) of the
mean-subtracted light curves of the comparison stars. Given N
comparison stars, one can estimate at most N components, and
thus we must have K N. As written in Equation (3), we
cannot separate s
i
(t) from
k
(t), and in general the principal
components will have contributions from both terms. If s
i
(t)
k
(t), the K principal components that contribute most to the
signal variance will be dominated by the perturbation signals,
but some projection of the
k
into the estimates of s
i
is to
be expected.
4.1.2. Selecting the Number of Principal Components
In our case, the number of components K is unknown a priori.
We need therefore to determine an optimal number of principal
components to describe the perturbation signal, taking into
consideration that there is noise present in the light curves of the
comparison stars and, thus, some of the principal components
obtained are mostly noise. There are several possibilities for
doing this depending on what we define as optimal. We will
determine the optimal number of components as the minimum
number of components that are able to achievethebest predictive
power allowed by the maximum set of N components available.
As a measure of predictive power we use a k-fold cross-
validation procedure (Hastie et al. 2007). k-fold cross-validation
is a procedure that estimates prediction error, i.e., how well
a model predicts out-of-sample data. The idea is to split the
datapoints into k disjoint groups (called folds). A “validation”
fold is left out, and a fit is done with the remaining “training”
folds, allowing us to predict the data in the validation fold that
was not used in this fitting procedure. This procedure is repeated
for all folds. Denoting the datapoints by y
i
and the values
predicted on the kth fold by the cross-validation procedure by
f
−k
i
, an estimate for the prediction error is
ˆ
CV =
1
N
N
i=1
L(y
k
i
− f
−k
i
),
where
L(·) is the loss function. Examples of loss functions are
the
L
1
norm (L
1
(x) =|x|)ortheL
2
norm (L
2
(x) = x
2
).
In our case, the light curves of the N comparison stars are
used to estimate l<Nprincipal components. These l principal
components, which are a set of light curves {s
i
}
l
i=1
, are our
estimates of the systematic effects, and we use the out-of-transit
part of the light curve of WASP-6 as the validation data by fitting
it with the {s
i
}
l
i=1
. In more precise terms, if y(t
k
) denote the time
series of the out-of-transit portion of the light curve of WASP-6,
we apply k-fold cross-validation by considering a model of the
form y(t
k
) =
l
i=1
α
i
s
i
(t
k
).
4.2. Joint Parameter Estimation for Transit
and Stochastic Components
In the past sub-sections we set up an estimation process for
the signal given in Equation (2) using PCA. It remains to specify
a model for the stochastic signal that we have termed noise, i.e.,
the (t) term in Equation (1). As noted above, the principal
components will absorb part of the (t), and so our estimate of
the noise may not necessarily accurately reflect the (t )termin
Equation (1) assuming that the model holds. Nonetheless, this
is of no consequence as we just aim to model the residuals after
the time series has been modeled with the {s
i
}
l
i=1
. While we
still call this term (t) in what follows, one should bear in mind
this subtlety. An important feature of the correlated stochastic
models we consider is that they can model trends. The {s
i
}
l
i=1
are obtained from the comparison stars, and while the hope is
that they capture all of the systematic effects, it is possible that
some systematic effectsunique to the target star are not captured.
The stochastic “noise” models considered below that have time
correlations can in principle capture remnant individual trends
particular to WASP-6.
We make use of Markov Chain Monte Carlo (MCMC;
see, e.g., Ford 2005) algorithms to obtain estimates of the
posterior probability distributions of our parameters, θ,α,η,
given a data set y, where we have introduced a new set of
parameters η characterizing a stochastic component (see below).
The posterior distribution p(θ,α,η|y) is obtained using a prior
distribution for our parameters p(θ,α,η) and a likelihood
function, p(y|θ,α,η). Following previous works (e.g., Carter &
Winn 2009; Gibson et al. 2012), we assume that the likelihood
function is a multivariate Gaussian distribution given by
p(y|θ,α,η) =
1
(2π)
n/2
|
η
|
1/2
exp
−
1
2
(y − g(θ,α))
T
× Σ
η
−1
(y − g(θ,α))
, (4)
where g(θ,α) is the function that predicts the observed data-
points and
η
is the covariance matrix that depends on the set
of parameters η. It is the structure of this matrix that defines the
type of noise of the residuals. Previous works have proposed
to account for time-correlated structure in the residuals using
flicker-noise models, where it is assumed that the noise follows
a power spectral density (PSD) of the form 1/f (Carter & Winn
2009), and Gaussian processes, where the covariance matrix is
4
