WASP-8b: CHARACTERIZATION OF A COOL AND ECCENTRIC EXOPLANET WITH SPITZER

When transiting exoplanets pass behind their host stars (a secondary eclipse), the observed flux drop provides a direct measurement of the planet's thermal emission and reflected light. Today, secondary-eclipse observations exist for nearly 30 exoplanets. The Spitzer Space Telescope (Werner et al. 2004) made most of these observations, capturing broadband photometric light curves in six near- and mid-infrared bands (3–24 μm). Each band probes a specific altitude range in a planet's atmosphere. With Bayesian fitting of model spectra, one can quantitatively constrain the atmospheric chemical composition and thermal profile of the planet's photosphere (Madhusudhan & Seager 2010). WASP-8b, with a time-averaged equilibrium temperature of 948 ± 22 K (T_eq, temperature at which blackbody emission balances absorbed energy, assuming zero albedo and efficient heat redistribution), is one of the coolest Jupiter-sized planets yet observed in eclipse, and thus serves as an end member to the set of measured hot-Jupiter atmospheres.

To classify the hot-Jupiter population, Fortney et al. (2008) proposed a separation between moderately and strongly irradiated planets. The higher atmospheric temperatures of the more strongly irradiated planets allow the presence of highly opaque molecules (like TiO and VO) at high altitudes. These strong absorbers produce hot stratospheres (thermal inversion layers). In contrast, for the moderately irradiated hot Jupiters, these absorbers condense and rain out to altitudes below the photosphere, thus presenting no thermal inversions.

In general, the observations agree with this hypothesis, but exceptions indicate that the picture is not yet completely understood. For example, secondary-eclipse observations of the highly irradiated WASP-12b (Madhusudhan et al. 2011; Crossfield et al. 2012), WASP-14b (Blecic et al. 2013), and TrES-3 (Fressin et al. 2010) do not show evidence of thermal inversions. Conversely, XO-1 has an inversion layer even though it receives a much lower stellar irradiation (Machalek et al. 2008). Photochemistry provides one explanation. The non-equilibrium atmospheric chemistry models of Zahnle et al. (2009) suggested that heating from sulfur compounds in the upper atmospheres of hot Jupiters could explain these inversions. Alternatively, Knutson et al. (2010) suggest that strong UV radiation from active stars destroys the high-altitude absorbers.

The Wide-Angle Search for Planets (WASP) Consortium discovered WASP-8b in 2008 (Queloz et al. 2010). The planet orbits the brighter component (WASP-8A) of a binary stellar system. The angular separation (483) with the secondary (WASP-8B) sets a minimum separation of 440 AU between the stars. WASP-8A is a G6 star, with effective temperature T_eff = 5600 K. Color and photometric analyses indicate that WASP-8B is a colder M star (Queloz et al. 2010). WASP-8b is a 2.18 Jupiter-mass (M_Jup) planet with 1.08 times Jupiter's radius (R_Jup) in a retrograde 8.16 day orbit. Its large eccentricity (e = 0.31) should make the planet's dayside temperature vary by hundreds of degrees along the orbit, possibly forcing an unusual climate.

The age of the host star (4 Gyr) is shorter than the planet's orbital circularization time (∼30 Gyr; see, e.g., Goldreich & Soter 1966; Bodenheimer et al. 2001); accordingly, WASP-8b has one of the most eccentric orbits among the ∼10-day-period exoplanets (Pont et al. 2011). The Kozai mechanism (Wu & Murray 2003) may explain the combination of high eccentricity and retrograde orbit orientation. The radial-velocity (RV) drift and the large eccentricity may also indicate a second planetary companion (Queloz et al. 2010).

We obtained six secondary-eclipse light curves of WASP-8b from four visits of the Spitzer Space Telescope, observing in the 3.6, 4.5, and 8.0 μm bands of the Infrared Array Camera (IRAC; Fazio et al. 2004). The eclipse depths determine the planet's dayside infrared emission. Our Markov-chain Monte Carlo (MCMC)-driven radiative-transfer code constrained the molecular abundances and temperature profile of WASP-8b's dayside atmosphere, testing for the expected absence of a thermal inversion and estimating the energy redistribution over its surface. We constrained the orbit of WASP-8b by determining the eclipse epochs and durations. We also modeled the thermal variations along the orbit of the planet to explore the effects of the eccentricity.

Section 2 presents the Spitzer observations of the WASP-8 system. Section 3 describes the photometric and modeling analysis of our secondary eclipse observations. Section 4 gives the orbital dynamical analysis. Section 5 presents our constraints on WASP-8b's atmospheric composition derived from the photometry. Section 6 discusses the effects of eccentricity on the orbital thermal variation of WASP-8b. Finally, Section 7 states our conclusions.

The Spitzer Space Telescope visited WASP-8 four times. From two consecutive eclipse observations, we obtained simultaneous light curves at 4.5 and 8.0 μm. Later, from two more consecutive eclipse observations during the Warm Spitzer mission, we obtained one light curve at 3.6 μm and one at 4.5 μm (see Table 1). The Spitzer pipeline (version 18.18.0) processed the raw data, producing Basic Calibrated Data (BCD).

During the initial minutes of our observations, the telescope pointing drifted ∼0.25 pixels before stabilizing. Throughout the observations, the pointing also jittered from frame to frame (∼0.01 pixel) and oscillated in an hour-long periodic movement (∼0.1 pixel amplitude).

The separation between the centers of WASP-8A and WASP-8B in the IRAC detectors is only 3.7 pixels. Consequently, the signal from the stars overlapped, demanding special care during the data analysis (see Figure 1). Table 2 shows the average and standard deviation of the flux ratio, separation, and position angle (PA) of the secondary star with respect to WASP-8A (derived from our centering routine, see Section 3.2). Our PA values agree with those of Queloz et al. (2010), but our separation values are consistently lower than theirs (483 ± 001).

Zoom In Zoom Out Reset image size

Our Photometry for Orbits, Eclipses, and Transits (POET) pipeline produces light curves from BCD images. Briefly, POET creates a bad pixel mask for each image, finds the center position of the target, executes interpolated aperture photometry, and fits a light-curve model that includes physical and systematic parameters.

3.1. POET: Initial Reduction

3.2. Centering

3.3. Photometry

3.4. Light-curve Modeling

The eclipse depths of WASP-8b are on the order of 0.1% of the system's flux, well below Spitzer's photometric stability criteria (Fazio et al. 2004). Thus, the eclipse light-curve modeling requires a thorough characterization of the detector systematics. Systematic effects have been largely observed and documented; they can have both temporal and spatial components, and vary in strength and behavior for each data set.

The main systematic at 3.6 and 4.5 μm is intrapixel sensitivity variation, M(x, y), where the measured flux depends on the precise position of the target on the array (Stevenson et al. 2012; Charbonneau et al. 2005). In addition, the detectors show a time-dependent sensitivity variation called the ramp effect, R(t), suspected to be caused by charge trapping (Agol et al. 2010) at 8.0 μm, but there are also reports of a ramp in the 3.6 and 4.5 μm bands (e.g., Campo et al. 2011; Nymeyer et al. 2011; Knutson et al. 2011; Deming et al. 2011; Blecic et al. 2013; Stevenson et al. 2010, 2012). The eclipse and both systematic variations entangle to produce the observed light curve. To account for their contributions, we modeled the light curves as

$$\begin{equation} F(x,y,t) = F_{s} M(x,y) R(t) E(t), \end{equation} \tag{ 1 }$$

where F_s is the out-of-eclipse system flux. We used the eclipse model, E(t), from Mandel & Agol (2002). The eclipse is parameterized by the eclipse depth, the mid-point phase, the duration, and the ingress and egress times. For the ingress/egress times we adopted a value of 18.8 minutes, derived from the orbital parameters of the planet. We used this value in all of our eclipse-model fits.

The strength and behavior of the ramp variations are specific to each data set. Many formulae have been applied in the literature (e.g., Deming et al. 2007; Harrington et al. 2007; Knutson et al. 2011; Stevenson et al. 2012). The models are formed with combinations of exponential, logarithmic, and polynomial functions. We tested dozens of equations; the best were:

$$\begin{eqnarray} {\rm risingexp}{:}& \quad R(t) = 1 - e^{-r_{0}(t-t_{0})}\\ {\rm logramp}_{q}{:}& \quad R(t) = 1 + r_{q} [\ln (t-t_{0})]^{q}\\ {\rm linramp}{:}& \quad R(t) = 1 + r_{1}(t-t_{c})\\ {\rm quadramp}_{1}{:}& \quad R(t) = 1 + r_{1}(t-t_{c}) + r_{2}(t-t_{c})^{2}\\ {\rm quadramp}_{2}{:}& \quad R(t) = 1 + r_{2}(t-t_{c})^{2}\\ {\rm loglinear}{:}& \quad R(t) = 1 + r_{1}(t-t_{c}) + r_{4} \ln (t-t_{0}) \end{eqnarray} \tag{ 2 }$$

where t_c is a constant value at the approximated mid-point phase of the eclipse (t_c = 0.515 for this planet). Slight changes in t_c do not significantly affect the fitted eclipse parameters.

We used our Bi-Linearly Interpolated Subpixel Sensitivity (BLISS) mapping technique (Stevenson et al. 2012) to calculate M(x, y). The BLISS method has been found to return a better result than a polynomial fit (Stevenson et al. 2012; Blecic et al. 2013).

To determine the best-fitting parameters of our model (Equation (1)), we used a χ² minimizer with the Levenberg–Marquardt algorithm. We used Bayesian posterior sampling via an MCMC algorithm to explore the phase space and estimate the uncertainties of the free parameters of the light-curve models. Our code implements the Metropolis random walk, which proposes parameter sets from a multivariate normal distribution centered at the current position in the chain, computes χ², and accepts (or rejects) the new set with greater probability for a lower (higher) χ². By generating millions of parameters sets, the algorithm samples the posterior distribution of the model parameters. As a necessary condition for chain convergence, we require the Gelman–Rubin statistic (Gelman & Rubin 1992) to be within 1% of unity for each free parameter between four MCMC chains.

The photometry routine uses the BCD uncertainty images to estimate the uncertainties of the light-curve data points, σ_i. However, since the Spitzer pipeline in general overestimates these uncertainties (it is designed for absolute photometry), we multiply by a constant factor (σ_i → f · σ_i), such that the reduced χ² = 1 in the light-curve fit. This is equivalent to estimating a single σ from the scatter of model residuals. Both methods account for red noise, but ours retains the (usually small) σ_i variations due to aberrant frames.

To determine the best raw light curve (i.e., by selection of photometry method and aperture radius), we calculated the SDNR of the light-curve fit (Stevenson et al. 2012; Campo et al. 2011). Poor fits or data with high dispersion increase SDNR; the optimum data set minimizes the SDNR value. Once we chose the best light curve, we compared the different ramp models according to the Bayesian Information Criterion (Liddle 2007),

$$\begin{equation} {\rm BIC} = \chi ^{2} + k\ln N, \end{equation} \tag{ 8 }$$

where k is the number of free parameters and N the number of data points. The best model minimizes the BIC. The probability ratio favoring one model over a second one is exp (− ΔBIC/2).

3.4.1. wa008bs11 Analysis

This observation started 2.9 hr before the eclipse's first contact. The telescope observed the target in sub-array mode, allowing a high cadence (Table 1). We discarded the initial 15 minutes of observation while the telescope pointing settled. Our data present both intrapixel and weak ramp systematics.

Table 2. WASP-8 System

Event	Flux Ratio		Separation (pixels)		Separation ('')		Position Angle (deg)
	average	stddev	average	stddev	average	stddev	average	stddev
wa008bs11	0.1420	0.0030	3.760	0.013	4.610	0.016	171.32	0.28
wa008bs21	0.1512	0.0017	3.734	0.007	4.541	0.009	171.00	0.15
wa008bs22	0.1600	0.0022	3.737	0.009	4.544	0.011	170.75	0.16
wa008bs23	0.1648	0.0039	3.726	0.017	4.497	0.020	170.78	0.29
wa008bs41	0.1718	0.0020	3.690	0.009	4.513	0.011	170.84	0.17
wa008bs42	0.1794	0.0023	3.686	0.010	4.506	0.012	170.96	0.16

Download table as:  ASCIITypeset image

The 2.25 pixel A aperture with B-subtract photometry minimized SDNR. Table 3 shows the five best-fitting models to the best wa008bs11 light curve. ΔBIC is with respect to the lowest BIC value. The quadramp₁ model is 2.4 times more probable than, and consistent with, the second-best model. The linear (11σ) and quadratic (4σ) terms of the quadramp₁ model (see Table 9) confirm the need for a ramp model. As a general remark, we noted that all the logramp_q models produce similar BIC and eclipse parameter values; therefore, we will refer only to the logramp₁ model in the future. Models with more free parameters do not improve BIC. Following Stevenson et al. (2012), we vary the bin size and the minimum number of data points per bin (mnp) of the BLISS map to minimize the dispersion of the residuals. We required at least 4 points per bin for any data set. The PSF-fitting position precision of 0.01 pixels sets our lower limit for the bin size. For wa008bs11, mnp = 5 and a bin size of 0.015 pixels optimized the fit.

Table 3. wa008bs11 Ramp Model Fits

R(t)	SDNR	ΔBIC	Ecl. Depth (%)
quadramp1	0.0061141	0.00	0.119
risingexp	0.0061148	1.73	0.106
logramp1	0.0061153	2.96	0.096
linramp	0.0061201	6.34	0.063
loglinear	0.0061141	11.10	0.119

Download table as:  ASCIITypeset image

Figure 2 shows the raw, binned and systematics-corrected wa008bs11 light curves with their best-fitting model. We considered the correlated noise in the residuals as well (Pont et al. 2006). Figure 3 shows the rms of the residuals versus bin size. The wa008bs11 rms curve deviates above the expected rms for pure Gaussian noise. Following Winn et al. (2008), to account for the correlated noise, we weighted the light-curve uncertainties by the factor β (the fractional rms excess above the pure Gaussian rms at the bin size corresponding to the eclipse duration). For wa008bs11 we found β = 2.4. We inspected all the pairwise correlation plots and histograms and found only unimodal Gaussian distributions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Alternatively, the residual-permutation (also known as prayer bead) algorithm is sometimes used to assess correlated noise in a fit. In this method, we cyclically shift the residuals from the best model by one frame, add them back to that model, and refit, repeating until we shift the residuals back to their original positions. This generates a distribution of values for each parameter, from which we estimate the parameter uncertainties. The eclipse-depth uncertainty is 0.021%, similar to the value found with the Winn et al. (2008) method (see Table 9). However, we are cautious. Although prayer bead has been broadly used for the analysis of exoplanet light curve and RV fits (e.g., Southworth 2008; Bouchy et al. 2005; Pont et al. 2005; Gillon et al. 2007; Knutson et al. 2008; Cowan et al. 2012), we have found no detailed description of its statistical properties in the literature.

3.4.2. wa008bs23 Analysis

With the same observing setup as wa008bs11, this observation started 3.3 hr prior to the eclipse's first contact. This data set also presented both intrapixel and ramp variations. Note that the intrapixel systematic is weaker than at 3.6 μm, attributed to the smaller degree of undersampling at larger wavelengths by Charbonneau et al. (2005) and Morales-Calderón et al. (2006). Even though the pointing stabilized only after the initial 20 minutes, the light curve did not deviate significantly; therefore, we included all data points in the analysis. We noted two sudden pointing and PA deviations near phase 0.519. After each incident, the telescope resumed its position within ∼10 s (Figure 4). Micrometeorite impacts on the telescope can explain the abrupt deviations. Simultaneously, we measured a slight increase in the background sky flux dispersion, which returned to normality shortly after. The target flux did not show any extraordinary fluctuations during these incidents. However, the points outside the normal pointing range were eliminated by the BLISS map's mnp criterion.

Zoom In Zoom Out Reset image size

The analysis is analogous to wa008bs11. The SDNR indicated clearly that the 2.25 pixel A aperture with B-Subtract photometry produced the lowest dispersion. The best-fitting ramp model is logramp₁, which is 21 times more probable than the rising exponential ramp (Table 4). The BLISS map is optimized at mnp = 4 and a bin size of 0.025 pixels.

Table 4. wa008bs23 Ramp Model Fits

R(t)	SDNR	ΔBIC	Ecl. Depth (%)
logramp1	0.0073830	0.00	0.0677
risingexp	0.0073832	6.10	0.0730
quadramp1	0.0073833	9.40	0.0777
loglinear	0.0073830	10.84	0.0685
linramp	0.0073846	12.19	0.0564

Download table as:  ASCIITypeset image

An initial MCMC run showed a significant linear correlation between the system flux and the r₁ parameter of the logramp, which prevented the MCMC chain from converging. We solved this problem by transforming the correlated parameters into an orthogonal set of parameters, rerunning the MCMC chain, and inverting the transformation on the resulting parameter values (Stevenson et al. 2012). Figures 2 and 3 show the wa008bs23 light curves with the best-fitting model and rms of the residuals versus bin size.

3.4.3. wa008bs21 and wa008bs41 Analysis

We simultaneously observed wa008bs21 and wa008bs41 in full-array mode. Prior to the eclipse observation, we exposed the detector (a "preflash" observation, Knutson et al. 2009) for 25 minutes to a bright H ii region, with coordinates α = 20^h21^m3928 and δ = +37°31'036, to minimize the ramp systematic variation. The secondary-eclipse observation started only 26 minutes before the first contact. The telescope pointing stabilized quickly, so fortunately we needed to remove only the initial four minutes of observation. Every 12 s, the detector recorded two consecutive images (two-second exposures) at 4.5 μm and one image at 8.0 μm. (Table 1).

The SDNR analysis of wa008bs21 showed that a 3.5 pixel A aperture with 1.6 pixel B-Mask photometry minimizes the dispersion (Figure 5). The ramp models indicated a negligible ramp variation. Accordingly, a fit without a ramp model yielded the lowest BIC. Table 5 shows the four best-fitting models for the best wa008bs21 data set. The no-ramp model is 15 times more probable than the quadramp₂ model.

Zoom In Zoom Out Reset image size

Table 5. wa008bs21 Ramp Model Fits

R(t)	SDNR	ΔBIC	Ecl. Depth (%)
no ramp	0.0036223	0.00	0.0718
quadramp2	0.0036195	5.76	0.1189
linramp	0.0036223	7.56	0.0714
quadramp1	0.0036197	13.14	0.1170

Download table as:  ASCIITypeset image

Because of the shorter out-of-eclipse observation, the system flux is less constrained for wa008bs21 than for the wa008bs11 or wa008bs23 events. Combined with a correlation between the eclipse depth and system flux (revealed by MCMC), the lower precision of the system flux translates into a larger eclipse depth uncertainty. Nevertheless, the wa008bs21 fit parameters were consistent among the different apertures (Figure 6). The optimum parameters of the BLISS map are mnp = 5 and a bin size of 0.025 pixels.

Zoom In Zoom Out Reset image size

The 8.0 μm detector did not present an intrapixel pattern like the 3.6 or 4.5 μm detectors. However, some of the raw light curves for different apertures and photometry methods showed large scatter and presented strong oscillations, producing implausible fit parameters. A pixelation effect (Anderson et al. 2011; Stevenson et al. 2012) might be responsible. As a consequence, we were unable to fit the eclipse parameters unambiguously for this data set alone. Normally we study the events individually to select the best aperture and photometry method, but in this case we used a joint fit with the best wa008bs21 data set and model to help constrain the 8.0 μm eclipse curve, sharing the eclipse duration and mid-point parameters. The 3.5 pixel A aperture with 1.6 pixel B-mask photometry for wa008bs41 minimized the joint SDNR (Figure 7).

Zoom In Zoom Out Reset image size

Table 6 compares the four best-fitting ramp models for the best wa008bs41 light curve. A linear ramp minimizes BIC, and is 20 times more probable than the next-best model. Figures 2 and 3 show the wa008bs21 and wa008bs41 light curves with their best-fitting models and rms of the residuals versus bin size, respectively.

Table 6. wa008bs41 Ramp Model Fits

R(t)	SDNR 4&2	ΔBIC	Ecl. Depth (%)
linramp	0.0030766	0.00	0.0931
quadramp1	0.0030744	5.94	0.1308
risingexp	0.0030754	6.33	0.1150
logramp1	0.0030758	6.50	0.1078

Download table as:  ASCIITypeset image

3.4.4. wa008bs22 and wa008bs42 Analysis

The observing setup of these events was identical to wa008bs21 and wa008bs41, including the preflash observation. The pointing of this observation drifted noticeably more than in the other observations, moving more than 0.4 pixels during the initial 30 minutes and stabilizing only during the eclipse. As a consequence, the illumination level of the individual pixels changed during the beginning of the eclipse. The ramp variation, which depends on the illumination (Knutson et al. 2008), was disrupted.

The wa008bs22 event, having a negligible ramp variation, was little affected by the telescope pointing shift. The SDNR calculation for wa008bs22 indicated the 1.8 pixel B-Mask photometry with 3.75 pixel A aperture as the best data set. A light-curve model without a ramp (Table 7) is 639 times more probable than the quadramp₂ model. The optimal BLISS map has mnp = 4 and a bin size of 0.02 pixels.

Table 7. wa008bs22 Ramp Model Fits

R(t)	SDNR	ΔBIC	Ecl. Depth (%)
no ramp	0.0025274	0.00	0.0814
quadramp2	0.0025253	12.94	0.1224
logramp1	0.0025267	14.22	0.0921
risingexp	0.0025274	15.04	0.0814

Download table as:  ASCIITypeset image

In contrast, we discarded the initial wa008bs42 light curve past the eclipse ingress due to the disrupted ramp variation. The eclipse model parameters are thus less constrained. By this point, we already had single-channel fits for the rest of the data, so we tuned the wa008bs42 analysis in a joint fit with all the other events, sharing the eclipse duration and mid-time. SDNR indicates the B-subtract method with 4.00 pixel A aperture as the best data set. Table 8 presents the four best-fitting models. The eclipse depth is consistent with the wa008bs41 depth.

Table 8. wa008bs42 Ramp Model Fits

R(t)	SDNR	ΔBIC	Ecl. Depth (%)
linramp	0.0032320	0.00	0.0932
quadramp1	0.0032312	6.19	0.0892
logramp1	0.0032321	6.67	0.0938
risingexp	0.0032330	7.08	0.0961

Download table as:  ASCIITypeset image

3.4.5. Final Joint-fit Analysis

From the three individual fits to the 4.5 μm observations we found eclipse depths of 0.072% ± 0.021%, 0.086% ± 0.022%, and 0.068% ± 0.007% for wa008bs21, wa008bs22, and wa008bs23, respectively. The weighted mean of the depths is 0.0692% ± 0.0065%. With a dispersion of 0.0062% around the mean, this is not larger than the individual uncertainties, thus we found no evidence for temporal variability. This dispersion corresponds to 10% of the mean eclipse depth. The consistency permitted a joint analysis of all observations. We used the best light curves and models found in the individual fits, where all events shared the eclipse duration, the three 4.5 μm events shared their eclipse depth, the two 8.0 μm events shared their eclipse depth, the simultaneous wa008bs21 and wa008bs41 events shared their eclipse mid-point phases, and the wa008bs22 and wa008bs42 events shared their eclipse mid-point phases. Table 9 shows the light-curve modeling setup and results. We used these joint-fit results for the orbital and atmospheric analysis. An electronic supplement contains the best light curves, including centering, photometry, and the joint fit.

Table 9. Best-fit Eclipse Light-curve Parameters

Parameter	wa008bs11	wa008bs21	wa008bs22	wa008bs23	wa008bs41	wa008bs42
Array Position ($\bar{x}$, pix)	14.74	20.76	20.39	14.65	19.11	18.72
Array Position ($\bar{y}$, pix)	15.07	233.30	233.30	15.12	230.27	230.20
Position Consistencya (δx, pix)	0.0072	0.0223	0.0220	0.0097	0.0273	0.0254
Position Consistencya (δy, pix)	0.0118	0.0228	0.0236	0.0101	0.0272	0.0274
A Aperture Size (pix)	2.25	3.5	3.75	2.25	3.5	4.0
WASP-8B Photometric Correction	subtract	1.6 mask	1.8 mask	subtract	1.6 mask	subtract
System Flux Fs (μJy)	144555.0(21.0)	91369.9(8.5)	90850.3(8.5)	87473.0(21.0)	32892.5(6.6)	34949.8(8.9)
Eclipse Depth (%)	0.113(18)	0.0692(68)	0.0692(68)	0.0692(68)	0.093(23)	0.093(23)
Brightness Temperature (K)	1552(85)	1131(35)	1131(35)	1131(35)	938(99)	938(99)
Eclipse Mid-point (orbits)	0.51428(34)	0.51446(37)	0.51468(41)	0.51536(28)	0.51446(37)	0.51468(41)
Eclipse Mid-point (MJDUTC)b	5401.4981(28)	4822.2301(31)	4814.0732(33)	5409.6656(23)	4822.2301(31)	4814.0732(33)
Eclipse Mid-point (MJDTDB)b	5401.4989(28)	4822.2309(31)	4814.0739(33)	5409.6663(23)	4822.2309(31)	4814.0739(33)
Eclipse Duration (t4-1, hr)	2.600(78)	2.600(78)	2.600(78)	2.600(78)	2.600(78)	2.600(78)
Ingress/Egress Time (t2-1, hr)	0.314	0.314	0.314	0.314	0.314	0.314
Ramp Equation (R(t))	quadramp1	None	None	logramp1	linramp	linramp
Ramp, Linear Term (r1)	0.0707(70)	⋅⋅⋅	⋅⋅⋅	0.000504(45)	0.205(22)	0.246(37)
Ramp, Quadratic Term (r2)	−3.17(75)	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
Ramp, Phase Offset (t0)	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	0.4917	⋅⋅⋅	⋅⋅⋅
BLISS Map (M(x, y))	Yes	Yes	Yes	Yes	No	No
Minimum Number of Points Per Bin	5	5	4	4	⋅⋅⋅	⋅⋅⋅
Total Frames	64320	2024	2024	64320	1012	1012
Frames Usedc	62203	1936	1879	64072	966	725
Rejected Frames (%)	3.29	4.35	7.16	0.39	4.54	28.36
Free Parametersd	6	4	4	5	5	5
BIC Value	80444.5	80444.5	80444.5	80444.5	80444.5	80444.5
SDNR	0.0053772	0.0036250	0.0035698	0.0073926	0.0030768	0.0032320
Uncertainty Scaling Factor	0.3075	1.0280	1.0077	1.0902	1.1187	1.1382
β Correction	2.4	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
Photon-limited S/N (%)	37.00	94.71	96.59	89.66	76.94	71.00

Notes. The values quoted in parentheses are the 1σ uncertainties. ^arms frame-to-frame position difference. ^bMJD = BJD − 2,450,000. ^cWe exclude frames during instrument/telescope settling, for insufficient points at a given BLISS knot, and for bad pixels in the photometry aperture. ^dIn the individual fits. Joint fit had 19 free parameters.

Download table as:  ASCIITypeset image

WASP-8b's high eccentricity (e = 0.31) implies that its separation from WASP-8A at periapsis (0.055 AU) is about half that at apoapsis. Given the argument of periapsis (ω = −8586), the secondary eclipse nearly coincides with the periapsis. The planet, therefore, receives over twice as much flux at eclipse as it would if the orbit were circular, explaining in part our high brightness temperature (see Table 9).

Secondary-eclipse times can refine estimates of ecos ω from RV data. The four eclipse events occurred at an average eclipse phase of 0.514695 ± 0.00018. After subtracting a coarse light-time correction of 2a/c = 80 s from this average phase, we calculated ecos ω = 0.02290 ± 0.00028 (see Equation (3) of Charbonneau et al. 2005). This is consistent with Queloz et al. (2010), and photometrically confirms the nonzero eccentricity of the planet's orbit (we fit ecos ω below without relying on the low e approximation).

The eclipse timings were combined with 130 available RV data points and with transit data from Queloz et al. (2010) using the method described by Campo et al. (2011) and Nymeyer et al. (2011). Forty-eight in-transit RV points were removed due to the Rossiter–McLaughlin effect.

Our fit presented a moderate improvement to the orbital parameters of WASP-8b (Table 10), except for the period. While Queloz et al. (2010) used several transits to measure the period, we used their published mid-point epoch (a single date); hence, our period is constrained mostly by our eclipses and the RV data, and thus has a larger uncertainty. By themselves, the secondary eclipses have a period of 8.158774 ± 0.00040 days and a midpoint epoch of BJD 2455409.6629 ± 0.0017 (TDB), not significantly ([5.9 ± 4.3] × 10⁻⁵ days) longer than the period found by Queloz et al. (2010). The transit and eclipse periods place a 9.8 × 10⁻⁵° day⁻¹ (3σ) upper limit on possible apsidal precession, nearly three orders of magnitude larger than the theoretical expectation for tidal effects (Ragozzine & Wolf 2009).

Table 10. Eccentric Orbital Model

Parameter	This Work	Queloz et al. (2010)
esin ω	−0.3078 ± 0.0020	−0.3092 ± 0.0029
ecos ω	0.02219 ± 0.00046	0.023 ± 0.001
e	0.309 ± 0.002	0.310 ± 0.0029
ω (°)	−85.00 ± 0.08	−85.73 ± 0.18
P (days)	8.158719 ± 0.000034	8.158715 ± 0.000016
T0 (MJDTDB)	4679.33486 ± 0.00057	4679.33509 ± 0.00050
K (ms−1)	221.9 ± 0.6	222.23 ± 0.8
γC (ms−1)	−1 565.9 ± 0.6	−1 565.8 ± 0.21
γH (ms−1)	−1 547.4 ± 0.4	−1 548.1 ± 0.6
${\dot{\gamma }}$ (ms−1 yr−1)	58.1 ± 1.2	58.1 ± 1.3
Reduced χ2	4.1	0.86

Download table as:  ASCIITypeset image

We use our IRAC observations of thermal emission from WASP-8b to constrain the thermal structure and composition of the dayside atmosphere of the planet. The Spitzer bandpasses at 3.6, 4.5, and 8.0 μm contain strong spectral features due to several carbon and oxygen-based molecules that are expected in hot-Jupiter atmospheres. Methane (CH₄) has strong spectral features in the 3.6 and 8.0 μm bands, carbon monoxide (CO) and carbon dioxide (CO₂) have features at 4.5 μm, while water vapor (H₂O) has features in all three bands (Madhusudhan & Seager 2010). The spectral features of the various molecules appear as absorption troughs or emission peaks in the emergent spectrum depending on whether the temperature decreases or increases with altitude, respectively. Consequently, strong degeneracies exist between the temperature structure and molecular composition derived from a spectral data set (e.g., Madhusudhan & Seager 2010). Nevertheless, photometric observations made with Spitzer have been successfully used to constrain chemical compositions and temperature structures in many exoplanetary atmospheres (e.g., Barman et al. 2005; Burrows et al. 2007; Knutson et al. 2008; Madhusudhan & Seager 2009; Stevenson et al. 2010; Madhusudhan et al. 2011).

We model the dayside emergent spectrum of WASP-8b using the atmospheric modeling and retrieval method of Madhusudhan & Seager (2009, 2010). The model computes line-by-line radiative transfer in a plane-parallel atmosphere assuming hydrostatic equilibrium, local thermodynamic equilibrium, and global energy balance. We assume a Kurucz model for the stellar spectrum (Castelli & Kurucz 2004) given the stellar parameters. The pressure–temperature (P – T) profile and molecular mixing ratios are free parameters in the model, which can be constrained from the data. The P – T profile comprises six free parameters and the mixing ratio of each molecular species constitutes an additional free parameter. Following Madhusudhan & Seager (2009), we parameterize the mixing ratio of each species as deviations from thermochemical equilibrium assuming solar elemental abundances (Burrows & Sharp 1999). We include the dominant sources of opacity expected in hot-Jupiter atmospheres, namely molecular absorption due to H₂O, CO, CH₄, and CO₂ (Freedman et al. 2008; R. S. Freedman 2009, private communication), and H₂–H₂ collision induced absorption (Borysow 2002). We explore the model parameter space in a Bayesian way using an MCMC sampler (Madhusudhan & Seager 2010, 2011). Given the limited number of observations (N_obs = 3), our goal is not to find a unique model fit to the data; instead, we intend to constrain the region of atmospheric parameter space that is allowed or ruled out by the data.

Our observations rule out a thermal inversion in the dayside atmosphere of WASP-8b. This is evident from the planet–star flux contrasts in the three IRAC bands at 3.6, 4.5, and 8.0 μm. In the presence of a thermal inversion, the planet–star flux contrasts in the 4.5 and 8.0 μm bands are both expected to be greater than the flux contrast in the 3.6 μm band (Burrows et al. 2008; Fortney et al. 2008; Madhusudhan & Seager 2010), due to spectral features of the dominant molecules appearing as emission peaks as opposed to absorption troughs. However, the low 4.5 and 8.0 μm flux contrasts relative to the 3.6 μm contrast require significant absorption due to H₂O and CO across the spectrum, and hence the lack of a thermal inversion in the atmosphere. Figure 8 shows model spectra of WASP-8b with no thermal inversion in the temperature profile. The observed 4.5 and 8.0 μm flux contrasts are explained to a good level of fit by a model without a thermal inversion and with solar abundance composition, as shown by the green curve in Figure 8. Our inference of the lack of a thermal inversion in WASP-8b is independent of any assumption about chemical composition or C/O ratio (e.g., Madhusudhan & Seager 2011). The lack of a thermal inversion in WASP-8b is not surprising, since it is amongst the cooler population of irradiated hot Jupiters, which are not expected to host inversion-causing species such as TiO or VO in their upper atmosphere (Fortney et al. 2008; Spiegel et al. 2009).

Zoom In Zoom Out Reset image size

Our models are unable to reproduce the high planet–star flux contrast observed in the 3.6 μm IRAC band, independent of the composition. The major sources of absorption in the 3.6 μm band are H₂O and CH₄. In principle, decreasing the CH₄ and/or H₂O abundances can lead to a higher 3.6 μm contrast. However, as shown by the red curve in Figure 8, such an increase also simultaneously increases the contrast in the 8.0 μm band, thereby worsening the fit overall. Another hindrance to fitting the observed 3.6 μm contrast is that it requires a hotter P – T profile, with T ≳ 1550 K in the lower atmosphere, predicts much higher fluxes in the 4.5 and 8.0 μm bands than observed. On the other hand, a cooler P – T profile than shown in Figure 8 would provide a better fit in the 4.5 and 8.0 μm bands, but would further worsen the fit in the 3.6 μm band. Consequently, we choose an intermediate P – T profile that provides a compromise fit to all three data points.

Although the one-dimensional (1D) models shown in Figure 8 output less energy than the instantaneous incident irradiation during the eclipse (concurrent with periastron passage), they output ∼20% higher energy compared to the time-averaged incident irradiation received at the substellar point. Considering that a pseudo-synchronous rotation should facilitate the redistribution of energy to the night side, the high emission measured suggests that WASP-8b is quickly reradiating the incident irradiation on its dayside hemisphere, i.e., nearly zero day–night redistribution. Such a scenario would lead to a large day–night temperature contrast in the planet which can be confirmed by thermal phase curves of the planet observed using warm Spitzer (e.g., Knutson et al. 2009). The high emergent flux also implies a very low albedo, as with most hot-Jupiter planets (Cowan & Agol 2011b).

As seen in the previous section, the 3.6 μm brightness temperature is anomalously higher than expected. The hemisphere-averaged equilibrium temperature for instantaneous reradiation (time-averaged around the orbit) is only 948 K; even the instantaneous equilibrium temperature at periapsis, 1128 K, is far lower than this observation. Thus, we modeled the orbital thermal variation due to the eccentricity to determine if such a high temperature is possible from irradiation alone.

Following Cowan & Agol (2011a), we solved the energy balance equation in a one-layer latitude–longitude grid over the planetary surface. The change in temperature of a cell with time, dT/dt, is determined by the difference between the absorbed flux from the star and the re-emitted blackbody flux,

$$\begin{equation} \frac{d T}{d t} = \frac{1}{c_{h}} \left[ (1-A) \sigma T^{4}_{\rm eff} \left(\frac{R_{*}}{r(t)} \right)^{2}\cos \psi (t)\ -\ \sigma T^{4} \right], \end{equation} \tag{ 9 }$$

where c_h is the heat capacity per unit area; T_eff and R_* are the star's effective temperature and radius, respectively; r(t) is the planet–star separation; cos ψ(t) = sin λ max(cos Φ(t), 0) is the cosine of the angle between the vectors normal to the planet surface and the incident radiation, with λ the latitude of the cell and Φ(t) the longitude from the substellar meridian; σ is the Stefan–Boltzmann constant.

Tidal interactions drive the planet's rotational angular velocity (ω_rot) toward synchronization with the orbital angular velocity (ω_orb). Hence, if the spin synchronization timescale (e.g., Seager & Hui 2002; Goldreich & Soter 1966) is shorter than the system age, we expect ω_rot = ω_orb. In the case of WASP-8b, the timescale for tidal synchronization is on the order of 0.05 Gyr, much shorter than the age of the star. However, a planet in an eccentric orbit, where ω_orb changes in time, is actually expected to reach a pseudo-synchronization state (e.g., Langton & Laughlin 2008; Hut 1981), in which the planet does not exchange net angular momentum with its orbit. The planet acquires then a constant rotational angular velocity close to the orbital angular velocity at periastron (ω_{orb, p}). In the literature we found different predictions for this equilibrium angular velocity, from 0.8 ω_{orb, p} (Hut 1981) to 1.55 ω_{orb, p} (Ivanov & Papaloizou 2007).

The tidal evolution drives the orbit of a planet toward zero obliquity in a timescale similar to the spin synchronization (Peale 1999). We thus adopted zero obliquity for our simulations. We also assumed A = 0, supported by the atmospheric analysis (Section 5). Beyond these assumptions, the parameters of interest that control Equation (9) are the radiative time $\tau _{\rm rad}=c_{h}/\sigma T^{3}_{0}$ (where T₀ is the substellar equilibrium temperature at periastron) and the rotational angular velocity of the planet ω_rot (which determines the substellar longitude of a cell through the equation dΦ(t)/dt = ω_rot − ω_orb(t)). With these definitions Equation (9) can be rewritten as

$$\begin{equation} \frac{d T}{d t} = \frac{T_{0}}{\tau _{\rm rad}} \left[ \left(\frac{a(1-e)}{r(t)} \right)^{2} {\rm max}(\cos \Phi (t),0)\ -\ \left(\frac{T}{T_{0}}\right)^{4} \right]. \end{equation} \tag{ 10 }$$

We derived the temperature of each cell as a function of time to study its thermal evolution. Assuming that each cell emits as a blackbody, we calculated the photometric phase curve of the planet by integrating over the hemisphere observable from Earth, weighted by the viewing geometry. Our simulations were for planets nearly in pseudo-synchronous rotation (ω_rot = 0.8, 1.0, and 1.5 ω_{orb, p}). We tested values of τ_rad between 1 and 10³ hr.

Figure 9 shows simulated brightness-temperature light curves of WASP-8b after reaching a periodic stationary state (after a few τ_rad). We noted that the higher irradiation at periastron is not the only contribution to a higher temperature. For ω_rot ⩾ ω_{orb, p}, the substellar angular velocity (dΦ/dt) is minimum during periastron, allowing the temperature to increase due to the longer exposure to the irradiation. For ω_rot < ω_{orb, p}, the substellar angular velocity is negative for an instant around periastron. Later, when the planet emerges from secondary eclipse the over-heated region becomes observable from Earth. As a result, the light curve shows a delayed maximum.

Zoom In Zoom Out Reset image size

Our models show that for large radiative timescales, the temperatures at secondary eclipse are lower than 1150 K, regardless of ω_rot. For radiative times shorter than ∼10² hr, the temperatures can be as high as 1400 K, similar to the 3.6 μm measurement (Figure 9, top panel). However, these models still cannot explain the observed brightness-temperature discrepancy with wavelength.

The study of eccentric hot-Jupiter atmospheric circulation by Kataria et al. (2012) hints at a resolution to this discrepancy. Their Figure 4 (top panel) shows that, as the planet passes through periapsis, the time that the peak temperature is reached varies as a function of pressure. This is typical of their simulations (T. Kataria 2012, private communication). If this differential response is significant in WASP-8b, it would introduce a discrepancy in the observations since the Spitzer bands sample different altitudes (see Figure 9, right panel).

Another possibility is to compare the radiative and advective timescales at the altitudes sampled by each band. Evaluating Equation (1) of Fortney et al. (2008) using WASP-8b's P – T profile indicates that τ_rad increases with depth between 0.1 and 1.0 bar, so there should be less longitudinal temperature contrast at depth. On the other hand, models of Kataria et al. (2012) show that wind speeds decrease with depth, and thus τ_adv also increases with depth. If the increase of τ_adv with depth is sharper than that of τ_rad, then one would expect less-homogenized temperatures at depth (but still above the photosphere). Hence, the rise in temperature (due to the increasing incident irradiation) near periapsis could be more pronounced at 3.6 μm than at longer wavelengths, given the weighting functions of Figure 8.

Spitzer observed secondary eclipses of WASP-8b in the 3.6, 4.5, and 8.0 μm IRAC wavebands. In our joint-fit model, we estimate eclipse depths of 0.113%  ±  0.018%, 0.069%  ±  0.007%, and 0.093%  ±  0.023% at 3.6, 4.5, and 8.0 μm, respectively. These depths correspond to brightness temperatures of 1552, 1131, and 938 K, respectively. Although the 3.6 μm eclipse depth is unexpectedly large, most of the ramp models had consistent depths (within 1σ), while those with inconsistent depths fit the data poorly.

Considering the P – T profile of WASP-8b, KCl, ZnS, Li₂, LiF, or Na₂S clouds could form (see Figure 2(a) of Lodders & Fegley 2006). In analogy to brown dwarfs, partial cloud coverage can cause photometric variability (Artigau et al. 2009); however, our three 4.5 μm observations, spanning 1.5 years, have consistent eclipse depths, suggesting no temporal variation at secondary eclipse above a hemispheric-mean level of ∼35 K (1σ). A moderate cloud layer at altitudes higher than those probed by Spitzer would produce a featureless planetary spectrum at wavelengths shorter than 2 μm (Pont et al. 2008; Miller-Ricci Kempton et al. 2012) and would block some of the stellar flux, decreasing the temperatures at levels probed by Spitzer. Yet, the observed temperatures, which exceed the time-averaged equilibrium temperature, challenge this idea.

Given the high eccentricity, spin-orbit misalignment, and observed RV drift in the WASP-8 system, Queloz et al. (2010) suggested the existence of an additional, unseen body in the system. Our orbital analysis is consistent with theirs. It also improves the orbital parameters and extends the baseline of sampled epochs. This constrains the long-term evolution of the orbit and aids the search for a second planet, for example through the study of timing variations (e.g., Agol et al. 2005).

The eclipse depths probe the dayside atmosphere of WASP-8b. Our results rule out the presence of a thermal inversion layer, as expected, given the irradiation level from the host star. A model with solar-abundance composition explains the 4.5 and 8.0 μm planet–star flux contrast; however, including the high 3.6 μm flux contrast requires models that output nearly 20% of the orbit-averaged incident irradiation, independent of the atmospheric composition. If the orbit were circular (and thus the irradiation steady-state), the high brightness temperatures would indicate a very low energy redistribution to the night side of the planet. For an eccentric planet, it at least indicates a short τ_rad (Figure 9).

By modeling the orbital thermal variations due to the eccentricity of the orbit, we determined that it is possible for WASP-8b to achieve temperatures as high as the 3.6 μm brightness temperature. However, the differing brightness temperatures in the other two bands remain puzzling. Neither the radiative-transfer model (Section 5) nor the phase-variation model (Section 6) embraces all the physics of the problem. The radiative transfer code is a 1D, steady-state model representing typical dayside conditions. The phase-variation model describes emission as a blackbody on a single-layer grid; it does not consider absorption or emission features from the species in the atmosphere. Clouds (e.g., Cushing et al. 2008), atmospheric dynamics (e.g., Showman et al. 2009), and photochemistry (e.g., Moses et al. 2011) are not directly considered by these models.

What we can say for certain is that the assumptions of our simple models have been violated, which is not surprising for this eccentric planet. While it may be possible to construct consistent, realistic models, model uniqueness may be elusive until more and better data are available.

Relatively few exoplanets with equilibrium temperatures below 1500 K have been observed at secondary eclipse (Cowan & Agol 2011b). The same is true for eccentric planets. The characterization of WASP-8b in this work thus addresses a particularly interesting, and challenging, region of the exoplanet phase space. Observation of other planets with similar equilibrium temperatures or eccentricities will help discover the physics that drive these unusual atmospheres.

We thank Ian Crossfield for his help with the PSF fitting algorithm. We thank contributors to SciPy, Matplotlib, and the Python Programming Language; the free and open-source community; the NASA Astrophysics Data System; and the JPL Solar System Dynamics group for software and services. P.C. is supported by the Fulbright Program for Foreign Students. N.M. acknowledges support from the Yale Center for Astronomy and Astrophysics through the YCAA postdoctoral fellowship. This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Partial support for this work was provided by NASA through awards issued by JPL/Caltech and by the Astrophysics Data Analysis Program, grant NNX13AF38G.