A Ground-based Optical Transmission Spectrum of WASP-6b

TL;DR: In this paper, a ground-based optical transmission spectrum of the inflated sub-Jupiter-mass planet WASP-6b was measured in 20 spectral channels from 480 nm to 860 nm using a series of 91 spectra over a complete transit event.

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.

Summary

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 exciting 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 transmission spectroscopy.
• 13 Fellow of the Swiss National Science Foundation.
• 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 scattering over the visible and near-infrared range, with the only detected feature being a narrow resonance line of Na (Pont et al. 2013).
• All of these studies have used highresolution spectroscopy or narrowband photometry to specifically target resonant lines of alkali elements.

2. OBSERVATIONS

• The transmission spectrum of WASP-6b was obtained performing 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.
• The seven comparison stars the authors used are listed in Table 1.
• 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 the authors used are shown for a typical exposure in Figure 1.
• 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.

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 estimated ci (one per row) was then subtracted from the central and middle regions in order to obtain a spectrum where only the stellar contribution remains.

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.
• 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.
• The authors 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 they are working with here.
• The authors 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.
• The authors checked the absolute deviation of each line from the fit and removed the most deviant one from their 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 procedure explained in the preceding paragraph served to wavelength-calibrate the first spectrum of the night closest in time to the NeHeAr lamps.

4. MODELING FRAMEWORK

• The largest departure from this idealized model in their observations will be given by systematic effects arising from atmospheric and instrumental effects, which are assumed to act multiplicatively on their signals.
• 15 Optimal extraction assumes that the profile along the wavelength direction is smooth enough to be approximated by a low-order polynomial.
• In particular, the authors 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.

4.1.1. Estimation of Systematic Effects via Principal Component Analysis of the Comparison Stars

• These perturbation signals have in common that they arise from the same physical and instrumental sources.
• The authors model this by assuming that a given perturbation signal is in fact a linear combination of a set of signals si(t), which represent the different instrumental and atmospheric effects affecting all of their light curves, i.e., Sk(t) = K∑ i=1 αk,isi(t).
• (2) Note that this model for the perturbation signal so far includes the popular linear and polynomial trends (e.g., si(t) = t i).
• A way of obtaining those signals is by assuming that the si(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.
• As written in Equation (3), the authors cannot separate si(t) from k(t), and in general the principal components will have contributions from both terms.

4.1.2. Selecting the Number of Principal Components

• The authors 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.
• The authors will determine the optimal number of components as the minimum number of components that are able to achieve the best predictive power allowed by the maximum set of N components available.
• K-fold cross-validation is a procedure that estimates prediction error, i.e., how well a model predicts out-of-sample data.
• This procedure is repeated for all folds.
• These l principal components, which are a set of light curves {si}li=1, are their estimates of the systematic effects, and the authors use the out-of-transit part of the light curve of WASP-6 as the validation data by fitting it with the {si}li=1.

4.2. Joint Parameter Estimation for Transit and Stochastic Components

• In the past sub-sections the authors 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 the authors have termed noise, i.e., the (t) term in Equation (1).
• Nonetheless, this is of no consequence as the authors just aim to model the residuals after the time series has been modeled with the {si}li=1.
• While the authors still call this term (t) in what follows, one should bear in mind this subtlety.
• In the present work the authors consider three different models: a white-noise model, where the covariance matrix is assumed to be diagonal; a flicker-noise model; and ARMA(p, q) models, where the structure of the covariance is determined via the parameters p and q (see Section 4.2.2 below for the definition of ARMA(p, q) models).

4.2.1. Flicker Noise Model

• Flicker noise is known to arise in many astrophysical time series (Press 1978).
• An efficient set of algorithms for its implementation in MCMC algorithms was proposed recently by Carter & Winn (2009).
• The basic idea of this implementation is to assume that the noise is made up of two components: an uncorrelated Gaussian process of constant variance and a correlated Gaussian process that follows this flicker-noise model.
• These two components are parametrized by σw and σr , characterizing the white and correlated noise components, respectively.
• A wavelet transform of the residuals takes the problem into the wavelet basis where flicker noise is nearly uncorrelated, making the problem analytically and computationally more tractable.

4.2.2. ARMA Noise Model

• ARMA models have been in use in the statistical literature for a long time with a very broad range of different applications (Brockwell & Davis 1991).
• Scargle 1981; Koen & Lombard 1993), these noise models have not been used so far for transit light curves to the best of their knowledge.
• In order to fit ARMA models to the residuals via an MCMC algorithm, the authors need the likelihood function of the model given that the residuals follow an ARMA(p, q) model.

4.2.3. Stochastic Model Selection

• Given the three proposed noise models for the stochastic signal (t), it remains to define which of the three affords a better description of the data, taking into account the trade-off between the complexity of the proposed model and its goodness of fit.
• There are several criteria for model selection; a comprehensive comparison between different criteria has been done recently by Vehtari & Ojanen (2012).
• The main conclusion is that, despite the fact that many model selection criteria have good asymptotic behavior under the constraints that are explicit when deriving them, there is no “perfect model selection” criterion, and there is a need to compare the different methods in the finite-sample case.
• Following this philosophy, the authors compare in this work the results of the AIC (“An Information Criterion”; Akaike 1974), the BIC (“Bayesian Information Criterion”; Schwarz 1978), the DIC (“Deviance Information Criterion”; Spiegelhalter et al. 2002), and the DICA, a modified version of DIC with a proposal for bias correction (Ando 2012).

5. LIGHT-CURVE ANALYSIS

• From the initial 10 comparison stars, only seven were used to correct for systematic effects.
• One star was eliminated on the grounds of having significantly less flux than the rest, and the other two due to not having the whole spectral range of interest recorded in the CCD.
• Given the seven comparison stars, the authors applied PCA to the mean-subtracted time series in order to obtain an estimate of the perturbation signals.
• The authors describe now the construction and analysis of the white-light transit light curve and the light curves for 20 wavelength bins.
• 17 ARMA(p, q) models have been considered recently in the modeling of radial velocity data (Tuomi et al. 2013).

5.1. White-light Transit Light Curve

• In order to obtain the white-light transit light curve of WASP-6, the authors summed the signal over the wavelength range 4718–8879 Å for the target and the comparison stars.
• The authors adopted priors for the white-light transit analysis are detailed in Table 2.
• The results of the MCMC fits assuming a white Gaussian noise model, an ARMA(2, 2) noise model, and a 1/f noise model for the residuals are shown in Figure 6, and a summary of the values of the information criteria for each of their MCMC fits is shown in Table 3.
• As opposed to deterministic components, the stochastic components cannot just be predicted given the times ti of the observations, as the authors only know the distribution of expected values once they know the parameters ({θ1, θ2, φ1, φ2, σw} for ARMA(2, 2), {σr, σw} for flicker noise, and σw for white Gaussian noise).
• The authors select the model parameters fitted using the 1/f noise model, which are quoted in Table 4, as the best estimates from now on.

5.2. WASP-6b Transmission Spectrum

• Priors were the same as the white-light analysis for parameters for μ1, σr , and σw, and the MCMC chains were set up similarly except that a thinning value of 103 was used.
• Boundaries that lie in the pseudo-continuum of the WASP-6 spectrum, as boundaries in steep parts of the spectrum such as spectral lines would in principle maximize redistribution of flux between adjacent bins under the changing seeing conditions that set the spectral resolution in their setup.
• For a given spectral bin, the number of principal components was selected separately because different systematics may be dependent on wavelength, and therefore the number of principal components needed may change.
• Figure 7 shows the baseline-subtracted data along with the best-fit transit model at different wavelengths, and Table 5 tabulates the transit parameters from the MCMC analysis for each wavelength bin.

5.3. Limits on the Contribution of Unocculted Stellar Spots

• As pointed out in several works (e.g., Pont et al. 2008; Sing et al. 2011b), stellar spots—both occulted and unocculted during transit—can affect the transmission spectrum.
• The decrease of flux during transit with respect to the out-of-transit flux F0 is given by (ΔF/F0) = k2 (neglecting any emission from the planet).
• The authors used the method described in Maxted et al. (2011) to look for periodic variations due to spots in the light curves of WASP6 from the WASP archive (Pollacco et al. 2006).
• The light curves typically contain ∼4500 observations with a baseline of 200 nights.

6. THE TRANSMISSION SPECTRUM: ANALYSIS

• The main feature of the transmission spectrum shown in Figure 8 is a general sloping trend with Rp/R∗ becoming smaller for longer wavelengths.
• The clear atmosphere models fail to give a better match to the spectrum due to the lack in the latter of evidence for the broad features expected around Na and K.
• The AIC for the scattering model assuming Gaussian noise given by the known error bars gives −115.3, while the values for the T = 1000 and T = 1500 clear atmosphere models are −97.2 and −90.9, respectively, providing a very significant preference for the scattering model.
• In order to assess the potential impact of correlations in the wavelength direction, the authors compute the partial autocorrelation function (PACF) for the residuals in the wavelength dimension.

7. DISCUSSION AND CONCLUSIONS

• The authors have measured the optical transmission spectrum for WASP-6b in the range ≈480–860 nm via differential spectrophotometry using seven comparison stars with IMACS on Magellan.
• The authors fit their transmission spectrum with three different models: two clear transmission spectra from Fortney et al. (2010) and a spectrum caused by pure scattering.
• The potentially prominent role of condensates or hazes in determining the transmission spectra of exoplanets has been apparent from the very first measurement (Charbonneau et al. 2002), and their transmission spectrum of WASP-6b is in line with what seems to be a building trend for transmission spectra with muted features in the optical.
• N.E. is supported by CONICYT-PCHA/Doctorado Nacional, and M.R. is supported by FONDECYT postdoctoral fellowship 3120097.

Figures

Table 2 MCMC Priors Used for the White-light Transit Analysis

Figure 4. Cross-validation error for the prediction of out-of-transit data using a different number of principal components with the 5-fold cross-validation procedure we adopted. Note that the minimum is at k = 7 (dashed lines indicate the value at the minimum and a value higher by 1σ ), but the value at k = 5 achieves similar error with lower degrees of freedom.

Figure 9. Predicted fractional change in k = (Rp/R∗) due to stellar spots that produce a rotation amplitude of = −4 mmag in the V band. The spots are assumed to have a temperature lower than that of the photosphere by ΔT = −500 K.

Table 1 List of Comparison Stars

Figure 1. Extracted spectra for WASP-6 and the seven comparison stars used in this work for a typical exposure.

Figure 5. Power spectral density of the residuals of the fit using white Gaussian noise (see Figure 6 to see the residuals). Note the preference for high power at small frequencies.

Table 3 Values for the Different Information Criteria (IC) for Each Noise Model Considered in Our MCMC Fits

Figure 6. Top: the circles show the baseline-subtracted light curves (i.e., light curves with the fitted perturbation signal subtracted) using the different noise models indicated. We also show the corresponding best-fit transit models (dashed line) and the best-fit transit models plus an estimate of the correlated noise component (solid line, only for the two rightmost light curves). The shaded regions indicate points that where used as out-of-transit data by the 5-fold cross-validation procedure that selected the number of principal components to use in the fits. Bottom: residuals between the best-fit transit model and the baseline-subtracted light curves (circles). The solid lines in the two rightmost set of points indicate estimates of the correlated components obtained by projecting the residuals into the best-fit correlated component model (see Section 5). The difference between the points and the solid lines (dashed line for the white Gaussian noise case) is the white Gaussian noise component, whose dispersion σw is indicated for each of the noise models considered and also illustrated with ±1σw bands.

Figure 8. Transmission spectrum of WASP-6b measured with IMACS.

Figure 7. Left: transits as observed in different wavelength channels along with the best-fit transit signal plus the stochastic 1/f noise signal. The obvious outlier close to t − t0 = 0.05 at the first, bluest, wavelength bin was not included in the fit. Right: residuals between the best-fit transit model and the baseline-subtracted light curves for each of the wavelength channels (circles). The solid lines indicate estimates of the correlated 1/f component obtained as described in Section 5. The best-fit parameters of the 1/f component are indicated over each set of residuals.

Figure 2. Raw light curves for WASP-6 and seven comparison stars used in this work as a function of time.

Table 4 WASP-6b Transit Parameters Estimated Using the White-light Transit Light Curve Using a 1/f -like Noise Model

Figure 10. Transmission spectrum of WASP-6b along with various models. Black dots with error bars indicate our measurements, while blue squares indicate the binned model of Fortney et al. (2010) with Teq = 1500 K, and red diamonds indicate the binned model using Teq = 1500 K. The green line and triangles indicate the best-fit line for a scattering model, which is the favored model in this case.

Table 5 Transit Parameters as a Function of Wavelength

Figure 3. Example of the structure that is expected in the power spectral density (PSD) of residual signals of the different types considered in this paper. The PSDs shown here are the mean of 10,000 realizations with the noise structures indicated in the figure. Note that the white-noise PSD is flat, while the flickernoise and the ARMA(0, 1) models cover low- and high-frequency ranges, respectively.

Full text

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