Table Of ContentAccepted for Publication in ApJ January 8th 2016
PreprinttypesetusingLATEXstyleemulateapjv.01/23/15
MARGINALISING INSTRUMENT SYSTEMATICS IN HST WFC3 TRANSIT LIGHTCURVES
H. R. Wakeford1,2
NASAGoddardSpaceFlightCenter,Greenbelt,MD20771,USA
D. K. Sing2 and T. Evans2
UniversityofExeter,Exeter,Devon,UK,EX44QL
D. Deming3
6 UniversityofMaryland,CollegePark,MD20742,USA
1
and
0
A. Mandell 1
2
n
Accepted for Publication in ApJ January 8th 2016
a
J
ABSTRACT
1
Hubble Space Telescope (HST) Wide Field Camera 3 (WFC3) infrared observations at 1.1-1.7µm
1
probe primarily the H O absorption band at 1.4µm, and has provided low resolution transmission
2
spectraforawiderangeofexoplanets. WepresenttheapplicationofmarginalisationbasedonGibson
]
P (2014)toanalyseexoplanettransitlightcurvesobtainedfromHSTWFC3,tobetterdetermineimpor-
E tant transit parameters such as Rp/R∗, important for accurate detections of H2O. We approximate
the evidence, often referred to as the marginal likelihood, for a grid of systematic models using the
.
h Akaike Information Criterion (AIC). We then calculate the evidence-based weight assigned to each
p systematic model and use the information from all tested models to calculate the final marginalised
- transitparametersforboththeband-integrated, andspectroscopiclightcurvestoconstructthetrans-
o
mission spectrum. We find that a majority of the highest weight models contain a correction for a
r
t linear trend in time, as well as corrections related to HST orbital phase. We additionally test the
s
a dependence on the shift in spectral wavelength position over the course of the observations and find
[ thatspectroscopicwavelengthshiftsδ (λ),bestdescribetheassociatedsystematicinthespectroscopic
λ
lightcurves for most targets, while fast scan rate observations of bright targets require an additional
1
level of processing to produce a robust transmission spectrum. The use of marginalisation allows
v
for transparent interpretation and understanding of the instrument and the impact of each system-
7
atic evaluated statistically for each dataset, expanding the ability to make true and comprehensive
8
comparisons between exoplanet atmospheres.
5
2 Subject headings: methods: data analysis, techniques: spectroscopic, planets and satellites: atmo-
0 spheres
.
1
0
1. INTRODUCTION solarcase,indicatingthatthetrendobservedinthemetal
6
abundanceofthesolarsystemgiantplanetatmospheres,
1 H2Oisthemostspectroscopicallydominantspeciesex- i.e. decreasingmetalenhancementwithincreasingmass,
: pected in hot Jupiter atmospheres, and a key molecules
v extends out to exoplanet atmospheres. Wakeford et al.
for constraining atmospheric compositions. In most
i (2013)werealsoabletoobtainastrongdetectionofH O
X lower atmosphere models of hot Jupiters, H2O is well in the atmosphere of the hot Jupiter HAT-P-1b. T2he
mixedthroughouttheatmosphere,andmostfeaturesbe-
r analysis made the use of a common spectral type star in
a tween 0.7 and 2.5µm are the result of H2O vibrationan- the field of view of the WFC3 detector and performed
rotation bands (Brown 2001). Hubble Space Telescope
differential spectrophotmetry using the companion spec-
(HST) Wide Field Camera 3 (WFC3) infrared observa-
trum, similar to methods used from the ground. The
tions at 1.1-1.7µm probe primarily the H O absorption
2 robust nature of the transmission spectral shape shown
band at 1.4µm, and has provided low resolution trans-
for analysis methods with and without the use of a com-
missionspectraforawiderangeofexoplanets(e.g. Berta
panion spectrum, shows that the emphasis needs to be
etal.2012;Gibsonetal.2012;Demingetal.2013;Wake-
placed on methods that not only effectively reduce the
ford et al. 2013; Knutson et al. 2014; Stevenson et al.
data, but ones that can be applied successfully to multi-
2014; ?). Using both transmission and emission spectra
pledatasetssimultaneouslyforatruecomparativestudy.
Kreidberg et al. 2014 measured the relative H O abun-
2 A wide range of reduction and analysis methods have
dance in the atmosphere of WASP-43b. They compare
been used on HST WFC3 transit datasets, making com-
the measured H O abundance with the solar system gi-
2 parisons between datasets and planetary atmospheres
antplanetsusingthemetallicityexpectedfora‘broadly’
difficult. Studies of multiple WFC3 datasets have been
conducted across multiple programs and observation
[email protected]
2 Wakeford et al.
modes. These have attempted to define a common de- riesobservationsisdependentonrobustlyaccountingfor
trendingtechniquetoapplytoalldatasets(Demingetal. anysystematiceffectsfromthetelescopeanddetectorfor
2013; Mandell et al. 2013; Ranjan et al. 2014; Kreidberg both ground and space based instruments. Each of the
et al. 2014, Kreidberg et al. 2014; Haynes et al. 2015). modelsusedintheliteratureattempttocorrectthearray
WFC3hastwomainobservingmodesthatarecommonly of systematics observed independently in each dataset.
used for transiting exoplanet spectra; stare and spatial However, not all datasets display the same combination
scan mode. Stare mode has been used for a majority of systematics, which appear to be highly dependent on
of HST observations and is useful when observing dim- the observational mode and set-up.
mer target stars where the photon counts/pixel/second Inthispaperweapplythemarginalisationmethodpro-
is low; observing brighter targets in this mode leads to posedbyGibson(2014)tonotonlyincorporatetheanal-
saturation. Stare mode maintains a constant pointing ysis of multiple systematic treatments on the lightcurve
ofthetelescopethroughouttheobservation,maintaining butalsotobeapplicabletomultipledatasetsallowingfor
thesamepixelpositiononthedetector. Spatialscanning a cross-comparison between different transmission spec-
mode was made available on WFC3 in Cycle 19 (2012) tra. In §2 we outline the current methods used in the
and is now implemented as the main mode of observa- literature to reduce and analyse current WFC3 transit
tions for transiting exoplanets, as targets observed for datasets. In§3weputforwardanddiscussthemarginal-
atmospheric follow-up with such instruments often orbit isation analysis technique, which incorporates a number
brighter target stars (V magnitudes brighter than 11). ofpreviouslypublishedmethodstoproducerobusttran-
WFC3 spatial scanning involves nodding the telescope sit parameters. We the discuss marginalisation in the
during an exposure to spread the light along the cross- context of our observations in §4 and highlight the key
dispersion axis, resulting in a higher number of photons aspects of marginalisation for transit datasets. We then
byafactoroftenperexposurewhileconsiderablyreduc- contrast this technique in §5, to assess the impact of
ing overheads. This also increases the time of saturation different analysis methods and on the computed transit
of the brightest pixels, which allows for longer exposure parameters.
times (McCullough 2011).
Mandelletal.(2013)conductedthefirstreanalysistest
2. SYSTEMATIC MODEL CORRECTIONS
of WFC3 data for WASP-19b and WASP-12b, with the
additionofthefirstanalysisofWASP-17b. Thereanaly- A transit lightcurve consists of N stellar flux mea-
sis incorporated a wavelength dependent systematic cor- surements observed at time t, collectively referred to as
rection over the methods used in the previous analysis. the data D or the lightcurve. To model each of the
Whilethisstudyproducedalmostphoton-limitedresults lightcurves we calculate a Mandel & Agol (2002) tran-
inindividualspectralbins,thespectralfeaturesobserved sit model T(t,λ) following a non-linear limb-darkening
in the transmission spectra were degenerate with vari- law as defined in Claret (2000). The model is then fit
ous models of temperature and compositions making in- tothedatabyallowingthebaselineflux F∗, Rp/R∗, and
terpretation difficult. Ranjan et al. (2014) conducted a centre of transit time to vary. We fixing the other plan-
study of four stare mode WFC3 transit lightcurves from etary system parameters such as inclination and a/R∗
the large HST program, GO12181. However, they were from previously published values as the phase coverage
unable to resolve any features in the transmission spec- of HST lightcurves do not well constrain these parame-
trum for three of the planetary atmospheres, and were ters on their own. The final form of the model fit to the
unabletoextractarobusttransmissionspectrumforone data becomes,
of the datasets, as different treatments to the data gave
F(λ,t)=F ×T(t,λ)×S(t,λ) (1)
moderately different results. ∗
Additionally, the analysis of the very hot Jupiter where S(t,λ) is the sytematics model normalised to
WASP-12b which was observed as part of GO12230, unity.
P.I. M. Swain is a good example of differences caused One of the issues encountered when analysing obser-
from analysis techniques. WASP-12 was observed in vational datasets is determining the impact that instru-
stare mode using WFC3 G141 slitless grism which con- mental systematics have on the resultant measurements.
tained the spectrum of both the target planetary host SincetheadventofWFC3’sapplicationtotransitingex-
star and the M-dwarf binary companions to WASP-12, oplanets, a number of systematic models have been used
which overlapped in the spectral response on the detec- toreduceG141spectroscopicdata(e.g. Bertaetal.2012;
tor. The most recent reanalysis of this data explores Gibsonetal.2012;Wakefordetal.2013;Lineetal.2013;
theeffectofsystematicmodelcorrectionontheabsolute Deming et al. 2013; Stevenson et al. 2014).
transit level measured for the atmospheric transmission Figure1showsthreeexamplesofsystematictrendsob-
spectrumfindingthattheabsolutetransitdepthissensi- servedinWFC3transitlightcurves: “HSTbreathing”ef-
tivetothesystematicmodelappliedtothedata(Steven- fects, visit-long slopes, and the ‘ramp’ effect. The“HST
son et al. 2014). breathing” effect shown in Fig. 1a displays a periodic
Spectroscopically resolving transmission absorption systematic across each orbit of data. This is attributed
features in transit lightcurves is important for determin- to the known thermal variations which occur during the
ingcompositionalinformationofexoplanetatmospheres. orbit of HST as it passes into and out of the Earths
Unlike Spitzer’s IRAC instrument, which is able to ob- shadow,causingexpansionandcontractionofHST.This
tain photometric points in H O absorption bands, HST can be most easily seen in the middle panel in relation
2
WFC3isvitallyabletospectroscopicallyresolvefeatures to the HST orbital phase. The “HST breathing” effect
in the atmospheres of exoplanets. Accurately determin- systematic has been noted and corrected for in multi-
ing the parameters of a transiting planet via time se- ple datasets with a variety of parametrised models (e.g.
Marginalising instrument systematics in HST WFC3 transit lightcurves 3
XO-1 HD 209458 WASP-17
Fig. 1.— Three of the main systematic effects observed in HST WFC3 transit datasets: a) “HST breathing” effect caused by the
temperaturevariationsintheorbitalperiodofHST.b)Visit-longSlope,alineartrendobservedacrosstheentireobservingperiodforall
transitlightcurvesobservedwithHSTWFC3. c)The‘ramp’effectwhichisthoughttobecausedbychargetrappingbetweenbufferdumps.
Lowerpanelsshowtheresidualsofeachdatasetwithrespecttodifferenttimeseriesparameters. Top: residualsintermsofplanetaryphase.
Middle: residualsintermsofHSTorbitalphase,whereeachHSTorbitofdataisoverplottedonsubsequentorbits. Bottom: residualsin
termsofexposurenumber.
Wakeford et al. 2013; Line et al. 2013; Stevenson et al. each group of exposures obtained between buffer dumps
2014) which remove systematics based on functions of referredtoasthe‘ramp’or‘hook’effect(e.g. Bertaetal.
the HST orbital time period and phase. 2012; Mandell et al. 2013). The dataset of WASP-17 in
Many groups have also reported a visit-long trend in Fig. 1cclearlydisplaysthiseffect,andthebottomresid-
WFC3 lightcurves. This trend can be seen clearly in the ual plot in terms of exposure number shows the highly
raw band-integrated lightcurve of HD209458 shown in repeatable aspect of the systematic. This is thought to
Fig. 1b, which displays a significant slope across the en- be caused by charge trapping on the detector and it has
tire observation period. This can be seen in relation to been found that the ‘ramp’ is, on average, zero when
both planetary phase and exposure number. This sys- the count rate is less than about 30,000 electrons per
tematic trend has not been correlated with any other pixel (Wilkins et al. 2014). This is commonly seen in
physical parameter of WFC3 observations. However, it stare mode observations where the count rapidly builds
has been shown to significantly affect the resultant sys- up over a small pixel range.
tem parameters obtained from the lightcurve. Strong wavelength dependent shift of the stellar spec-
In addition to orbital phase trends, both in planetary trum across the detector throughout the course of the
and HST space, a number of lightcurves have been dom- observation can also affect the spectroscopic lightcurves
inated by a systematic increase in the intensity during and therefore the measured transit parameters (Deming
4 Wakeford et al.
TABLE 1
TableofexoplanetsobservedusingWFC3G141grism,modethedataisobservedin,thesystematicmodelusedtode-trend
the transmission spectral data, author, HST program, and year of observation.
HST program PI Cycle Planet Year Mode Model Paper
GO11740 F. Pont Cycle 17 HD189733 2010 stare GP Gibson et al. (2012)
GO12181 D. Deming Cycle 18 WASP-17 2011 stare θ+doot+δ Mandell et al. (2013)
λ
HD209458 2012 scan θ+δ (λ).A(λ) Deming et al. (2013)
λ
XO-1 2011 scan θ+δ (λ).A(λ) Deming et al. (2013)
λ
HAT-P-12 2011 stare θ+eψ/τ +φ Line et al. (2013)
WASP-19 2011 stare θ+doot Huitson et al. (2013)
θ+doot+δ Mandell et al. (2013)
λ
WASP-4 2010 stare doot Ranjan et al. (2014)
TrES-2 2010 stare doot Ranjan et al. (2014)
TrES-4 2010 stare θ+doot Ranjan et al. (2014)
CoRoT-1 2012 stare θ+doot Ranjan et al. (2014)
GO12230 M. Swain Cycle 18 WASP-12 2011 stare θ Swain et al. (2013)
θ+doot+δ Mandell et al. (2013)
λ
θ+φN Sing et al. (2013)
θ+eψ/τ +φ Stevenson et al. (2014)
GO12251 Z. Berta Cycle 18 GJ1214 2010 stare doot Berta et al. (2012)
GO12449 D. Deming Cycle 19 HAT-P-11 2012 scan θ+δ (λ).A(λ) Fraine et al. (2014)
λ
GO12473 D. Sing Cycle 19 WASP-31 2012 scan θ+φN +δN ?
λ
HAT-P-1 2012 scan θ+φN Wakeford et al. (2013)
GO12881 P. McCullough Cycle 20 HD189733 2013 scan θ McCullough et al. (2014)
GO13021 J. Bean Cycle 20 GJ1214 2012-13 scan θ+eψ/τ +φ Kreidberg et al. (2014)
GO13064 D. Ehrenreich Cycle 19 GJ3470 2013 stare doot Ehrenreich et al. (2014)
GO13338 K. Stevenson Cycle 21 GJ436 2013 scan θ2+eφ+φ Stevenson et al. (2014)
GO13467 J. Bean Cycle 21 WASP-43 2013 scan θ+eψ/τ +φ Kreidberg et al. (2014)
GO13501 H. Knutson Cycle 20 HD97658 2014 scan θ+eψ/τ Knutson et al. (2014)
GP - Gaussian Process; doot - Divide-oot; θ - Visit-long slope; φ - HST phase; δ - wavelength shift; eψ/τ - model ramp
λ
et al. 2013). Large shifts in the wavelength direction ramp in the form of the equation
on the WFC3 detector may introduce sub-pixel to pixel
F
sized variations in each spectroscopic bin when dividing orb =(C+Vθ+Bφ)(1−Reψ/τ), (2)
the spectrum into wavelength channels. This is likely Fcor
the result of fast scanning rates used for very bright tar-
where F /F are the lightcurve residuals, θ is the
gets where additional motion during a rapid scan rate orb cor
planetary phase, φ is the HST orbital phase, and ψ is
can introduce variations between spectral exposures, as
the phase over which the ramp feature occurs, which ac-
fastscanratesspreadthestellarspectrumacrossalarger
counts for the visit-long slope V, the orbit-long slope B,
range of the detector making it much harder for the fine
and a vertical offset C applied to the whole lightcurve.
guidance system to hold a fixed wavelength position on
The exponential model for the ramp has an additional
the detector.
two parameters: the ramp amplitude R, and the ramp
A combination of observational and phase dependent
timescale τ. This method is used by a number of groups
instrument systematics has been observed in all WFC3
when a ramp is observed in the raw band-integrated
transit datasets. We outline the main systematic mod-
lightcurve (Kreidberg et al. 2014; Kreidberg et al. 2014;
els used to correct for these effects, with the full list of
Knutson et al. 2014).
published systematic correction models presently used
Incaseswheretheorbitaltimescalematchesthephase
outlined in Table 1.
overwhichtherampoccursasimplificationoftheramp-
model can be used, as seen in Stevenson et al. (2014),
2.1. Exponential model ramp S(t,λ)=[1+r0θ+r1θ2]×[1−er2φ+r3 +r4φ], (3)
ThefirstmethodputforwardtocorrectWFC3system- where r are free parameters and the phase ψ over
0−4
atics, as an analytical model whose parameters attempt which the ramp feature occurs is now equal to the HST
torepresentthephysicalprocessesoftheinstrument,was orbitalphaseφ. Thisisdisplayedas“θ+eφ+φ”inTable
outlinedbyBertaetal.(2012)definedastheexponential 1. Fraine et al. (2014) find that a simple linear visit-
model-ramp. Theexponentialmodelsapplyanexponen- long correction and a single exponential ramp in HST
tial ramp over sets of exposures, corrects for orbit-long phase can correct for the systematics observed in the
andvisit-longslopes,andisintendedtomodelthecharge transit lightcurve of the hotNeptune HAT-P-11 without
trapping. Lineetal.(2013)showtheexponentialmodel-
Marginalising instrument systematics in HST WFC3 transit lightcurves 5
the need for a squared term in time. sure within an orbit and divides the in-transit orbits
Both of these methods rely on the timescales of each by the template created. This requires each of the in-
ramp to be the same in each orbit and for each orbit to transit exposures to be equally spaced in time with the
have the same repeating systematic. This therefore de- out-of-transitexposuresbeingusedtocorrectthem,such
pends heavily on the scheduling of each exposure within that each corresponding image has the same HST phase
aHSTorbit, whichanumberofHSTWFC3datasetsdo andthatadditionalsystematiceffectsarenotintroduced.
not meet. While this does not rely on knowing the relationship be-
tween measured photometry and the physical state of
2.2. Polynomial models thecameraitdoesrequiretheretobeanevennumberof
An additional polynomial method to correct for “HST exposures equally spaced from orbit to orbit where the
breathing” effects in the data is discussed in Wakeford exposures occupy the same HST phase space. This is so
et al. (2013), which assumes that the effects fit a poly- thatthesystematicsinducedbytheknown“HSTbreath-
nomialfunctionratherthanbeingexponentialinnature. ing” trend caused by temperature variations in its orbit
Thismethodalsoseekstoremovethevisit-longslopeob- can be effectively eliminated. The divide-oot method
served in each WFC3 transit dataset using a linear time relies on the cancellation of common-mode, wavelength
trend in planetary phase in addition to the HST phase independent, systematic errors by operating only on the
corrections, data themselves using simple linear procedures, relying
ontrendstobesimilarinthetimedomainoveranumber
n of orbits. This is listed as ‘doot’ in Table 1.
(cid:88)
S(t,λ)=T θ+ p φi (4)
1 i
i=1 2.4. Spectral template corrections
where θ is planetary phase representing a linear slope AsomewhatsimilartechniquewasadoptedbyDeming
over the whole visit when multiplied by the free or fixed etal.(2013)andMandelletal.(2013)fortheiranalysisof
parameter T1, φ is HST phase accounting for “HST WFC3datarelyingoncommontrendsinthewavelength
breathing” effects when multiplied by p1−n, which are domain. Toaccountforthissub-pixelchangeinthespec-
eitherfreeparametersorfixedtozerotofitthemodelto tral wavelength solution for the stellar spectrum across
the data. each exposure, Deming et al. (2013) introduced a spec-
In addition to “HST breathing” trends, functions in traltemplatetechniquetotheextractedspectrumofeach
wavelength shift on the detector were needed to correct image. The template spectrum is constructed from the
for systematics seen in the lightcurve of WASP-12 (Sing average observed spectrum from exposures within one
et al. 2013). The systematic model then takes the form, hourof1stand4thcontact. Thetemplateisthenusedto
fitthewavelengthsolutiontoeachexposurespectrumby
n n
(cid:88) (cid:88) stretchingandshiftingitinwavelength,steppingthrough
S(t,λ)=T θ+ p φi+ l δj (5)
1 i j λ small increments to determine the best fit solution. The
i=1 j=1 individualspectraarethendividedbythetemplatespec-
trum to create residuals effectively removing any addi-
where δ is the array of the shift in the wavelength (x)
λ
tionalbackgroundcontributionsandthewavelengthshift
direction on the detector over the visit, and l are
1−n
on the detector. This technique (‘stretch and shift’) also
fixed to zero or free parameters, similar to that used for
allowsforthecancellationofcommon-modesystematics,
the “HST breathing” correction. δ is created by cross-
λ
similar to the divide-oot method, however, this requires
correlating the stellar spectrum for each exposure with
the common-mode systematics to be consistent in wave-
a template stellar spectrum across the whole wavelength
length rather than in time. This technique is labeled
range. This is shown in Table 1 by a combination of the
different systematics corrected for e.g. θ+φN +δN. as ‘δλ(λ)’ in Table 1. However, this method only pro-
λ duces a relative depth measurement for the atmospheric
Importantly, this model does not require each of the
signatures and does not provide absolute planet-to-star
orbits to have the same number of exposures, or consis-
radius ratio values needed for comparative studies and
tent repeating systematics in each orbit or between each
combined datasets.
buffer-dump making it a robust method to apply to any
transit dataset from WFC3. However, due to the large
number of potential free parameters in each systematic 3. MARGINALISATION
model, and the potential to go to high orders of poly- Thus far it is not clear which of the systematic mod-
nomial for each systematic, using this method can po- els that have been employed to correct WFC3 transit
tentially introduce additional systematics if the correct lightcurves is the best for each individual dataset. How-
model is not initially chosen. ever, it is clear that none of these simple systematic
models will work for all datasets. We attempt to rec-
2.3. Divide out-of-transit
tify this by performing a marginalisation over a grid of
Berta et al. (2012) also suggested a separate method, systematic models which incorporate corrections for the
calleddivide-oot, forcorrectingthesystematic‘ramp’or differentsystematicsobservedacrosstheexoplanettran-
‘hook’observedinanumberofdatasets. Thedivideout- sit datasets. In this way, model selection does not have
of-transitmethod(divide-oot)reliesonthehooksystem- such a large influence on the final result.
atic being “extremely” repeatable between orbits, HST Use of the Bayesian Information Criterion (BIC) to
phase, in a visit. select between different systematic models, which takes
Divide-oot uses the out-of-transit orbits to compute Occam’s Razor into effect by penalising models with in-
a weighted average of the flux evaluated at each expo- creasing complexity, has been previously used for WFC3
6 Wakeford et al.
TABLE 2
Table of all parametrised systematic models applied to the lightcurves showing the combination of visit-long trends as a
function of planetary phase (θ), functions of HST orbital phase (φ), and functions dependent on wavelength shifts (δ ) in
λ
the data. In addition to these models we apply the two exponential orbital phase models outlined in Stevenson et al.
(2014) .
No. θ φ φ2 φ3 φ4 δ δ2 δ3 δ4 No. θ φ φ2 φ3 φ4 δ δ2 δ3 δ4
λ λ λ λ √ λ λ λ λ
0 25
√ √ √
1 26
√ √ √ √ √
2 27
√ √ √ √ √ √ √
3 28
√ √ √ √ √ √ √ √ √
4 29
√ √ √
5 30
√ √ √ √ √
6 31
√ √ √ √ √ √ √
7 32
√ √ √ √ √ √ √ √ √
8 33
√ √ √ √ √ √ √ √ √ √ √
9 34
√ √ √ √ √
10 35
√ √ √ √ √ √ √
11 36
√ √ √ √ √ √ √ √ √
12 37
√ √ √ √ √ √ √ √ √ √ √
13 38
√ √ √ √ √ √ √ √ √ √ √ √ √
14 39
√ √ √ √ √ √ √
15 40
√ √ √ √ √ √ √ √ √
16 41
√ √ √ √ √ √ √ √ √ √ √
17 42
√ √ √ √ √ √ √ √ √ √ √ √ √
18 43
√ √ √ √ √ √ √ √ √ √ √ √ √ √ √
19 44
√ √ √ √ √ √ √ √ √
20 45
√ √ √ √ √ √ √ √ √ √ √
21 46
√ √ √ √ √ √ √ √ √ √ √ √ √
22 47
√ √ √ √ √ √ √ √ √ √ √ √ √ √ √
23 48
√ √ √ √ √ √ √ √ √ √ √ √ √ √ √ √ √
24 49
50 θ ×(1-eφ)+φ 51 (θ +θ2)×(1-eφ)+φ
Fig. 2.—TheevidencebasedontheAICplottedagainstthesystematicmodelnumbercorrespondingtoTable2. Thisisanexampleof
theevidencecomputedforWASP-31wherethebest-fittingsystematicmodelandtherawlightcurveevidenceareindicatedwitharrows.
observations (e.g. Wakeford et al. 2013; Sing et al. 2013; statistical likelihood of each model. This ensures that a
?;Haynesetal.2015),andhasbeenappliedtoarangeof varietyofsystematicmodelsaretakenintoaccountwhen
datasets (e.g. Huitson et al. 2013; Crossfield et al. 2013; measuring the R /R without having to choose between
p ∗
Nikolovetal.2014). Herewetakethisonestepfurtherby models.
adoptingthemethodologyproposedbyGibson(2014)of Before marginalisation, the overall systematic models
marginalising over multiple systematic models to calcu- that are going to contribute to the final weighting must
laterobusttransitparameters. Thisallowsustoquantify bedecidedupon. Weuseagridof49polynomialmodels
thedegeneracybetweenourphysicalparametersofinter- whichincorporateallknownidentifiedsystematictrends
est and our choice of systematic model. In this case we in the data (see Table 2), with the addition of the two
want to determine the value and associated uncertainty exponential models outlined by Stevenson et al. (2014),
the planet-to-star radius ratio (R /R ) for each of our and the uncorrected lightcurve. All 52 models are then
p ∗
lightcurves after correcting for the systematics inherent placed into an array of varying free parameters to be fit
in the data. For each systematic model used to correct or fixed in turn and looped over for each lightcurve fit
the data we calculate the evidence of fit, defined as the to calculate the weighting assigned to each in turn. We
probability that the data would be produced given the assume equal priors on all the systematic models tested,
systematic model, which is then used to apply a weight where no model is preferentially preferred over another.
to the parameter of interest measured using that model. It is also important to note that marginalisation re-
A weighted average of R /R is then calculated which lies on the fact that at least one of the models being
p ∗
takes into account the individual weights of each fit and marginalised over is a good representation of the sys-
Marginalising instrument systematics in HST WFC3 transit lightcurves 7
tematics in the data. Our grid of parametrised mod- sation to fit the data. We then determine the best sys-
els includes all combinations of factors up to the fourth tematic model used based on the AIC and rescale the
order in both HST phase, to correct for “HST breath- lightcurve uncertainties such that it has a reduced χ2
ing” effects, and up to the fourth order in wavelength of 1. We then re-run each of the systematic models with
shift, in addition to the visit-long linear trend noted by thelightcurvespriortotomarginalisation. Typicallythis
all groups. By also incorporating the Stevenson et al. rescales the errors by a factor of ∼1.1 to ∼1.2 times the
(2014) exponential HST phase models, with a linear and theoretical photon noise limit of WFC3. Once we have
squared planetary phase trend, we make the assumption applied this inflation to our errorbars the approximated
that this condition of marginalisation is satisfied. Under evidence is modified to incorporate the likelihood func-
thisapproach,weeffectivelyaveragetheresultsobtained tion.
from a suite of systematic models in a principled man- Toapplythistoourdatasetweexpandtheabovelike-
ner. In doing so, we marginalise over our uncertainty lihood function,
as to which systematic model is actually the “correct”
(cid:34)N (cid:35)
mhaovdeeel.quTahlliysiwseelslpfietctiianlglysyimstpeomratatincts,wahseinssoefvteenratlhmeocdaseels. ln[P(D|α∗,Sq)]=ln (cid:89) √1 e−σri22 (8)
σ 2π
i=1
3.1. Evidence and weight
whereσ istheuncertaintyonthedata,andr represents
i
To calculate the weighting assigned to each of the sys- the model residual for the ith datapoint.
tematic models and subsequently the final marginalised
ppalarnaemtaerteyrtrfaonrsitth,eweplfiarnsetth-taov-esttaordreatedrimusinreatthioeefovrideeancche ln[P(D|α∗,Sq)]=(cid:88)N ln(cid:20) √1 e−2rσi22(cid:21), (9)
σ 2π
thateachsystematicmodelhaswhenfittothedata. The
i=1
evidenceE offitassignedtoeachsystematicmodelS is
q q
givenbytheprobabilityofthedataD giventhemodelq =(cid:88)N ln(cid:20) √1 (cid:21)− 1(cid:16)ri(cid:17)2, (10)
andisoftenreferredtoasthemarginallikelihood. Inthe
σ 2π 2 σ
absenceofaccuratepriorsonwhichtoplacealikelihood, i=1
we can use an approximation for the evidence (Gibson (cid:16) (cid:17) 1
2014), such as, =−Nln σ(2π)21 − χ2, leading to (11)
2
1 1
lnE =lnP(D|S )≈− BIC=ln[P(D|α ,S )]− MlnN, N 1
q q 2 ∗ q 2 (6) ln[P(D|α∗,Sq)]=−Nlnσ− 2 ln2π− 2χ2 . (12)
where the BIC is the Bayesian Information Criterion, Substituting Eq 11 into Eq 6, we arrive at,
which is equated to the logarithmic probability of
1 1
the data given the parameter and systematic model lnE =−Nlnσ− Nln2π− χ2−M (13)
(ln[P(D|α ,S )]) minus the number of free parameters q 2 2
∗ q
M multiplied by the log number of data points being fit This gives us the final form of the evidence function
N. for each of our systematic models applied to the data.
The BIC is the most commonly used criterion to se- This now needs to be transformed into a weighting so
lect between models in the current exoplanet literature. that each of the systematic models is assigned a per-
Alternatively to the BIC, the Akaike Information Crite- centage of the overall probability and, when normalised,
rion (AIC) can be calculated, which does not penalise (cid:80) P(S |D)=1.
themodelasstronglyforaddedcomplexitygivenalarge q q
The individualweight(W ) foreachsystematicmodel
number of data points, q
is calculated by
1
lnEq =lnP(D|Sq)≈−2AIC=ln[P(D|α∗,Sq)]−M. (cid:88)Nq
W =P(S |D)=E / E . (14)
(7) q q q q
BoththeAICandBIChavestrongtheoreticalfounda- q=0
tions and can be used for model selection (Burnham & Where N is the number of models fit, α is the
q m
Anderson 2004). As the number of data points in each marginalised parameter, and α is the measured param-
q
datasetgreatlyexceedsthatofthenumberoffreeparam- eter for each model.
eters in our most complex model we choose to minimise The weighting assigned to each model due to the ev-
the AIC to give our best-fitting model with the largest idence parameter can then be used to calculate the
evidence. This is also favoured in Gibson (2014) as it weighted mean of all the parameters (α) of interest
allowsformoreflexiblemodelsintothelikelihood,which
typically leads to more conservative error estimates. Nq
(cid:88)
Theevidencecalculatedforeachmodeladditionallyre- α = (W ×α ), (15)
m q q
lies upon the uncertainty placed on the data (σ), which
q=0
is dominated by photon noise in spectral extraction
pipelines. To ensure appropriate uncertainties, we start and the uncertainty (σα) on that parameter can be de-
by applying photon noise errorbars to each point and terminedfromσαq i.e. theuncertaintyontheparameter
including our systematics models when running MPFIT
(Markwardt2009)whichusesL-Mleast-squaresminimi-
8 Wakeford et al.
HAT-P-1 WASP-31
XO-1 HD 209458
Fig. 3.— The evidence based on the AIC plotted against the
evidencebasedontheBIC(top-squares),theχ2(middle-triangles),
theRp/R∗(bottom-circles)forthe52modelscomputedforWASP-
31b. Thebest-fittingmodelisthatwiththehighestAICevidence
parameter, while a similar trend is seen by reducing the ∆χ2 the
AICpenalisesthemodelsforaddedcomplexity.
WASP-17
α determined from the qth model, Fig. 4.— Single exposure ima output image for each observed
(cid:118) exoplanet host star. Each spectral exposure is outlined in green
(cid:117) N with the center of the spectrum marked with a dotted line. The
(cid:117)(cid:88)
σ(α)=(cid:116) (Wq[(αq−αm)2+σα2q]). (16) rliengei.onusedforbackgroundsubtractionisboxedbyadashedyellow
q=0
This assumes that the AIC is a good approximation for 4. OBSERVATIONS
the evidence calculated for each systematic model and Observationsofsingletransiteventsareconductedus-
uses a Gaussian approximation for the posterior distri- ing Hubble Space Telescope’s (HST) Wide Field Cam-
butions (?). era 3 (WFC3) in the infrared (IR) with the G141 spec-
Fromthemarginalisationoverall52systematicmodels troscopic grism. Our observations span two HST large
ontheband-integratedlightcurve,thebest-fittingmodel programs and two observation modes, stare mode and
canrevealthedominantcontributingsystematicstoeach spatial scan mode, acquired between 2011–2012 over 25
dataset. Figure2showsthecalculatedevidencebasedon HST orbits.
the AIC approximation for all 52 models when fit to an Table 3 outlines the observational and planetary
example dataset, where the systematic model with the system parameters for each of the exoplanet hot
highest overall weighting for the marginalisation is indi- Jupiters studied here: HAT-P-1b, WASP-31b, XO-1b,
catedbyanarrowalongwiththerawlightcurvefit. This HD209458b, and WASP-17b. Each of these exoplanet
clearly shows the sample of systematic models that are transmission spectra has been previously analysed and
favouredwhencorrectingthisdatasetandemphasisesthe published: HAT-P-1b(Wakefordetal.2013),WASP-31b
need for through systematic model analysis. It is some- (?), XO-1b (Deming et al. 2013), HD209458b (Deming
timesthecasethatanumberofmodelswillhaveastrong etal.2013),andWASP-17b(Mandelletal.2013). Across
weightingontheband-integratedlightcurve, wheretheir these five exoplanet transmission spectra there are a va-
weight>10%. Often these models correct for the same riety of measured features over the expected H O ab-
2
combination of systematic trends assigning a different sorption bands, from full amplitude features extending
order to the polynomial used for correction to within 1 several scale heights in the atmosphere, to muted and
order. absent features. The study here is intended to demon-
In Fig.3 we demonstrate the statistical correlation be- stratethemostuniformanalysisandcomparisontodate,
tween different factors used to select models, comparing such that differences between spectra can be more eas-
the evidence based on the AIC approximation to that of ily tied to the planets themselves and changes between
the evidence based on the BIC approximation, and to different reduction techniques can be minimised.
the ∆χ2 for all 52 systematic models. We additionally We use the “ima” outputs from WFC3’s Calwf3
show the difference in the measured Rp/R∗ computed pipeline. For each exposure, Calwf3 conducts the fol-
for each model relative to the weight applied to that fit, lowing processes: bad pixel flagging, reference pixel sub-
again demonstrating the importance of using the correct traction,zero-readsubtraction,darkcurrentsubtraction,
systematic model when correcting transit data. non-linearity correction, flat-field correction, and gain
Marginalising instrument systematics in HST WFC3 transit lightcurves 9
TABLE 3
Table of the observation parameters for the five planetary transits measured with WFC3. The planetary system
parameters used for each of the datasets is also listed along with the band-integrated limb-darkening parameters
calculated using 1D Kurucz stellar models.
HAT-P-1 WASP-31 XO-1 HD209458 WASP-17
GO Program 12473 12473 12181 12181 12181
Date 2012–07–05 2012–05–13 2011–09–30 2012–09–25 2011–07–08
Mode Scan Scan Scan Scan Stare
NSAMP 4 8 9 5 16
Subarray size 512 256 128 256 512
Exposure Time (s) 46.695 134.35 50.382 22.317 12.795
Peak Pixel Count 25,000 38,000 38,000 44,000 64,000
No. Exposures 111 74 128 125 131
Scan Rate (pix/s) 1.07 0.15 0.43 7.43 -
Rp (R ) 1.319 1.55 1.209 1.38 1.932
J
Mp (M ) 0.525 0.48 0.942 0.714 0.477
J
T (K) 1322 1570 1210 1459 1755
eff
Period (days) 4.46529976 3.405886001 3.94163 3.52474859 3.7354833692
Epoch (MJD) 53979.43202 56060.69042 55834.3419 56195.7595 55750.2973239
inclination (◦) 85.634 84.670 88.92 86.59 86.917160
a/R* 9.91 8.19 11.24 8.859 7.03
T* (K) 5980 6250 5750 6117 6550
[Fe/H]∗ 0.130 -0.2 0.02 -0.02 -0.25
Best-fit band-integrated limb-darkening coefficients
c1 0.58115522 0.49271309 0.57505270 0.55407671 0.46913939
c2 0.08416305 0.32553801 0.07068500 0.15318674 0.38874433
c3 -0.1886652 -0.5299505 -0.1131334 -0.2741835 -0.6396845
c4 0.05957531 0.21191308 0.01649965 0.09583456 0.26591443
andphotometriccalibration. Theresultantimagesarein
units of electrons per second. A single exposure for each
of the five targets are shown in Fig. 4, with the target
spectrumoutlined. Thisfigurealreadydemonstratesthe
difference in the individual observation strategies used
for WFC3. The larger subarrays used for HAT-P-1 and
WASP-17 include both the 0th and 1st order spectrum,
while the smaller subarrays only contain the 1st order
spectrum. There is also a clear difference between the
brightest target, HD209458 (V mag = 7.7), and dimmer
targets (V mag >10.3), where a much larger scan area Fig. 5.—Left: exampleexposureinspatialscanmodefromthe
is needed for bright targets so that the detector is not “ima” output. Right: The spectral trace in the cross-dispersion
direction from one pixel column is shown in black, with a series
saturated during the exposure. We also note that each
of aperture sizes marked out by the colored brackets. The best
HAT-P-1exposureincludesthespectraltraceofthecom-
aperture size is determined by minimising the standard deviation
panion star to the target exoplanet host star. ofthelightcurveresidualswhenfitwithastandardMandel&Agol
For each of our five exoplanetary transit datasets we (2002)transitmodelfromthebaseplanetaryparameters.
extract each exposure spectrum with our custom IDL
routine spextract, which optimises the aperture over ther constraints cannot be placed with these individual
whichthetargetspectrumisexposedoneachimage(see datasets.
Fig. 5). Wethencomputetheband-integratedlightcurve For each of the spectroscopic transit observations we
bysummingthefluxoverallexposedwavelengthstoob- computetheband-integratedlightcurveandtransmission
tain the planet-to-star radius ratio (R /R ) and centre spectrum from 1.1–1.7µm. Here we present the results
p ∗
oftransittimebyfittingaMandel&Agol(2002)transit fromeachofthehotJupiterexoplanettransitlightcurves
model, created with non-linear limb-darkening parame- together,discussingthedatawithrespecttotheanalysis
ters, to our data using the IDL code MPFIT. Due to technique and individual observation strategies.
the limited phase coverage of the HST transit observa-
tions,whenfittingtheband-integratedandspectroscopic 4.1. Band-integrated lightcurve marginalisation
lightcurvesforR /R wefixthesystemparameterssuch
p ∗ Weuseoursystematicgridtodeterminetheimpacting
as inclination and a/R from a joint fit between HST
∗ systematicsforeachtransitdataset. Thegridconsistsof
STIS, WFC3, and Spitzer IRAC using a Markov Chain
52systematicmodelsusingacombinationofsystematics
MonteCarlo(MCMC)analysis(Singetal.2016),asfur-
in the polynomial model and the two different exponen-
10 Wakeford et al.
Fig. 6.— The evidence based on the AIC approximation for each systematic model applied to the band-integrated lightcurve for each
exoplanet transit. Top to Bottom: HAT-P-1 (blue), WASP-31 (orange), XO-1 (red), HD209458 (purple), WASP-17 (green). The best-
fittingsystematicmodelfortheband-integratedlightcurveisfilledinwithablacksquareforeachplanetineachplot. Anumberofpoint
ontheHD209458plothavebeenartificiallyshiftedonthescaleshownwiththeiroriginalvalueslistedabovetheshiftedpointsastheyare
stronglydisfavouredbythedata.
tial models (see Table 2). We then perform a marginal- troduction of an additional visit-long linear correction
isation on the computed planet-to-star radius ratio, and makes these models most favourable. The systematic
centreoftransittime,anddeterminethebest-fittingsys- model correction for XO-1 favours higher order polyno-
tematicmodelwiththemaximumweightfromtheband- mials across HST phase and wavelength shift with little
integrated lightcurve. difference incurred with or without the visit-long slope.
For each of our datasets we use marginalisation as a HD209458 strongly disfavours corrections with only the
tool to compute robust transit depths from the band- “HST breathing” effect, with the evidence based on the
integrated lightcurve. This also allows us to determine AIC approximation several hundred points below the
the main systematics impacting the lightcurves from the other systematic correction fits (these points have been
different spectral targets in a common and directly com- artificially shifted on the scale shown with their origi-
parative manner. Figure 6 shows the evidence based nal values listed above the shown points). This makes
on the AIC parameter for each systematic model used the other systematic models relatively stable in that no
to correct the band-integrated lightcurve following the model appears to be greatly favoured above any others.
model numbers in Table 2. We highlight the model For WASP-17b it can clearly be seen in the raw
favoured by this criterion corresponding to the system- lightcurve that there is a strong hook feature present
atic corrections with the highest weighting. From each (Fig. 1), with an orbit-to-orbit repeating pattern in
of the fitting statistics it becomes clearer to see the dif- the residuals. The grid of systematic models used to
ferences across the datasets. For HAT-P-1 there is lit- correct the white lightcurve do not accurately account
tle difference between a large number of the system- for this hook. To effectively correct for this systematic
atic models with distinct groups of models which cor- in the white lightcurve we use the divide-oot routine
rect for only wavelength shifts and no “HST breathing” discussed in section 2.3, which removes common- mode
effectsbeingdisfavouredinthefinalmarginalisation. In- time-dependentsystematics. Wethenadditionallyapply
terestingly, the band-integrated lightcurve of WASP-31 amarginalisationoverthesystematicgridtoremoveany
does not favour systematic models which only correct additional time- independent systematics and determine
for the known “HST breathing” effect, while the in- theplanet-to-starradiusratio,wheremostmodelsincor-
