THERMAL PHASE VARIATIONS OF WASP-12b: DEFYING PREDICTIONS

Occultations. Transits are modeled using the IDL implementation of Mandel & Agol (2002) with fixed nonlinear limb darkening. To determine the limb darkening coefficients, we fit a four-parameter Claret (2000) model to a Kurucz stellar model with [Fe/H] = 0.3, T_eff = 6250 K, and log g = 4.5 (Kurucz 1979, 2005). We set the eccentricity to zero and fix the orbital period to the value from Maciejewski et al. (2011). The impact parameter, b, and the geometrical factor, a/R_*, are allowed to vary freely. Eclipses are modeled using the uniform-disk version of the Mandel & Agol (2002) expressions, and both eclipses are modeled simultaneously with the same parameters. This means that we do not look for eclipse depth variability (searches for such variability have so far only resulted in upper limits: Agol et al. 2010; Knutson et al. 2011). We account for light-travel time within the system, but this is only a matter of 22.5 s and does not affect our analysis. By the same token, we neglect eclipse mapping effects for the planet (Williams et al. 2006; Rauscher et al. 2007; Agol et al. 2010), since we are insensitive to the resulting offset in-eclipse time of less than a minute, let alone the detailed ingress/egress morphology.

Diurnal phases. The planet's temperature and hence brightness vary as a function of local stellar time. This inhomogeneous intensity is modeled with three parameters: the orbit-averaged planet/star flux ratio, 〈F_p/F_*〉, the semi-amplitude of thermal phase variations, A_therm, and the offset of the phase peak from superior conjunction, α_max (α_max < 0 corresponds to a peak prior to superior conjunction and therefore to an eastward-advected hot spot and super-rotating winds). The phase variations have a sinusoidal shape, corresponding to a sinusoidal longitudinal brightness profile for the planet (Cowan & Agol 2008).¹² Note that in the limit of poor recirculation, the longitudinal temperature profile should be more akin to a half-sine (uniformly dark on the night side), leading to thermal phase variations similar to the Lambert phase function: broader minimum, briefer maximum. We neglect reflected starlight, which does not contribute appreciably at these wavelengths.

Ellipsoidal variations. Because of the planet's small semimajor axis and inflated radius, it is likely that it—and possibly its host star—is ellipsoidal in shape rather than spherical. This leads to changes in cross-sectional area throughout the planet's orbit. To a good approximation, these variations are sinusoidal with a period half of the orbital period; the maxima occur at quadrature, when we are seeing the long axes of the star and planet, and minima at superior and inferior conjunction, when we are seeing the short axes of the two bodies.¹³

Li et al. (2010), Leconte et al. (2011), and Budaj (2011) all predict that the projected area of WASP-12b should vary by approximately 10% (peak to trough) due to its prolate shape, whether it is modeled as a prolate ellipsoid or a partially filled Roche lobe. Given the mid-infrared planet/star flux ratio of F_p/F_* ≈ 4 × 10⁻³ (Campo et al. 2011), we expect to see ellipsoidal variations in the planet at the level of ΔF/F_* ≈ 4 × 10⁻⁴ (or a semi-amplitude of 2 × 10⁻⁴).

The presence of a massive companion should also produce ellipsoidal variations in the star, as seen at optical wavelengths in the system HAT-P-7 (Welsh et al. 2010). Using the expressions given in Faigler & Mazeh (2011), we estimate the semi-amplitude of these variations to be ∼4 × 10⁻⁵ at all wavelengths. We therefore expect that at thermal wavelengths, the ellipsoidal variations of the system should be dominated by the shape of the planet and not that of its host star.

We experimented with different functional forms for the planet's phase variations in an effort to reduce correlations between astrophysical variables. Our best model in this regard is

$$\begin{equation} \frac{F_{p}}{F_{*}} = \left\langle \frac{F_{p}}{F_{*}}\right\rangle + A_{\rm therm}\cos (\alpha -\alpha _{\rm max}) - A_{\rm ellips}\cos (2\alpha), \end{equation} \tag{ 1 }$$

where A_therm and A_ellips are the semi-amplitudes of diurnal and ellipsoidal phase variations, respectively, and α is the phase angle (α = 0 at superior conjunction, α = π at inferior conjunction).

The secondary eclipse depth is related to the model variables by

$$\begin{equation} \frac{F_{\rm day}}{F_{*}} = \left\langle \frac{F_{p}}{F_{*}}\right\rangle + A_{\rm therm}\cos \alpha _{\rm max} - A_{\rm ellips}, \end{equation} \tag{ 2 }$$

and to first order the associated uncertainties can be propagated as

$$\begin{eqnarray} \sigma _{F_{\rm day}/F_{*}}^{2} = & \sigma _{\left\langle F_{p}/F_{*}\right\rangle }^{2} + \cos ^{2}\alpha _{\rm max} \sigma _{A_{\rm therm}}^{2} \nonumber \\ &+ A_{\rm therm}^{2}\sin ^{2}\alpha _{\rm max} \sigma _{\alpha _{\rm max}}^{2}+\sigma _{A_{\rm ellips}}^{2}. \end{eqnarray} \tag{ 3 }$$

In practice, this is an overestimate of uncertainty, because 〈F_p/F_*〉 and A_therm are anti-correlated.

IRAC channels 1 and 2 exhibit well-known IPSVs: photons hitting certain parts of a pixel lead to more electron counts than photons hitting other parts of the same pixel (e.g., Charbonneau et al. 2005; Morales-Calderón et al. 2006).¹⁴ In general, the sensitivity to photons is lowest near the edges of a pixel and greatest near its center (see bottom panels of Figure 2). On its own, IPSV would not be a problem for precision time-resolved relative photometry. But over the course of observations, the point-spread function (PSF) of the target star moves on the detector. Even though the PSF spans many pixels, the IPSVs do not average out, because most of the flux falls in the PSF core.

Zoom In Zoom Out Reset image size

Since they are ultimately caused by changes in the PSF position, attempts to correct for IPSVs rely on accurate centroiding, described in Section 2. The centroiding is shown in Figure 2 for 3.6 μm (left) and 4.5 μm (right). The top panels show the centroid jitter and drift over the course of the observations; the second panels show the 2D path of the centroid; the bottom panels show the intra-pixel sensitivity map of the four central pixels of the detector, constructed by applying the Ballard et al. (2010) point-by-point decorrelation to our mapping data (see Section 3.2.2).

Both x and y centroids exhibit fast jitter (period of roughly half an hour) with peak–trough amplitude of 0.05–0.1 pixels. This is the same jitter that used to have a period of roughly an hour: it is related to the heater cycling on Spitzer. The cycling frequency was doubled in fall 2010, which doubled the centroid jitter frequency and roughly halved the amplitude of the jitter.¹⁵ The smaller amplitude of the jitter is undoubtedly an improvement, while the higher frequency may or may not be a nuisance, depending on the duration of ingress/egress for a given planet. In our data, the 3.6 μm flux exhibits clear 1% peak-to-trough flux variations on the centroid jitter timescale, while the 4.5 μm flux does not.

There is also a longer-term centroid drift, which is greatest in the y-direction: 0.5 pixels of motion over the ∼1 day observation at both 3.6 and 4.5 μm. The x-direction shows almost no long-term drift (0.05 pixels, comparable to the faster jitter). Looking at the bottom panels of Figure 2, it is easy to understand why the IPSVs are worse at 3.6 μm than at 4.5 μm: at the shorter wave band, the PSF drifted up a steep slope in sensitivity from a pixel corner toward a center; at the longer wave band, the PSF contoured below a ridge in sensitivity.

The telescope takes a few hours to settle after pointing at a new target, resulting in larger PSF excursions for the first few data cubes of each light curve. It is difficult to correct for IPSVs in poorly sampled regions of the detector, so we remove the first 0.05 days of both the 3.6 and 4.5 μm light curves (3.8% of our data in each wave band).

The crux of the data reduction process is correcting for IPSVs, because our astrophysical signals (eclipses and phase variations) have an expected amplitude of 0.4%, comparable to—or smaller than—these systematics. We tried a variety of techniques, of which we describe the most promising below. We first present two methods for removing the IPSVs before fitting our astrophysical model of the system. These techniques are attractive because they allow one to produce a systematics-corrected light curve independent of any astrophysical model assumptions. We then present an attempt to simultaneously fit the IPSV and the astrophysical brightness variations.

3.2.1. Gaussian Decorrelation

We follow Ballard et al. (2010) in using point-by-point positions and fluxes to generate an intra-pixel sensitivity map with x and y Gaussian smoothing lengths of σ_x and σ_y, respectively. Stevenson et al. (2011) recently introduced a similar—but faster—method using bilinear interpolation. Since we are not attempting to iteratively fit the astrophysical and IPSV model, we use the simpler Ballard et al. (2010) method. We experimented with different smoothing lengths and chose the combinations that yield the smallest χ² value for the final model fit; the astrophysical parameters are not very sensitive to changes in the smoothing length. At 3.6 μm we use σ_x = 0.017 and σ_y = 0.0043, as in Ballard et al. (2010); at 4.5 μm we use σ_x = σ_y = 0.015. The resulting pixel maps are shown in the top panels of Figure 3.

Zoom In Zoom Out Reset image size

We then divide the raw photometry by the weight function and fit our astrophysical model. The second panels of Figure 3 show the corrected data with best-fit astrophysical model and residuals. The bottom panels show the scatter in the residuals as a function of binning, along with a red line indicating the Gaussian-noise limit of root-mean-squared scatter scaling as $\sqrt{N}$. The normalization of this theoretical curve is based on the Poisson error for our electron counts (see first the raw photometry in Figure 1). The best-fit astrophysical parameters are listed in Table 2 (3.6 μm) and Table 3 (4.5 μm).

Table 2. 3.6 μm Parameters

Calibration Method	χ2R	b	a/R*	(Rp/R*)2	Fday/F*	2Atherm	αmax	Aellips
Point-by-point decorrelation	1.392a	0.5(1)	2.8(2)	0.0125(4)	0.0038(4)	0.0004(3)	0(29)°	1(1) × 10−4
Polynomial in x and yb	1.384a	0.3(2)	3.1(2)	0.0123(3)	0.0033(4)	0.0038(6)	−53(7)°	1(1) × 10−4

Notes. ^aWhen comparing these values, it is worth remembering that with approximately 52,000 degrees of freedom, a model is significantly better if it improves χ²_R by at least 0.004. ^bOur fiducial analysis.

Download table as:  ASCIITypeset image

Table 3. 4.5 μm Parameters

Calibration Method	χ2R	b	a/R*	(Rp/R*)2	Fday/F*	2Atherm	αmax	Aellips
Point-by-point decorrelation	1.326a	0.5(1)	2.9(2)	0.0112(4)	0.0039(3)	0.0019(3)	−12(6)°	1.1(1) × 10−3
Polynomial in x and yb	1.324a	0.5(1)	2.9(2)	0.0111(4)	0.0039(3)	0.0040(3)	−16(4)°	1.2(2) × 10−3

Notes. ^aWhen comparing these values, it is worth remembering that with approximately 52,000 degrees of freedom, a model is significantly better if it reduces χ²_R by at least 0.004. ^bOur fiducial analysis.

Download table as:  ASCIITypeset image

3.2.2. Mapping Data

3.2.3. Polynomial IPSV Model

Here, we model the IPSVs as a polynomial in the centroid x and y. We simultaneously fit our astrophysical model and the IPSVs by treating the x and y centroids as independent variables in our function. We model the IPSVs as

$$\begin{equation} \frac{F_{\rm obs}}{F_{\rm astro}} = \frac{1 + \sum _{i=1}^{n} \left[a_{i}(x-\bar{x})^{i} + b_{i}(y-\bar{y})^{i}\right]}{\langle 1 + \sum _{i=1}^{n} \left[a_{i}(x-\bar{x})^{i} + b_{i}(y-\bar{y})^{i}\right]\rangle }, \end{equation} \tag{ 4 }$$

where F_obs is the observed flux, F_astro is the astrophysical model, and $\bar{x}$ and $\bar{y}$ are the mean centroid positions. Formally, cross-terms are necessary to describe an arbitrary 2D function, but we find that including cross-terms does not significantly improve the χ² or affect the astrophysical parameters. This is a testament to the fact that the bulk of the PSF motion is in the y-direction.

We experiment with polynomials up to sixth order. To test whether each additional pair of coefficients (one each for x and y) improved the fit, we use the Bayesian Information Criterion (BIC; Schwarz 1978), which imposes a penalty term on the χ² for additional free parameters: BIC = χ² + kln N, where k is the number of free parameters and N ≈ 52, 000 is the number of data. Since ln (52, 000) ≈ 11, an additional parameter is acceptable if it improves the χ² by at least 11. The BIC for our 3.6 μm data improved with the addition of parameters up to and including fifth order. The BIC for our 4.5 μm data was not improved by the addition of terms beyond cubic.

Unlike some previous studies, we do not include a linear ramp in time. Since the bulk of the PSF motion is a monotonic drift in the y-direction, we found the ramp in time to be highly correlated with the linear and quadratic coefficients of the y sensitivity, and the fit was not significantly improved.¹⁶

We show the resulting fits in Figures 4 and 5. For each figure, the top panel shows: the systematics-corrected light curve and best-fit astrophysical model (top inset), the residuals after subtraction of the best-fit thermal phases, along with the best-fit ellipsoidal variations model (middle inset), and the residuals after removing ellipsoidal variations (bottom inset). (Ellipsoidal variations of the planet do not affect in-eclipse data since the planet is hidden from view; hence we remove the in-eclipse data from this panel.) The bottom left panel shows the weight function used to correct the data. The bottom right panel shows the scatter in the residuals as a function of binning; the red line shows the photon noise limit; the vertical dotted line denotes the timescale of ingress/egress.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Pont et al. (2006) proposed a simple method to account for red noise by considering how the scatter in residuals decreases with bin size (bottom panels of Figures 3–5). At the left end of the plots, where we are considering point-to-point scatter in the residuals, the scatter is only 10%–20% greater than the Poisson counting limit shown in red ($\propto \sqrt{M/N(M-1)} \approx \sqrt{N}$, where N is the number of observations per bin and M is the number of bins). But the observed scatter does not follow the theoretical relation as the data are binned. The most important timescale for transit and eclipse parameter estimation is the duration of ingress/egress, which we denote by a vertical dotted line in those panels (21 minutes for WASP-12b). The scatter on this timescale determines the accuracy we can expect to achieve for transit or eclipse depths.

Following Winn et al. (2007, 2008), we define the factor β as the actual scatter (black line) divided by the theoretical Poisson limit (red line) on the 21 minute timescale (vertical dotted black line); our residuals have β = 2–3. To account for red noise in parameter uncertainties, we simply inflate the L–M parameter uncertainties by β. (Inflating the photometric error bars by β and recomputing the covariance matrix using L–M takes longer and produces slightly smaller parameter uncertainties.) The binning of residuals method turns out to be the most conservative error estimate for transit- and eclipse-specific model parameters (b, a/R_*, (R_p/R_*)², etc.).

We also estimate parameter uncertainties using two resampling techniques: boot-strap and prayer-bead MCs. Both of these techniques use the scatter in the residuals of our best-fit model as an estimate of photometric uncertainty. In both cases the residuals are shifted in time and added back to the best-fit model to produce a new instance of the light curve, which is then fit using the L–M.¹⁷ The standard deviation in the sequence of model parameters is our estimate of their 1σ uncertainty. Note that—by construction—resampling techniques cannot improve parameter estimates: the best-fit parameters are those determined by fitting the original time series.

For the boot-strap MC, the residuals are randomly shuffled so—much like the L–M and MCMC techniques—the boot-strap error analysis is insensitive to the ordering of residuals. Such error estimates will therefore only be accurate insofar as residuals are uncorrelated in time. Indeed, we found that error estimates from the boot strap were comparable to those from the L–M or MCMC analyses.

The prayer-bead analysis maintains the relative ordering of the residuals and simply shifts them all by the same amount (wrapping around the start/end of the data), so that correlated noise present in the residuals is preserved. A prayer-bead analysis is not appropriate if the nature of the noise is expected to change throughout the observations (e.g., the 8 μm ramp seen in cryogenic Spitzer observations). But there is no evidence for such changes in behavior in Warm Spitzer data in general, or in our time series in particular. The prayer bead can only have as many iterations as there are data points, but this is not a problem for the current study given our approximately 52,000 images; we run the prayer bead for 10,000 iterations with randomly chosen offsets and verify that the uncertainties have converged by comparing to 100- and 1000-iteration prayer-bead MCs. It is also worth noting that with such a long data set, we are more likely to observe rare instances of bad behavior in the detector, making the prayer-bead technique particularly conservative. Indeed, prayer-bead MC provides the largest error bars for the phase variation parameters: A_therm, α_max, A_ellips.

Our values for the impact parameter and geometrical factor are broadly consistent with published values. Note that we simultaneously fit the transits and eclipses, and eclipses are notoriously bad at constraining these orbital parameters.¹⁸

Combining our transit times for epochs 925 and 947 (BMJD of 55,518.0407(4) and 55,542.0521(4), respectively) with those of Chan et al. (2011) and the ephemeris of Maciejewski et al. (2011), we obtain a BJD discovery epoch transit center of t₀ = 2454508.9768(2) and period P = 1.0914207(4) days. Since our transits occurred slightly earlier than predicted, our best-fit orbital period is 3.5σ shorter than that of Maciejewski et al. (2011), but it is notoriously difficult to compare transit times at different wavelengths analyzed by different groups because of subtle differences in star-spot coverage, treatment of limb darkening, etc. (Désert et al. 2011b).

More importantly, we are simultaneously fitting an entire orbit's worth of data including many different astrophysical effects, while the optical transit data were fit on their own, irrespective of longer-term astrophysical trends. In order to make a more appropriate comparison, we try fitting the transits independently of the rest of the data (considering only those data within 0.15 days of the transit center). This yields later transit centers, leading to an ephemeris more consistent with that of Maciejewski et al. (2011): t₀ = 2454508.9767(2), P = 1.0914210(4).

Because of WASP-12b's high temperature and inflated radius, the variations in transit depth with wavelength should be 3 × larger than for "typical" hot Jupiters (using the scaling relation from Winn 2010). Based on observations and models of HD 189733b (Fortney et al. 2010), one might therefore expect the transit depth at 4.5 μm to be ∼10⁻³ deeper than at 3.6 μm. This is borne out by the dotted black line in Figure 6, which shows a Burrows et al. (2007, 2008) model transit spectrum assuming solar composition and a day-like temperature–pressure (T-P) profile. If the planet's terminator has a night-like T-P profile (shorter scale height; solid black line in Figure 6), the difference in transit depth could be as small as 2 × 10⁻⁴—but always with a transit depth greater at 4.5 μm than at 3.6 μm.

Zoom In Zoom Out Reset image size

In Figure 6, we compare our transit depths at 3.6 and 4.5 μm with previously published optical values: 0.0138(2) (B and z' band, Hebb et al. 2009), 0.01380(16) (R band, Maciejewski et al. 2011), and 0.0125(4) (V band, Chan et al. 2011), yielding a three-band transmission spectrum of the planet. Curiously, the transit depth at 4.5 μm is considerably shallower than at 3.6 μm. If we adopt the larger Maciejewski et al. (2011) optical transit depth, then our mid-IR transit depths could be indicative of hazes or optical absorbers in the atmosphere of WASP-12b, as has been inferred for HD 189733b (Pont et al. 2008; Sing et al. 2011). The larger radius at 3.6 μm as compared to 4.5 μm indicates a higher atmospheric opacity at the shorter wave band, which is difficult to reconcile with current models.

In an attempt to explain its peculiar eclipse spectrum, Deming et al. (2011) hypothesized that CoRoT-2b might have equal opacity at 3.6 and 4.5 μm (and lower opacity at 8 μm) due to a haze of—as yet unknown—μm-sized particles. One may similarly explain the unusual WASP-12b transit spectrum in terms of a haze of slightly smaller particles, such that the opacity drops from 3.6 to 4.5 μm. Given the large difference in transit depth between the two wave bands, this hypothesis also requires a large atmospheric scale height at the planet's terminator.

At 3.6 μm, the eclipse depth is ∼1σ lower using the polynomial fit as compared to the decorrelation. If we fit the two eclipses individually using the polynomial IPSV fit, we obtain depths of 0.0030 (highly correlated residuals) and 0.0038 (incomplete egress), respectively. Given the prior measurement of 0.0038(1) by Campo et al. (2011), we adopt the larger value from the Gaussian decorrelation for our analysis (this choice does not significantly affect any of our conclusions).

At 4.5 μm, we obtain comparable χ² values and eclipse depths regardless of our treatment of systematics. In both cases they are consistent with the Campo et al. (2011) value: 0.0038(2).

In all cases our error estimates are somewhat larger than the Campo et al. (2011) estimates. We observed the same planetary system with the same instrument, so one expects the same eclipse depths and uncertainties. The only significant change in the detector is that our observations occurred after the fall 2010 change in heater cycling, as mentioned in Section 2. This means that the short-term telescope jitter for our observations has a period of 30 minutes rather than 1 hr, and the amplitude of the centroid excursions—and hence flux variations—is reduced by a factor of two (for reference, the ingress/egress time for WASP-12b is 21 minutes, and the total transit/eclipse duration is 3 hr).

It is conceivable that the near coincidence between the centroid jitter half-period and the ingress/egress timescale leads to greater residual systematics in our data than in the Campo et al. (2011) time series. However, our eclipse depths are based on two occultations, and we have a much longer baseline of observations to help us correct for detector systematics, characterize noise properties, and estimate uncertainties (see Section 4.2). Our MCMC error estimates are similar to those of Campo et al. (2011), but residual-binning and prayer-bead analyses produce eclipse depth uncertainties more than 2 × larger than the MCMC. This means that there is still red noise present in our residuals, and the prayer-bead analysis is a more realistic estimate of our parameter uncertainties.

5.2.1. Dayside Emergent Spectrum

If one assumes solar composition, the relative eclipse depths at 3.6 and 4.5 μm can probe the temperature versus pressure (T-P) profile of the planet. Water vapor absorbs less at 3.6 μm than at 4.5 μm, so the shorter wave band will have a higher brightness temperature if the temperature is locally dropping with height or vice versa (e.g., Burrows & Orton 2010).

If the eclipse depth is measured at sufficiently many wavelengths, one may hope to simultaneously constrain a planet's atmospheric composition and T-P profile (e.g., Madhusudhan & Seager 2009). The high C/O chemistry invoked by Madhusudhan et al. (2011) was the result of an abnormally low eclipse depth at 4.5 μm in the Campo et al. (2011) data, which was interpreted as being due to CO absorption in the planet's relatively cool upper atmosphere. More recently, Kopparapu et al. (2011) used a photochemical model to study the disequilibrium chemistry of WASP-12b, confirming that CO would be enhanced in a high C/O composition atmosphere.

In Figure 7, we compare the near- to mid-IR broadband spectrum of WASP-12b to various one-dimensional (1D) radiative transfer models (Burrows et al. 2007, 2008). We vary the abundance of CO as a proxy for varying the C/O ratio. Our eclipse depths are consistent with those of Campo et al. (2011), so we still favor models with enhanced CO (10 × solar) and a weak inversion for this planet, in agreement with Madhusudhan et al. (2011). We also find that we can obtain an equally good fit to the data by reducing H₂O to 1% solar abundance and partitioning carbon evenly between CO and CH₄. Given the caveat that these models are in radiative—but not chemical—equilibrium, the crux of fitting WASP-12b's unique dayside spectrum is to suppress H₂O with respect to CO.

Zoom In Zoom Out Reset image size

Our larger error bars, however, make these composition statements marginal: the solar composition model with modest inversion (the green line in Figure 7) is only worse by Δχ² = 8 for 8 degrees of freedom. Furthermore, neither the standard composition nor the enhanced CO scenario is consistent with the relative transit depths at 3.6 and 4.5 μm (see Figure 6 and previous section). It may be possible to reconcile these two measurements if the atmospheric composition is grossly different at the day–night terminator than near the sub-stellar point.

The best-fit ellipsoidal variations at 3.6 μm are consistent with the predicted amplitude of 2 × 10⁻⁴, but as one can see from the relative uncertainty (or from glancing at Figure 4), they are not robustly detected: the χ², residuals and remaining astrophysical parameters do not change significantly if A_ellips is set to zero.

As shown in Figure 5, however, we clearly detect power at 4.5 μm in the second cosine harmonic, cos (2α), consistent with the prediction of ellipsoidal variations due to the prolate shape of the planet (Li et al. 2010; Leconte et al. 2011; Budaj 2011). The semi-amplitude of the variations, however, is 1.2(2) × 10⁻³ using either decorrelation or polynomial IPSV-removal, approximately 6 × the predicted value.

Since ellipsoidal variations are primarily a geometrical effect, it is difficult to understand how the measured amplitude could be so different at 3.6 and 4.5 μm. The upper layers of the atmosphere should be more distorted; the deeper layers more spherical. The relative strengths of ellipsoidal variations at the two wave bands imply that the 4.5 μm flux is originating from much higher up in the planet's atmosphere than the 3.6 μm flux. The simplest way to do this is for the atmosphere to have a greater opacity at 4.5 μm than 3.6 μm. But the relative transit depths indicate exactly the opposite, as described above and shown in Figure 6.

Alternatively, it is possible that detector systematics still present after IPSV-removal attenuate the ellipsoidal signal at 3.6 μm or enhance it at 4.5 μm. The raw photometry shown in Figure 1 implies that the 4.5 μm light curve is more trustworthy of the two. We therefore begin by assuming that the ellipsoidal signal at 4.5 μm is entirely astrophysical in nature, then consider the opposite scenario.

5.3.1. Interpreting the Ellipsoidal Variations at 4.5 μm

If we take the 4.5 μm ellipsoidal variations at face value, they have surprising implications for the planet's shape. The dimensions of a prolate planet may be described by its short and long radii, denoted by R_p and R_long, respectively. The third dimension (parallel to the system's angular momentum vector) is assumed to be R_p because we are neglecting rotational effects, which tend to produce oblate planets. Note that WASP-12b is on a very short period orbit, so—if tidally locked—it has a rotation rate only a factor of two slower than Jupiter or Saturn. As a result of this rotation, Budaj (2011) estimates the planet's polar radius to be 2.5% shorter than its lateral equatorial radius,¹⁹ so strictly speaking WASP-12b is a triaxial ellipsoid, but this does not affect the analysis below because of the planet's edge-on orbit, and our relatively short baseline of observations is insensitive to the spin precession of the planet (Carter & Winn 2010b).

At conjunction—either transit or eclipse—we are seeing the planet's smallest projected area (πR²_p, in the case of a perfectly edge-on orbit). The projected area of an ellipsoid on an edge-on orbit varies as (e.g., Vickers 1996)

$$\begin{equation} A_{p}(\alpha) = \pi R_{p}^{2} \sqrt{\cos ^{2}\alpha + \left(\frac{R_{\rm long}}{R_{p}}\right)^{2}\sin ^{2}\alpha }. \end{equation} \tag{ 5 }$$

For R_long/R_p ≈ 1, the changes in projected area follow a A_p∝cos (2α) shape and this is in fact how we modeled them; for more severe elongations, the peaks become broader and the troughs narrower, until a limiting case of A_p∝|sin α|.

If we interpret all of the power in the second cosine harmonic as being due to the changing cross-sectional area of the planet (the red curves in Figure 8), we may estimate the planet's aspect ratio as

$$\begin{equation} R_{\rm long}/R_{p} \approx 1 + 2A_{\rm ellips}\left\langle F_{*}/F_{\rm day}\right\rangle = 1.8(1). \end{equation} \tag{ 6 }$$

The predicted aspect ratio for the planet is R_long/R_p = 1.1, while the aspect ratio for a Roche lobe is 3/2 = 1.5.

Zoom In Zoom Out Reset image size

It is worth noting that a planet with an aspect ratio of 1.8 would have a 13% larger projected area at the start and end of transit as compared to transit center, resulting in a w-shaped transit, in the absence of stellar limb darkening (note that this differs from changes in transit morphology due to an oblate planet discussed in Carter & Winn 2010a). Since the transit depth of WASP-12b is approximately 1%, this shape-induced transit effect would come in at the 1.3 × 10⁻³ level and might be detectable in our current data. In Figure 9, we estimate the expected transit morphology by treating the planet as a sphere with variable cross-sectional area. The red line shows our fiducial model: a spherical planet with nonlinear limb darkening of the star. The green line shows the effect of the planet's changing cross-sectional area (but not its changing shape). The blue line shows how limb darkening partially washes out this signal. Figure 9 implies that our data rule out the most extreme prolate toy model, but clearly the treatment of limb darkening is important here.

Zoom In Zoom Out Reset image size

If we instead interpret the phase curve as being caused entirely by longitudinal brightness variations on a spherical planet (the blue curves in Figure 8), it can be inverted into a longitudinal intensity map of the planet, following Cowan & Agol (2008). Because of the low-pass filtering that occurs in the map → light curve convolution, the deconvolution will enhance the highest frequency terms present in the light curve (e.g., Figure 2 of Cowan & Agol 2009). In the current case, the resulting map has two prominent temperature peaks: one at the dawn terminator and one at the dusk terminator. More importantly, the only way to simultaneously fit the bright terminators and dark night side is by having negative intensity at the anti-stellar point, clearly an unphysical solution.

Another argument against a spherical planet is that while the thermal phase variations of a spherical planet will in general contain power in the second harmonic, there is no reason for it to all appear in the cosine rather than sine term.²⁰ When we run fits allowing for the phase of the cos (2α) term to vary, the offset is consistent with zero at the 2σ level, with the largest offset being less than 6°. This indicates that the flux is peaking at quadrature, as expected for ellipsoidal variations. To our knowledge, there is no reason to expect such a temperature profile if the planet is spherical.

If the planet is prolate, however, one might expect something akin to the gravity darkening/brightening seen in some binary stars: the sub-stellar and anti-stellar points on the planet are farther from the center of the planet than the terminator, leading to lower surface gravity. The lower surface gravity at the sub-stellar and anti-stellar regions might lead to cooler temperatures, all things being equal (von Zeipel 1924). But this would only affect the intrinsic component of the planet's power budget, expected to be insignificant for hot Jupiters.

Finally, it is possible that the planet's shape affects the circulation of its atmosphere in more subtle ways, leading to relatively hot regions at the dawn and dusk terminators. This brings us to our favored astrophysical interpretation of the ellipsoidal variations. The planet may have an aspect ratio of R_long/R_p ≈ 1.5 (somewhat easing the transit morphology constraints described above), with an additional enhancement of the cos (2α) power because of a relatively hot day–night terminator (the green curve in Figure 8). Since the atmospheric dynamics on severely non-spherical planets has not yet been addressed in the literature, it is difficult to say whether this scenario is reasonable.

5.3.2. The Null Hypothesis at 4.5 μm

We model thermal phases using only first harmonics (cos α and sin α), while our parameterization of ellipsoidal variations is a second harmonic (cos 2α). These functions are by definition orthogonal, so it is not surprising that our conclusions about thermal phase variations (A_therm, α_max) described below are not significantly affected by the presence or absence of ellipsoidal variations discussed above.

That said, the amplitude and offset we obtain for the thermal phase variations depends on which IPSV-removal scheme we use. Unfortunately, it is not clear how to perform a direct model comparison between decorrelation and polynomial fits using the BIC. The Gaussian decorrelation can be thought of as having a large number of free parameters and a somewhat smaller number of additional constraining equations, but it is not obvious how to estimate its degrees of freedom. As discussed by Ballard et al. (2010), the point-by-point decorrelation does a great job of removing the short-term jitter and long-term detector drift, but may also remove any longer-term astrophysical signal. (This did not interfere with their goal of searching for transits.) Insofar as diurnal phase variations are the most gradual astrophysical signal in our study, it is not surprising that it is the most dependent on IPSV-removal.

The polynomial fit leads to χ² values at least as good as—and sometimes significantly better than—the Gaussian decorrelation. More importantly, the scatter in the residuals on the critical 21 minute timescale is lower for the polynomial fits, indicating that this method is doing a better job of removing correlated noise.

Our 4.5 μm pixel sensitivity map obtained from Gaussian decorrelation (top right panel of Figure 3) shows a valley at y ≈ 14.75; this does not follow the usual pattern of sensitivity decreasing toward pixel edges. (Note that we observe the same ripples in sensitivity as Ballard et al. 2010, but those occur on a much smaller spatial scale—and exhibit a much smaller amplitude—than the y ≈ 14.75 valley.)

Finally, the middle left panel of Figure 3 suggests that the decorrelation method has overcorrected the 3.6 μm system flux near superior conjunction: the eclipse bottoms—which should be flat since the planet is hidden from view—slope upward toward the central transit.

We therefore argue that the decorrelation has filtered out much of the phase variations along with the systematics and adopt the polynomial-fitted thermal phase parameters in what follows.²¹

Following Cowan & Agol (2011b), we estimate the hemispheric effective temperatures to be T_day = 2928(97) K and T_night = 983(201) K, where the uncertainties include an estimate of systematic errors in going from brightness temperatures to effective temperatures.²² This implies a Bond albedo of A_B = 0.25 (presumably due to Rayleigh scattering and/or reflective clouds) and very low heat recirculation efficiency, ε < 0.1, at 1σ (see Figure 10).²³ Note that our 1D radiative transfer models used for interpreting transit and eclipse spectra are cloud-free and have an albedo lower than this inferred value.

Zoom In Zoom Out Reset image size

Space-based optical secondary eclipse depths have been measured for a handful of hot Jupiters. If the planet's equilibrium temperature is sufficiently low, such measurements provide an unambiguous estimate of geometric albedo, A_g: 0.04(5) for HD 209458b (Rowe et al. 2008), 0.32(3) for Kepler-7b (Demory et al. 2011), 0.10(2) for Kepler-17b (Désert et al. 2011a), 0.30(8) for KOI-196b (Santerne et al. 2011), and 0.025(7) for TrES-2b (Kipping & Spiegel 2011). Acknowledging that WASP-12b is more than 1000 K hotter than any of those planets, a Bond albedo of 0.25 is well within the observed range for hot Jupiters.

Assuming gray albedo and a Lambertian scattering phase function (A_B = 3/2A_g; Hanel et al. 1992), the reflected-light secondary eclipse of WASP-12b should have a depth of 2.4 × 10⁻⁴, comparable to current ground-based precision for this target. In practice, most planets exhibit an opposition surge (making them disproportionately bright at superior conjunction) and Rayleigh scattering, so the actual contrast ratio in blue optical wave bands is likely more favorable than this estimate.

Assuming solar atmospheric composition and hence opacity, the 3.6 μm thermal flux should originate from deeper in the atmosphere than any other mid-IR wave band. Insofar as radiative times increase monotonically with pressure, we may therefore expect the 3.6 μm phase variations to be muted compared to the 4.5 μm phase variations, and the phase offset should be greater at the shorter wavelength (e.g., Figure 9 of Burrows et al. 2010).

The hot spot offset is 53(7)° east of the sub-stellar point at 3.6 μm. While not as extreme as υ-Andromeda b (Crossfield et al. 2010), it is difficult to reconcile our large phase offset and large amplitude at 3.6 μm.²⁴ It is worth noting, however, that there are highly correlated residuals near the purported peak (∼0.35 day after transit) which may be partially responsible for the large offset.

On the other hand, we find that the 4.5 μm phase amplitude and offset, 16(4)°E, are consistent with a Cowan & Agol (2011a) model with heat transport efficiency of ≡ τ_rad/τ_adv ≈ 0.1. To put this in context, phase variations and eclipse timing offset of HD 189733b at 8 μm indicate ≈ 0.7 (Agol et al. 2010).

5.4.1. Implications of Thermal Phase Variations