INFRARED ECLIPSES OF THE STRONGLY IRRADIATED PLANET WASP-33b, AND OSCILLATIONS OF ITS HOST STAR

We observed secondary eclipses of WASP-33b on UT 2011 October 10 and 16, using the FLAMINGOS infrared HgCdTe imager at the Cassegrain focus of the 2.1 m telescope at Kitt Peak National Observatory. The sky on both nights was cloudless, with excellent photometric conditions prevailing, especially on October 10. The observations used a K_s filter, and we defocused the telescope so that the diameter of stellar images—measured at the half-intensity level—was 30 pixels (18''). This substantial defocus improves the photometric precision and photon-collection efficiency for this bright star (WASP-33 has V = 8.3). We compensated for drift in the telescope defocus using manual updates at approximately 30 minute intervals, based on a known formula that relates focus position to temperature and zenith distance. In that way, we attempted to maintain the greatest possible image stability.

The 20 × 20 arcmin field of view of FLAMINGOS (2048 × 2048 pixels) provided five comparison stars imaged simultaneously with WASP-33. We obtained a nearly continuous sequence of 20 s exposures on each night, amounting to 725 exposures on October 10 and 574 exposures on October 16. Including the overhead of reading the detector and writing FITS files, the observational cadence was 45 s per exposure. We verified that the times written to the FITS headers were free of clock errors (to about 1 s precision) by comparing the header values to manual timing made using a web-displayed UT clock. In addition to WASP-33, we observed multiple sky exposures with position offsets that were used to construct a flat-field frame by median-combining the offset sky exposures. We also observed dark frames using the same exposure times as for the stellar and sky-flat observations.

All of our WASP-33 observations used a quad-detector off-axis guide camera, sensitive to optical radiation, and producing real-time pointing corrections for the telescope. The guider was very effective at damping image motion on timescales of minutes, but differential refraction between the optical and the infrared leads to a slow drift in stellar positions, amounting to about 1'' over 60° of zenith distance. This slow positional drift has only a small effect on our photometry.

During the observations, we noted a significant instability in the infrared signal from the FLAMINGOS detector, which occurred in response to the changing position of the telescope. In order to diagnose and correct for these instabilities, we obtained a sequence of K-band exposures during the day, with the dome closed and dark, and the primary mirror cover closed. These 2 s exposures measured the thermal emission and scattered light from the telescope mirror cover, and we moved the telescope to different positions because the signal instabilities seemed to be a function of telescope position, as described in Section 3.2.

During our WASP-33 observations, we also noticed relatively prominent intensity oscillations of the star in K_s band. Hence, we also observed 651 consecutive exposures of WASP-33 on the night of UT October 14, when no eclipse (or transit) of the planet occurred. This time series was observed using J band and helps to establish the degree to which the amplitude of stellar oscillations in intensity may depend on wavelength.

Warm Spitzer observed one eclipse of WASP-33 at 3.6 μm, for 9.6 hr on 2011 March 26, as well as observations of equal duration spanning one eclipse at 4.5 μm on 2011 March 30. The observations were made under the Cycle-6 Target of Opportunity program (PI: J. Harrington). Both observations used subarray mode to collect 1264 data cubes, each data cube comprising 64 exposures of a 32 × 32 pixel section of the detector. The exposure time was 0.4 s per exposure, and the observations were not interrupted by any re-pointing.

In both K_s and J bands, we subtract a dark frame from each image, and divide by the sky flat for that wavelength. Each sky flat is normalized to have an average value near unity; this facilitates conversion of the observed images from data numbers to electrons using the known gain of the detector electronics (4.9 e/DN). In the limit of large defocus, the stars are essentially images of the telescope's primary mirror; they are not well approximated using Gaussians or other functions commonly used for centroiding. Therefore, we determine the center position of each star by exploiting the sharp edges that are characteristic of the pupil images. Summing the image of each star in the X-coordinate produces intensity as a function of Y, I(Y). The sharp edges of the pupil image will produce peaks in the derivative ∂I/∂Y, one peak at each edge of the image. We find those peaks in the spatial derivative, and adopt an average of the Y-positions of derivative peaks on each side of the I(Y) profile as being the center of the star in Y, and vice versa for X. With the center of each star determined for each exposure, we calculate the flux in a circular aperture of a given radius centered on that star.

In addition to WASP-33, we measure 4–6 comparison stars on each image. The identification of the comparison stars are given in Table 1, with their J and K magnitudes. The subset of comparison stars that were actually used varied from night to night, due to using the J-band wavelength on October 14, and K-band sky transparency and thermal background that varied between October 10 and 16. The subset of comparison stars used for each night are noted in the caption of Table 1.

For both the target and comparison stars, we determined the value of the sky background by constructing a histogram of pixel values in a 100 × 100 pixel box surrounding each star. The sky background dominates the peak of those histograms, and we fit a Gaussian to each histogram, and thereby determine the sky background value for each star as the centroid of each Gaussian. Multiplying the background value per-pixel for each star times the area of the aperture containing that star yields the background contribution for that star. Subtracting the background from each aperture measurement yields the stellar intensities of WASP-33 and the comparison stars for that image.

We repeat our photometry by varying the radius of the synthetic aperture from 18 to 30 pixels, in 1 pixel increments. The 13 different aperture radii produce 13 different realizations of photometry for each night. Our rationale for varying the aperture radius is to optimize the signal to noise (S/N) and the match between the target star and the comparison stars. Although our images are strongly defocused, there is still scattered light beyond the edges of the defocused pupil image for each star. Variable seeing and errors in aperture centering cause the total flux intercepted by a chosen aperture to vary. Measurements with larger apertures are less sensitive to these variations, but include more background noise. The optimum aperture radius is approximately 20 pixels, but we determine it separately for each night, as described in Section 4.1.

As noted in the Introduction, we experienced instabilities in the photometric signals. These instabilities were discovered using quick-look evaluation of the data during each night of observations. The nature of the instabilities is that sharp increases or decreases in stellar intensities occur typically 4–6 times each night. These changes in intensity are large compared to our photometric precision: Sudden changes as large as 4% were seen. Stars close together in angular extent experienced similar, but not identical, effects. For our observations, the comparison stars were typically several hundred pixels distant from WASP-33, so the instabilities were not precisely common mode. The sky background also exhibited this instability, but to a much smaller relative degree—barely detectable.

During our observing run, we noticed that the signal instabilities never occurred for stars at negative declination, and they tended not to occur at large hour angles. We concluded that these puzzling instabilities were triggered by motion of the telescope, and that motivated our diagnostic observations described in Section 2. The root cause of the signal instabilities is now known: following our observing run, Dick Joyce of the Kitt Peak scientific staff disassembled FLAMINGOS, and discovered that a 4-40 screw (probably from the filter box) had come loose and fallen onto the first camera lens, producing field-dependent vignetting as it shifted position. We were initially tempted to discard all of these "screwy" data, but the excellent photometric quality of the nights (especially on October 10) motivated us to develop a methodology to remove the effect of the rolling screw.

Fortunately, the nature of this effect—a sudden shift in position of the screw followed by periods of stability—makes it amenable to robust correction using only the data themselves. We used our test observations with the dome closed to validate our correction methodology. We produce aperture photometry from these frames using apertures of 20 pixel radius, and no background subtraction (the background is the signal for the test observations). We center the apertures at the same locations as WASP-33 and the comparison stars, and extract time series signals. Figure 1 (top panel) shows the time series at the WASP-33 position. The large and sudden changes in signal level are obvious, and they correspond to times of telescope motion. The telescope motion for these test observations was done in a series of 0.5 hr increments, versus continuous tracking for the actual stellar observations. Nevertheless, the results are very similar, and we conclude that we have successfully reproduced the signal instabilities.

Zoom In Zoom Out Reset image size

We correct these instabilities by operating on the time derivative of the signal. We calculate the distribution of the ensemble of the numerical time derivative values for each time series. This distribution is primarily Gaussian, with extreme outliers that correspond to the times of signal instability. We fit to the Gaussian core of the distributions, to determine the unbiased standard deviation of the signal derivative at the position of each star. We adopt a threshold value (typically 3σ) and we correct the outlying points in each series of derivative values, where the absolute value of the derivative exceeds this threshold. This threshold value (in σ) is an adjustable parameter in our correction procedure. After identifying the derivative points beyond the threshold, one option is to set these outlying derivative points to zero. However, in practice we replace derivative outliers with a 5 minute smoothed version of the derivative. We choose the 5 minute smoothing time because it represents the typical time for significant changes in stellar intensity, e.g., as caused by the telluric atmosphere. After correcting the outlying points in the derivatives, we integrate the corrected derivatives to yield corrected time series.

Figure 1 (upper panel) shows the corrected time series for our dome test observations at the WASP-33 position. The derivative of the time series is shown in the middle panel, with the rejected outliers marked. The lower panel of Figure 1 shows the result of dividing the corrected time series at the WASP-33 position by the total of the corrected series at the comparison star positions. In the absence of our correction procedure, this quotient would show considerable noise at the several-percent level. The correction procedure successfully produces a quotient whose short-term variations are of order 200 parts per million (ppm), with long-term variations approaching 700 ppm. We do not regard the long-term variations as meaningful to our correction procedure because the background in the telescope dome is not guaranteed to be stable at this level. We conclude that our procedure to correct for the signal instabilities is effective, and we apply it to our stellar photometry.

Photometry of Warm Spitzer data is now a familiar exercise (Hebrard et al. 2010; Beerer et al. 2011; Deming et al. 2011a; Demory et al. 2011; Desert et al. 2011; Todorov et al. 2012), and we used well-tested procedures applied to the BCD files produced by version S18.18.0 of the Spitzer pipeline. We performed aperture photometry on each frame, in each data cube, for both eclipses. We found the centroid of the stellar image by fitting a two-dimensional Gaussian, and computed the flux in a circular aperture centered on the star, for aperture radii from 2.0 to 5.0 pixels in 0.5 pixel increments. We accounted for the background contribution in each frame; the per-pixel background was measured as the centroid of a Gaussian fit to a histogram of pixel values for that frame. Multiplying the background per-pixel times the area of the aperture yielded the background contribution that was subtracted from the flux in the aperture. Aperture radii near 3.0 pixels produced the highest S/N photometry, as judged by the point-to-point scatter. Accordingly, we adopted a 3.0 pixel radius for our photometry at both 3.6 and 4.5 μm.

4.1.1. Optimization of Photometric Parameters

4.1.2. Normalization Using Comparison Stars

Following selection of the best photometric aperture radius, and threshold for instrumental correction, we further correct the WASP-33 time series using the comparison stars. Normally this would be accomplished by dividing WASP-33 by the sum of the comparison stars. However, we find improved results using a slightly different procedure: Instead of a straight sum of the comparison stars, we use a weighted sum. We denote the intensity of WASP-33 at time index i as W_i, and we write

$$\begin{equation} W_i = \sum _{j=1}^{N} \alpha _j c_{ij}, \end{equation} \tag{ 1 }$$

where c_ij is the intensity of the jth comparison star at time index i, and the α_j coefficients—one for each of the N comparison stars—are determined by linear regression (matrix inversion). The linear regression seeks to produce the best overall match between the weighted sum of the comparison stars and WASP-33, i.e., to make Equation (1) be exact. Because of noise and intrinsic variations in WASP-33, Equation (1) can never be solved exactly, only for the set of α_j that produces the best approximation. Having solved for those α_j, we divide the W_i by the right-hand side of Equation (1). This division by the weighted sum of the comparison stars removes instrumental and telluric effects, but leaves noise and the intrinsic variations of WASP-33 itself.

Equation (1) is a generalization of the usual methodology of ratioing the target star to the sum of the comparison stars. Using an equally weighted sum of the comparison stars can be an imperfect divisor, for several reasons. First, the comparison stars usually differ in spectral type from the target star, and the integral of their flux over our broad filter produces a slightly different effective wavelength for each star. The extinction in the stellar signals caused by the terrestrial atmosphere can vary strongly with wavelength, so different effective wavelengths can degrade the correlations between WASP-33 and the comparison stars. Second, the comparison stars lie at quite different locations on the detector, and higher order effects can be different at these different locations. For example, our data exhibit a slow drift in the position of all stars over the course of a night (about 1'' total), caused by differential refraction between the optical guider and the infrared wavelength used for imaging. That small drift, combined with small inaccuracies in the flat-field calibration, can produce different slow variations for each comparison star. Also, slight changes in the degree of defocus can result from inaccurate compensation of thermal effects and mechanical flexure (Section 2.1), and can produce subtle changes in the defocused images—hence in the aperture photometry—as a function of field position. For these reasons, we believe it is good practice to optimally weight the comparison stars using linear regression.

We compared our optimally weighted photometry to photometry that used an unweighted sum of the comparison stars. The point-to-point scatter in the two photometric data sets was about the same, but the baseline was much flatter using the optimal technique. Specifically, we found that the unweighted sum produced overall slopes of 0.6% and 0.8% for the nights of October 16 and 10, respectively. Using the optimal technique reduced these baseline slopes to 0.03% on both nights.

The linear regression used to solve Equation (1) can potentially affect the measured eclipse, because the regression attempts to remove all fluctuations in the target star, and the eclipse is a fluctuation. In the hypothetical case where one of the comparison stars happened to exhibit a noise-like decrease at the expected time of the WASP-33 eclipse, the regression would overweight that star and the result would be to weaken or remove the eclipse. We avoid that possibility by solving for the α_j using only out-of-eclipse portion of the data (adopting the central phase as 0.5, and setting the duration of eclipse to equal the duration of the transit).

Figure 2 shows the photometry for WASP-33 in the K_s band on both October 10 and 16, and the weighted sum of the comparison stars (right-hand side of Equation (1)) is overplotted as a blue line. The overall rise and fall of these signals is due to the airmass dependence of telluric absorption, and shorter term fluctuations due to telluric effects are also visible. Note also that the weighted sum of the comparison stars shows more short-term noise than does WASP-33b because most of the comparison stars are 1–2 mag fainter than WASP-33 (see Table 1). For that reason, we smooth the comparison star sum over 10 points (about 7 minutes in time) before using it to correct WASP-33. That degree of smoothing decreases the point-to-point noise in the ratio, while still following most telluric fluctuations.

Zoom In Zoom Out Reset image size

A dominant effect in Spitzer photometry at both 3.6 and 4.5 μm is the presence of intra-pixel sensitivity variations. The photometric intensity of a star will depend on its position on the detector, and therefore will vary with time because of a pointing oscillation in the telescope (Carey et al. 2010; see Todorov et al. 2012 for a recent example). The Spitzer project recently implemented software updates that decrease the amplitude of the pointing oscillation, and also decrease its period from 1 hr to about 40 minutes.¹¹ This reduces the impact of the intra-pixel sensitivity variations on the photometry. Indeed, our observations of WASP-33 show the intra-pixel effect to a much less degree than many previous observations.

Many previous investigations have established that the intra-pixel effect is more dependent on the Y-coordinate than on the X-coordinate, and it is stronger at 3.6 μm than at 4.5 μm (Knutson et al. 2009; Beerer et al. 2011; Deming et al. 2011a; Todorov et al. 2012). We see the intra-pixel effect in our 3.6 μm photometry, and we corrected it using a quadratic fit to the photometry as a function of the Y-coordinate of the image, with a linear term as a function of the X-coordinate (the variation with X is weaker than for Y). However, we cannot detect any intra-pixel sensitivity variations in our 4.5 μm photometry. Plots of the measured intensity of the star versus both the X- and Y-coordinates of the image centroid are essentially scatter plots (not illustrated), with no significant trends. Pearson correlation coefficients of intensity versus X- and Y-coordinates have values near zero. We also calculated the Pearson correlation coefficients for temporal subsets of the 4.5 μm data, chosen using a moving boxcar window of various widths. We can find no significant correlations between our 4.5 μm photometry and spatial coordinates, during any time period. As an additional check, we repeated our photometry using an alternative method to determine the position of the star on the detector (center of light). This method also failed to reveal correlation between position and intensity. We therefore use our 4.5 μm photometry for eclipse analysis without applying any spatial decorrelation. Note that we do see very small variations in the photometry at the 40 minute period corresponding to the telescope oscillation (see Section 4.4.2). We cannot rule out the possibility that the spatial-intensity correlation is being obscured by the stellar intensity oscillations.

In the case of the 3.6 μm photometry, we also used the spatial decorrelation method described by Ballard et al. (2010). That weighting function method commonly uses a time threshold to zero-weight points that lie close in time to any given point. When the only intrinsic temporal variation is an eclipse or transit, the time threshold is straight forward to implement. However, when continuous stellar oscillations are part of the desired signal, the weighting function time threshold can be problematic. We applied a weighting function to the 3.6 μm photometry, using a time threshold of zero, i.e., not excluding any points in the calculated weights. Fitting to this alternative version of the 3.6 μm photometry produced results that differed insignificantly from our quoted results (0.34σ and 0.03σ differences in central phase and eclipse depth, respectively). Our quoted results (see Section 4.4.2) are based on the polynomial decorrelation described above because we believe that method is better suited to the nature of these data.

After decorrelating (or not) the intra-pixel effect, we omit some of the initial data for both Spitzer wavelengths. The 3.6 μm data exhibit an initial transient drift in the image position, amounting to about 0.17 pixels, about four times as large as the peak-to-peak variation caused by the 40 minute telescope oscillation. This positional instability is accompanied by a similar large transient effect in the photometry. We are familiar with the nature of these large transient effects based on seeing them in other Spitzer data. The relation between intensity and X, Y position is not the same during the transient as during the stable portion of the time series because the image traces a different region of the pixel. Therefore, we omit the first 83 minutes of the 3.6 μm data—this being the time for the position of the image to stabilize. We see no obvious transient effects in the 4.5 μm image position, but the first 27 minutes of the photometry are anomalously low. As a precaution, we omit those initial data from our 4.5 μm analysis.

Figures 3, 4, and 5 show the ground-based WASP-33 photometry after all corrections. WASP-33 is known to be an oscillating δ-Scuti star (Collier-Cameron et al. 2010; Herrero et al. 2011). Herrero et al. (2011) find a dominant oscillation period of 68 minutes, but oscillating stars can exhibit different oscillation modes that rise and fall in amplitude. We have therefore decomposed the photometry on each night using wavelets (Torrence & Compo 1998) to account for possible intra-night changes in the oscillatory spectrum. Because the eclipse can significantly affect the wavelet analysis, we remove the eclipse prior to computing the wavelet power spectra. For the removal, we adopt the best-fit eclipse parameters as derived in Section 4.4. (This is not a circular procedure because the fits for the eclipse parameters converge relatively independently of the starting estimates for the oscillation periods.)

Zoom In Zoom Out Reset image size

Our results show the stellar oscillations very prominently, but (interestingly) to varying degrees. The J-band results are shown in Figure 3; the bottom panel shows the corrected photometry and the overlaid blue line shows part of the wavelet decomposition—just the lowest order components to guide the eye to the general trend. (Note that we do not use the blue lines on Figures 3– 5 in our analysis, but we do use the results of the wavelet power spectra, shown in the upper panels.)

The power spectrum in Figure 3 peaks at 146 minutes (frequency 9.9 cycles day⁻¹), close to twice the period of the 68 minute oscillation (21.2 cycles day⁻¹) seen by Herrero et al. (2011). There was no mention of this lower frequency mode in the work of Herrero et al. (2011), but δ-Scuti stars are known to exhibit multiple periods of oscillation (Balona & Dziembowski 2011), wherein the non-radial modes cluster around different radial modes (e.g., n = 1, 2, 3, ...; Breger et al. 2009). Hence, detection of a different mode frequency than Herrero et al. (2011) is not surprising. Figure 4 shows the results in similar format for the K_s-band observations on October 10. In this case a very prominent oscillation is seen with a power spectral peak at 71 minutes (20.3 cycles day⁻¹), in good agreement with essentially the same oscillation period (68.56 minutes) that was prominent in the results of Herrero et al. (2011). Note also that it can be clearly discerned in the photometry even prior to the correction using the comparison stars—see Figure 2. The amplitude of this oscillation—about 0.2% in flux—is about twice what was observed by Herrero et al. (2011) in the optical, and we discuss possible reasons for this difference in Section 5.

Zoom In Zoom Out Reset image size

Figure 5 shows the analogous result for the K_s band observed on October 16. In this case we see two peaks in the power spectrum at 54 and 126 minutes that are probably due to the stellar oscillations as seen on October 10 and 14. For each night, and also for the Spitzer data, we estimate the total stellar oscillation amplitude from the power present in the oscillatory component of the Markov Chain Monte Carlo (MCMC) eclipse fit (see Section 4.4), and we tabulate those amplitudes in Table 2. All of these amplitudes are in reasonable agreement with the 0.1% oscillation seen in the optical by Herrero et al. (2011), except in the K_s band where the amplitude is about twice as great, consistently on both October 10 and 16.

Zoom In Zoom Out Reset image size

Variability in stellar oscillation amplitudes is sometimes observed in δ-Scuti stars (e.g., Breger & Pamyatnykh 2006), but is not the most likely explanation for our results, because of our J-band data. Although the J band shows a prominent oscillation near phase 0.82 on Figure 3, quantitative calculation of the average amplitude over that entire night gives 0.16% (Table 2), only marginally higher than the optical value. A close-to-normal oscillation amplitude in J band, between the two nights of higher amplitude in K_s band, seems too much of a coincidence to attribute to variability. While we cannot strictly reject mode variability, we hypothesize that oscillation amplitude is associated with wavelength. We discuss possible wavelength-related explanations in Section 5.2.

4.4.1. K_s-band Eclipse

The eclipse of the planet is visible near phase 0.5 in the lower panels of Figures 4 and 5. No eclipse is seen in Figure 3, as expected for that range of orbital phase. The stellar oscillations are a very significant source of confusion as regards measuring the properties of the eclipse. That confusion is mitigated by sufficiently long duration of our observed sequences (7.3 hr on October 10 and 5.4 hr on October 16) compared with the periods of oscillation, and by the fact that the stellar oscillations will not be in phase on the two nights, whereas the eclipse should repeat at the same orbital phase. In order to exploit the latter advantage, we combine the data for October 10 and 16 into one binned data set, using a bin width of 0.001 in phase, and we fit models to these binned data. Figure 6 shows the combined and binned data, with best-fit eclipse plus oscillation curves overplotted.

Zoom In Zoom Out Reset image size

To obtain the best-fit eclipse depth, we fit a combination of eclipse curve plus two sinusoidal oscillations, using an MCMC method, with Gibbs sampling (Ford 2005). We fit to the combined and binned data of Figure 6, in order to gain the advantage of canceling the stellar oscillation as much as possible. We generate four chains of 10⁶ steps each, and we discard the first 2 × 10⁵ steps in each chain. We initialize the MCMC fit using sinusoids with periods of 68 and 146 minutes, guided by our wavelet power spectra in Figures 3– 5. The MCMC fit is insensitive to the exact choice of initial periods. As long as the initial values contain one period in the ∼70 minute range, and one period in the ∼140 minute range, the Markov chains rapidly converge to the best-fit values for that particular eclipse (this is also true for the Spitzer eclipses; Section 4.4.2).

In order to derive accurate MCMC posterior distributions, it is essential to re-scale the error bars to yield a best-reduced χ² near unity. We find that errors about four times larger than the photon noise (for the binned data) are necessary to obtain a best-fit reduced χ² of unity. Error re-scaling factors almost as large as three have been required in other investigations, even for very precise spectrophotometry (Bean et al. 2011). An even larger error re-scaling in our case is probably a consequence of the complexity of the stellar oscillation spectrum that is only approximated using two oscillation periods. In other words, small amplitude residual oscillations may masquerade as noise.

We experimented with adding additional oscillation modes to the fit, and indeed this does reduce the necessity for error re-scaling. However, we can only justify initializing the fit using the two periods that we can objectively identify in our wavelet power spectra. Other modes, if present, cannot be resolved as clear oscillations in our ground-based data, instead they appear only as extra noise. We also verified that adding a third oscillation mode to our MCMC fit does not significantly alter the best-fit eclipse depth.

We vary nine parameters in the MCMC chains: an amplitude, period, and phase for two independent sinusoids, an eclipse depth and central phase, and an additive constant. We generate the eclipse curve using a new version of the Mandel & Agol (2002) methodology (see Deming et al. 2011b). We adopt system parameters that determine the shape and duration of the eclipse, from Collier-Cameron et al. (2010), except for the orbital period where we use the value from Smith et al. (2011). Observations of recent transits in the Czech Exoplanet Transit Database¹² show that the original period from Collier-Cameron et al. (2010) is too short; the slightly longer period derived by Smith et al. (2011) is more consistent with recent transit times in the Czech database.

The top panel of Figure 6 overplots the best-fit curve, including both the eclipse and oscillatory component. The lower panel removes the oscillatory component from the data and shows the comparison with the eclipse curve alone. Our best-fit eclipse has a depth of 0.0027 ± 0.0002, with central phase of 0.4995  ±  0.0010. Those errors are implied by our MCMC posterior distributions. They include possible degeneracies between fitted parameters, but there can be external sources of error than are not represented in the MCMC chains. The most obvious source of such errors is the presence of unmodeled stellar oscillation structure, as noted above. The most realistic check on the errors is to fit each of the two independent nights separately, and compare the independent results. Those fits yield eclipse depths of 0.0024 ± 0.0002 and 0.0033 ± 0.0004 for October 10 and 16, respectively. The best-fit central phases for the separate nights are 0.502 ± 0.001 and 0.495 ± 0.001, so the difference between two independent nights exceeds the formal errors, especially in the case of central phase. We adopt the best-fit values determined for the combined and binned data, but we assign the errors based on the differences in the two independent nights. The difference between two independent values drawn from a normal error distribution is, on average, twice the error associated with the average of those two values. So we estimate the errors appropriate to our average K_s-band eclipse depth and central phase to be half the difference between the best-fit values for the individual nights. Our best-fit eclipse depth and central phase is given in Table 3, together with the Spitzer eclipse results.

Table 3. WASP-33b Secondary Eclipse Depths, Time of Central Eclipse, and Orbital Phase

Wavelength	Eclipse Depth	Time as BJD(TDB)	Phase	Comment
(μm)
2.15	0.0027 ± 0.0004	2455844.8156 ± 0.0040	0.4995 ± 0.0035	Average of October 10 and 16; Figure 6
3.6	0.0026 ± 0.0005	2455647.1978 ± 0.0001	0.50041 ± 0.00008	Warm Spitzer; Figure 7
4.5	0.0041 ± 0.0002	2455650.8584 ± 0.0005	0.5012 ± 0.0004	Warm Spitzer; Figure 8

Notes. For the ground-based eclipse at 2.15 μm the barycentric time of central eclipse is based on the average phase for both October 10 and 16, then converted to the central BJD value for the October 10 eclipse.

Download table as:  ASCIITypeset image

4.4.2. The Eclipse in Warm Spitzer Bands

At both Spitzer wavelengths, we solve for the eclipse depth and central phase also using an MCMC method, with Gibbs sampling. We allow for an exponential baseline ramp (Harrington et al. 2007), because we find that a linear ramp is inadequate. Indeed, detailed analysis of very high S/N Spitzer 8 μm data (Agol et al. 2010) indicates that even a second exponential term can be warranted. However, we use a single exponential for two reasons. First, the single exponential ramp alone is more complex than is normally required for these specific Spitzer bands (Knutson et al. 2009; Todorov et al. 2012)—exponential ramps are normally only required at the longest wavelength Spitzer bands. Second, our data do not attain sufficient S/N to justify the inclusion of a second exponential term.

As in the case of our ground-based data, we calculate wavelet power spectra at both Spitzer bands, and these are shown as the upper panels of Figures 7 and 8. Like the ground-based data (perhaps a coincidence), we see two significant oscillatory spectral peaks after removing the best-fit eclipse. Figure 7 (3.6 μm) shows one of these peaks at 68 minutes. However, Figure 8 (4.5 μm) shows two peaks, 39 and 68 minutes. The former is likely to be associated with the Spitzer telescope pointing (Section 4.2), while the latter is clearly due to the stellar oscillation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our first MCMC fits at both Spitzer wavelengths vary 12 parameters: two sinusoidal terms (each having amplitude, period, and phase), an eclipse of variable depth and central phase (i.e., two parameters), the exponential ramp (three parameters), and a zero-point constant in intensity. We discard the first 20% of each 10⁶-step Markov chain and locate the best-fit values from the minimum in χ² (we keep track of χ² during the evolution of the chain). Subtracting the best fit from the data, we find that the residuals have a quasi-sinusoidal shape. These residuals represent the unmodeled portion of the data. They can be due to red noise in the data created by uncorrected instrumental effects, or by imperfect modeling of the stellar oscillations. Unlike the case of our ground-based data, the Spitzer data have sufficient stability and S/N to define these small amplitude imperfections in the fit. Both the oscillating intensity of the star and (to a lesser extent) the subtlety of red noise caused by Spitzer's improved performance require some new methodology in their analysis. We now introduce that new methodology.

Following the initial fit using a single MCMC chain of 10⁶ steps, we re-fit by including an additional term to explicitly account for the imperfections in the initial fit. We calculate residuals for the first MCMC fit by subtracting that best fit from the data. We then approximate those residuals using a wavelet decomposition, using N coefficients for Morelet wavelets, and we vary N in our subsequent analysis. For each choice of N, we multiply the wavelet decomposition of the residuals by an adjustable amplitude, and include that amplitude as a fit parameter in subsequent MCMC chains. In the limit where N equals the number of data points, this procedure would be trivial, because the wavelet decomposition would reproduce the residuals exactly, and including those residuals in the fit would simply be a fudge, wherein we arbitrarily subtract that part of the data not accounted for by the model. In the real analysis, this procedure is not trivial. With not-too-large N, the wavelet decomposition approximates only the major features of the residuals, not the point-to-point noise. Because the MCMC chain varies the amplitude of those major features, correlating the chained values of eclipse depth and phase with the amplitude of the unmodeled features indicates the degree to which the best-fit eclipse depth and phase depend on unmodeled aspects of the data. We determine the optimal N by examining the noise in the final fits, requiring that it be close to, but not less than the photon noise, and have no detectable red noise component. This dictates N = 9 for both our 3.6 and 4.5 μm data. Our final values are based on four independent MCMC chains at N = 9, each having 10⁶ steps, with the first 20% of each chain being ignored.

Using this procedure, we find that the Spitzer eclipse depths and central phases are reasonably robust in the sense that they do not vary greatly with N. This is especially true for the central phase at 3.6 μm, and the eclipse depth at 4.5 μm. In those cases, the variation in eclipse depth and central phase, as we vary N, are consistent with the errors implied by the MCMC posterior distributions. However, for the converse set of parameters—the 3.6 μm eclipse depth and 4.5 μm central phase—the agreement is not as good. In those cases, the differences in best-fit values as we vary N are larger than the errors from the MCMC posterior distributions, and we adopt the larger errors consistent with the fluctuations in the best-fit values as N varies.

The derived best-fit eclipse depths, central phases, and errors are included in Table 3. We also extract the amplitude of the stellar oscillation in the Spitzer bands (see Table 2), based on the total amplitude of the oscillatory portion of the MCMC fits, but with the 40 minute portion subtracted because we attribute that portion to the Spitzer observatory.

We interpret our observed planet–star flux contrast of WASP-33b using plane-parallel models of the dayside atmosphere of the planet. We use the atmospheric modeling and retrieval methodology of Madhusudhan & Seager (2009, 2010). The model computes line-by-line radiative transfer for a plane-parallel atmosphere with the assumptions of hydrostatic equilibrium and global energy balance. The composition and pressure–temperature (P-T) profile of the atmosphere are free parameters in the model. Since there are as yet insufficient data to constrain abundances in WASP-33b, we adopt solar composition for all elements except carbon, which we allow to increase in abundance, following Madhusudhan et al. (2011b). We explore a variety of temperature profiles, especially temperature profiles that are consistent with nearly complete absorption of stellar irradiance. The models include all the major opacity sources expected in hot hydrogen-dominated atmospheres, namely H₂O, CO, CH₄, CO₂, TiO, and VO, and collision-induced absorption (CIA) due to H₂–H₂. Our molecular line lists are obtained from Freedman et al. (2008), R. S. Freedman (2009, private communication), Rothman et al. (2005), Karkoschka & Tomasko (2010), and E. Karkoschka (2011, private communication). Our CIA opacities are obtained from Borysow et al. (1997) and Borysow (2002). A Kurucz model (Castelli & Kurucz 2004) is used for the stellar spectrum, and the stellar and planetary parameters are adopted from Collier-Cameron et al. (2010).

The first and most obvious result from our observations is the existence and prominence of the stellar oscillations. WASP-33 was already known to exhibit oscillations, but the amplitude observed by Herrero et al. (2011) in the optical (Johnson R band) was about 0.001. Our results (Figures 3–5 and Figures 7– 8) show oscillation amplitudes in agreement with the optical, except for K_s band where the amplitude is about twice the optical value (2.15 μm, Table 2), as noted in Section 4.3. Because the largest difference with the optical amplitude is seen in our ground-based (K_s-band) data, we contemplated whether the difference could be attributed to errors in our ground-based result. An argument against that possibility is the prominence of the stellar oscillation in the raw photometry (e.g., upper panel of Figure 2). We therefore explore whether properties of the stellar atmosphere that may be unique to K_s band could cause the oscillations to have greater amplitude at that wavelength.

We here consider the possibility that the larger K_s-band oscillation amplitude as compared with the optical is due to the different height of formation for continuum radiation in the stellar atmosphere, in concert with height-dependent variations in the mode amplitudes. Due to the increase in atomic hydrogen bound-free continuous opacity in the infrared, our K_s-band observations of WASP-33 sample a greater height in the stellar atmosphere than that observed by Herrero et al. (2011).

The upward propagation of a pressure-mode oscillation in a stellar atmosphere can in principle cause the mode amplitude to increase. As a propagating mode encounters lower mass density, the wave velocity—and hence the temperature perturbation in the compression—increases. However, propagation is strongly affected by the stratification of the stellar atmosphere (Marmolino & Severino 1991). Frequencies less than the acoustic cutoff frequency will not propagate, and their velocity amplitude decreases with height. To the extent that the temperature amplitude scales with velocity, it too will decrease with height for non-propagating modes. The acoustic cutoff frequency is c/2H, where c is sound speed and H is the pressure scale height. We calculated the acoustic cutoff frequency and other parameters for WASP-33, using a T_eff/log (g)/[M/H] = 7500/4.5/0.0 model atmosphere from Castelli & Kurucz (2004). We find that the acoustic cutoff corresponds to an oscillatory period of about 1 minute. The much longer period oscillations we observe for WASP-33 are therefore not propagating, and their amplitudes should decay with height.

The dominant opacity due to atomic hydrogen (Menzel & Pekeris 1935) is higher at 2.15 μm versus the optical by a factor of about 1.6, translating to a height difference of about 60 km. That is not sufficient to account for significant changes in the mode properties even for propagating modes and, as noted above, these modes do not propagate. Moreover, if height dependence of the stellar continuous spectrum were significant to our observed amplitudes, then we would expect even larger amplitudes in the Spitzer bands than at 2.15 μm, which are not observed.

However, there is one unique feature of the K_s band that may well be responsible for a higher oscillation mode amplitude. The Brackett-γ line at 2.165 μm is centered in the K_s bandpass. The strong opacity in that line for an A5V star could have a large potential effect on oscillation amplitudes. The impact of strong oscillations in the line, when diluted over the broad K_s band, could potentially be calculated using techniques beyond the scope of this paper (i.e., radiation hydrodynamics). A more direct method would be to obtain infrared spectroscopy of the star and directly measure the oscillations in the infrared hydrogen lines.

Our measured times of central eclipse can in principle determine ecos ω for the planet's orbit (e.g., Knutson et al. 2009). Considering the 25 s of light travel time across the planet's orbit, we expect to find the eclipse at phase 0.50024 if the orbit is circular. Weighting both our ground-based and Spitzer eclipses (Table 2) by the inverse of their formal variances, we find an average eclipse phase of 0.50044 ± 0.00008, totally dominated by the 3.6 μm eclipse. Thus, we find that ecos ω—approximated as π/2 times the phase offset from 0.5—is 0.0003 ± 0.00013. The central phase measurement for these eclipses is complicated by the stellar oscillations, so more-than-usual caution is needed in the interpretation of the measured central phase. Moreover, our value for ecos ω differs from zero by less than 3σ, so our results provide little evidence for a non-circular orbit.

Figure 9 shows our result for the eclipse of WASP-33b in comparison to two models of its atmosphere, both of which agree with the available measurements to date. One of these models has a temperature inversion with solar composition, and one has a non-inverted atmospheric structure with a carbon-rich composition (Madhusudhan et al. 2011b). Their temperature profiles are shown in the bottom panel of Figure 9, with the approximate formation depths of the four bandpasses overplotted as points. The K_s bandpass is relatively devoid of strong molecular absorption features, and probes the relatively deep planetary atmosphere (pressure, P ∼ 0.6 bars), relatively independent of the composition of the atmosphere (asterisks on Figure 9). (The 3.6 μm bandpass also peaks relatively deep in the atmosphere, but has significant contribution from higher altitudes, having more molecular absorption than the K_s band.) The large eclipse depth we observe in K_s band (brightness temperature ≈3400 K) thus indicates a hot atmosphere at depth, and a high effective temperature for the planet. Both models illustrated in Figure 9 have hot lower atmospheres (both above 2500 K) and both have inefficient longitudinal energy re-distribution. These models are consistent with the observed tendency for the most strongly irradiated planets to exhibit the least longitudinal re-distribution of heat (Cowan & Agol 2011; Perna et al. 2012). We find that these very hot models are necessary to reproduce our results as well as the result of Smith et al. (2011), and we conclude that WASP-33 strengthens the Cowan & Agol (2011) result.

Zoom In Zoom Out Reset image size

The two models we show are representative of a larger set of solutions that explain the data with and without thermal inversions. Given that there are 10 model parameters (Madhusudhan & Seager 2009, 2010) and only four data points, it is not possible to derive a unique model fit to the data. We ran large MCMC chains (of ∼10⁶ models) with and without thermal inversions, and identified regions of composition space in each case that are favored by the data (Madhusudhan & Seager 2010).

Both models on Figure 9 are unusual as compared with, for example, the well-observed archetype HD 189733b (Charbonneau et al. 2008). One model in Figure 9 adopts solar composition but with an inverted temperature structure (temperature rising with height), while the other model has temperature declining with height, but requires a carbon-rich composition. We integrate the fluxes of each planetary model, and the Kurucz model stellar atmosphere, over the observational bandpasses, and ratio those integrals. These band-integrated points are shown as squares in Figure 9. The χ² values for the models as compared to all four observed points (our three measurements, plus Smith et al. 2011) are 7.9 for the inverted solar-composition model and 2.8 for the non-inverted carbon-rich model. The K_s-band point at 2.15 μm favors the non-inverted model. Although the difference is not sufficiently significant to rule out the inverted model at this time, additional eclipse observations in the K_s band would be helpful for ruling out this type of atmospheric structure.

In our population of models with thermal inversions, several models with slightly different inverted temperature structure fit the data almost equally well. However, none of them fit the K-band point to within the 1σ errors while also fitting the remaining points. The Figure 9 model is the best among this set of inverted models.

In our models without thermal inversions, the best-fit model requires a carbon-rich composition (i.e., C/O ⩾ 1). However, at the 2σ level of significance per wavelength point, several solar composition models (not illustrated) provide an acceptable fit to the data. So, although the carbon-rich composition is favored, a solar abundance composition cannot be absolutely ruled out.

Our results illustrate the limitations of eclipse photometry in broad bands, especially for challenging cases like planets orbiting oscillating stars. Once we admit the possibility of non-solar compositions (because we are largely ignorant of true exoplanetary compositions), the range of models that can fit broadband photometry can be large, in this case extending to inverted and non-inverted models with drastically different temperature structure. The degeneracy is exacerbated by the relatively small range of atmospheric pressures probed by the four bandpasses we analyze (points on the bottom panel o Figure 9). Fortunately, future observations can break this degeneracy using Hubble Space Telescope/WFC3 spectroscopy near 1.4 μm (Berta et al. 2012). The water band near 1.4 μm is sufficiently strong that eclipse observations with the Hubble WFC3 grism should be feasible. Although water absorption in the carbon-rich model is suppressed by C/O > 1, the water emission in the solar-composition inverted model is predicted to be significant, readily detectable near 1.4 μm (see Figure 9). Moreover, spectroscopic techniques can potentially probe a larger depth range in exoplanetary atmospheres than does photometry because the cores of resolved spectral features have strong opacities. Our results therefore illustrate the complementary value of acquiring both broadband and spectroscopic observations of transiting exoplanets at secondary eclipse.