Thermal emission from WASP-33b, the hottest known planet

We report ground-based observations at 0.91 microns of the occultation of the hot Jupiter WASP-33b by its A5 host star. We measure the planet to be 0.109 +/- 0.030 per cent as bright as its host star at 0.91 microns. This corresponds to a brightness temperature, T_B = 3620 +200 -250 K, significantly higher than the zero-albedo equilibrium temperature for both isotropic re-radiation (2750 +/- 37 K) and uniform day-side only re-radiation (3271 +/- 44 K), but consistent with the zero-redistribution temperature (3515 +/- 47 K). This indicates that the heat redistribution from the day-side of WASP-33b to the night side is inefficient, and further suggest that there is immediate re-radiation, and therefore little or no redistribution, of heat within the day-side. We also detected the stellar pulsations of WASP-33, which we model as the sum of four sinusoids, with periods of between 42 and 77 minutes and amplitudes of 0.5 to 1.5 mmag.


INTRODUCTION
Extra-solar planets which transit their host stars are also likely to exhibit occultations (secondary eclipses), in which the planet passes out of sight behind the star. Observing an occultation allows measurement of the planet's emergent flux. This flux consists of thermal emission and reflected light, but is strongly dominated by thermal emission in the infrared wavelength regime. The first detections of thermal emission from an exoplanet were made using the Spitzer space telescope (Charbonneau et al. 2005;Deming et al. 2005), since when the occultations of many more planets have been observed at wavelengths of 3.6 µm and greater with Spitzer.
Complementary to Spitzer observations at 3.6, 4.5, 5.8, 8 and 24 µm, are ground-based observations made in the near infrared. Such observations are important because they extend the observed planetary spectral energy distribution (SED) to lower wavelengths, and to the peak of the SED, ⋆ Based on service observations made with the William Herschel Telescope operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofsica de Canarias. † E-mail: amss@astro.keele.ac.uk which is around 1 µm for the shortest-period, hottest planets.
Measurement of a planet's occultation depth allows the brightness temperature to be determined, which leads to an evaluation of how efficiently heat is redistributed from the day-side to the night-side of the planet.
A planet's observed SED (from occultation observations at several different wavelengths) allows the chemical composition of the atmosphere to be inferred (e.g. Swain et al. 2010). Detection of a thermal inversion in the atmosphere is also possible; observations to date indicate two classes of hot-Jupiter atmosphere, those with and those without such an inversion. One possible explanation for the existence of these two classes is that provided by Knutson, Howard & Isaacson (2010) who argue that chromospherically active stars host planets with no temperature inversion.
Occultation observations also allow refinement of the orbital parameters of a planetary system, particularly the orbital eccentricity, which is related to both the timing of the occultation and its duration.
WASP-33 (HD 15082) was first identified as a transiting planet candidate by Christian et al. (2006), but was identified as fast rotator, which rules out the use of precision radial velocity measurements to confirm its plane-tary nature and to help characterise the system. Instead, Collier , hereafter C10, used lineprofile tomography during transit to confirm WASP-33b as the first known planet to orbit an A-type star.
Given the nature of its host star and its very short orbital period (1.22 d), WASP-33b has the largest equilibrium temperature (T eql ≈ 2700 K, assuming zero-albedo and uniform heat redistribution; C10) of any known exoplanet. WASP-33b is an excellent target for ground-based occultation observations, due to the large predicted planet-to-star flux ratio. Planets with similarly large brightness temperatures include WASP-19b (T eql ≈ 2000 K;Hebb et al. 2010), occultations of which have been observed from the ground in H-band and K-band (Anderson et al. 2010a;Gibson et al. 2010) and WASP-12b (T eql ≈ 2500 K; Hebb et al. 2009), observed from the ground in the KS-, H-and J-bands by Croll et al. (2011).
C10 also reported that the stellar line profiles show evidence for non-radial pulsations, and suggest that the star may be a γ Dor-type variable. Further, photometric, evidence for these pulsations is presented in Herrero et al. (2011), hereafter H11, who suggest that WASP-33 is a δ Scuti-type variable.
In this paper, we present ground-based observations of the occultation of WASP-33b at 0.9µm.

OBSERVATIONS AND DATA REDUCTION
We observed the occultation of WASP-33b on 2010 October 29/30 using the Auxiliary Camera (ACAM) of the William Herschel Telescope (WHT). The ACAM CCD covers a fieldof-view of 8 ′ ×8 ′ with 2048×2048 pixels, giving a plate scale of 0.25 ′′ per pixel. The target was positioned on the CCD such that three relatively bright comparison stars are present in the images. A total of 203 images were taken from 22:43 to 05:55 UT.
The observations were conducted using a narrow-band S[III] filter (ING filter #S9077), which has a central wavelength of 9077Å, and a FWHM of 54Å. The telescope was defocussed to spread the light from the target over a large number of pixels, thus reducing noise associated with imperfect flat-fielding and allowing longer integration times. Exposure times of 100, 120 and 150 s were used. The conditions were not photometric; thin cirrus cloud affected the observations for the duration of the night.
Aperture photometry was performed on the flat-field and bias corrected images, using the photometry pipeline of Southworth et al. (2009), which is based around the astrolib / aper routine 1 . The pointing was monitored by cross-correlating each image with a reference image, and the apertures shifted to follow the stars.
Photometry was performed using a range of different aperture radii, from 8 to 65 pixels. A radius of 50 pixels was chosen to maximise the target signal-to noise, and ensure that a large fraction (> 95 per cent) of the target flux is included in the aperture. The flux from the three comparison stars was combined to form an ensemble reference star, 1 The astrolib subroutine library is distributed by NASA. For further details see http://idlastro.gsfc.nasa.gov approximately 0.5 magnitudes fainter than WASP-33. The resulting differential light curve is shown in Fig. 1.

DATA ANALYSIS
In addition to an apparent eclipse-like signal (centred about JD ≈ 2455499.58 as expected, assuming a circular orbit), the raw light curve ( Fig. 1) exhibits large systematic effects, specifically (i) a series of sine-like oscillations, which we attribute to the stellar pulsations for which evidence was presented by C10 and H11; and (ii) an increase in flux across the second half of the light curve. We analyse the new occultation data alongside the SuperWASP, Keele and JGT light curves and the radial velocities of C10, and the Montsec Observatory light curve of H11. This global analysis is performed using the Markov Chain Monte Carlo (MCMC) code described in Collier Cameron et al. (2007), Enoch et al. (2010) and Anderson et al. (2011b).
In addition to fitting the occultation, we also simultaneously fit functions describing the stellar pulsations and systematic effects present in the raw light curve (Fig. 1). The de-trending within the MCMC code is done using the svdfit routine (Press et al. 1992). Initially, we neglected the stellar pulsations, and attempted only to remove the instrumental effects. We tried de-trending with several parameters: airmass, time, sky background amplitude, and the position of the target on the CCD; we also tried de-trending using various combinations of parameters and various functional forms for the fit to each parameter. We used the Bayesian information criterion (BIC) to discriminate between the different models. In determining the best de-trending model we did not re-scale the error bars on any of the data. The lowest BIC value was produced when the data were de-trended with a quadratic function of the sky background, where f sky is the sky background flux and a0, a1, and a2 are the coefficients we fit for.
The light curve resulting from our de-trending with sky background, while an improvement upon the raw light curve, still exhibits the sine-like variations present in the raw light curve which we attribute to stellar pulsations. It is apparent the signature of these pulsations is multi-periodic; a model consisting of a single sine curve produced a very poor fit to the data. We conducted a Fourier analysis on the out-ofoccultation data using the Period04 time-series analysis software package (Lenz & Breger 2005). This reveals three frequencies which are deemed significant (with a signal-tonoise ratio greater than 4). These frequencies correspond to periods of 53.99±0.29 min, 76.90±0.62 min and 41.85±0.31 min (see periodogram, Fig. 2). The 68-min period of H11 is present in our periodogram, but with a slightly lower amplitude than the 54-min period. When prewhitening with the 54-min period, the the significance of the 68-min period is much reduced, suggesting that one of these two periods may be an alias of the other. When, instead we used the H11 period and amplitude as the starting point for a Period04 analysis yields a different set of frequencies, but a poorer overall fit to our data. We also tried a period search of our data combined with the existing photometry, but the results from this were inconclusive.
We then fitted for either a single pulsation period (that of H11), three pulsation periods (from our Period04 analysis) or four pulsations periods (our periods and that of H11) with an MCMC code (see below). Comparing the BIC values of the resulting models we found that the 4period model (∆BIC = 0) is preferred to both the 3-period (∆BIC = 391) and the 1-period (∆BIC = 1057) models.
The four periods and uncertainties were input to a version of our MCMC code, modified to include the four pulsation periods as proposal parameters. The periods were, therefore, allowed to vary but Gaussian priors were applied to prevent the periods from drifting too far from the fitted periods. The priors were applied by means of a Bayesian penalty added to our merit function (χ 2 ), given by where Pi,0 and σ 2 P i,0 are the initial value and uncertainty of a particular period, Pi, respectively. The amplitudes and phases of each of these sine waves were fitted using svdfit, along with the coefficients of Equation 1. The function fitted is of the form where g(f sky ) is Equation 1, t is time, and bi and ci are the coefficients we determine in the fit. This is equivalent to f sky is the sky background flux, z is the zenith distance of the target, and δt is the time since the start of the observations.
where Ai are the amplitudes, and δi the phase offsets of each sine term. We then again tried different forms of Equation 1 to detrend the data, again finding that a quadratic in sky background results in the lowest BIC value (Table 1). Once we had selected the functional form of our trend model, we performed a further MCMC analysis, this time scaling the error bars of each photometric dataset so as to obtain a reduced χ 2 of unity.
We initially fitted for the orbital eccentricity and the resulting solution is significantly eccentric (e = 0.174 +0.058 −0.039 ), but e cos ω is close to, and consistent with, zero (e cos ω = 0.0027 +0.0079 −0.0033 ). The argument of periastron, ω, in this solution is −89.0 +3.0 −1.2 • , i.e. the major axis of the orbit is aligned almost exactly with our line-of-sight. We suggest that this is an improbable configuration, and is caused by a degeneracy in fitting the occultation and the stellar pulsations. e cos ω is constrained by the timing of the occultation, whereas the duration of the occultation informs us about e sin ω.
Furthermore, in the case of WASP-33, the shape of the radial velocity curve places almost no constraint on the eccentricity (because the RVs are imprecise due to the nature of the star). This allows the fitted occultation duration to deviate from the transit duration to a large extent. The light curve, after de-trending for sky background and stellar pulsations still clearly contains correlated ('red') noise, which is likely to be caused by pulsation activity at further periods (Fig. 4). It is this red noise that mimics the ingress and egress of the occultation and results in an occultation duration significantly shorter than the transit duration, and an argument of periastron close to −90 • .
The eccentric solution requires that the planet be moving more slowly during transit than for the circular solution, and this in turn requires that the stellar radius is smaller in the eccentric case. This results in a stellar density (ρ * = 0.75 +0.16 −0.12 ρ⊙) which is greater than that of a mainsequence star with the mass of WASP-33 (ρ * ≈ 0.5 ρ⊙; Gray 2008). The stellar density of our circular solution, however, lies very close to that predictionρ * = 0.44 ± 0.05 ρ⊙.
Although short-period planets have been found in significantly eccentric orbits (e.g. WASP-14b, Joshi et al. 2009;Husnoo et al. 2010), the orbit of WASP-33b is expected to have been circularised by tidal interactions with the star. These are anticipated to be particularly strong because of the relatively large radii of planet and star, and the small orbital separation.
Using Equation (1) of Jackson, Greenberg & Barnes  Table 2, and MP = 4.1MJup (the upper limit established by C10, which will maximise the calculated circularisation time-scale). We obtain where QP and Q * are the tidal dissipation parameters for the planet and star respectively. Adopting QP = 10 5.5 and Q * = 10 6.5 (the best-fitting values from the study of Jackson, Greenberg & Barnes 2008), we find τe = 1.72 Myr. This time-scale is an order of magnitude less than the 25 Myr main-sequence age estimate of C10, and a very small fraction of the 500 Myr upper limit to the age. Although this suggests that there has been ample time for the orbit to have been circularised, τe is affected by large uncertainties in the Q parameters, so we cannot say for certain that this is the case. C10 determined that the orbital axis of WASP-33b is misaligned with respect to the stellar spin axis, such that the planet orbits in a retrograde sense. This suggests that the system may not exist in a dynamically relaxed state, and the orbit may be eccentric. The time-scale for co-planarisation of the orbit, however, is very much longer than that for its circularisation. Using Equation (14) of Winn et al. (2005) we obtain the characteristic time-scale, τ ψ , for significant change in the angle ψ between the planetary orbital angular momentum vector and the stellar rotation angular momentum vector, τ ψ = 11.4 Q * 10 6.5 Gyr.
Here, we have assumed values of 0.01 for the apsidal motion constant, √ 0.1 for the dimensionless radius of gyration, and a planet mass of 4.1 MJup. The time-scale increases with decreasing planet mass. It is therefore perfectly reasonable to suppose that the orbit of WASP-33b is circular, despite the fact that it is retrograde.
Furthermore, the empirical evidence provided by other planetary systems lends support to the circular hypothesis. Pont et al. (2011) present an analysis of the eccentricities of a large number of hot Jupiters, determining that the only known planets with significantly detected, large eccentricities are massive and have relatively long-period orbits. Several planets with large spin-orbit misalignments are known to have near-circular orbits, such as WASP-17b which has a retrograde orbit (Anderson et al. 2010b), and an orbital  eccentricity close to zero, e = 0.028 +0.015 −0.018 (Anderson et al. 2011c).
All of the above evidence leads us to rule out the eccentric solution, and to take the approach (advocated by e.g. Anderson et al. 2011a) of adopting a circular solution. The best-fitting parameters for this circular model are presented in Table 2. The fit to the pulsations is shown in Fig. 3 and the fit to the occultation is shown in Fig. 4.
We assessed the presence of correlated noise in the occultation light curve residuals, by plotting the RMS of the binned residuals to the best-fitting model as a function of bin width (Fig. 5). The RMS of the binned residuals deviates significantly from the uncorrelated noise expectation, indicating significant levels of red noise are present in the light curve. This is likely to be due to the presence of pulsa- tion signals in addition to those we fitted. The results of a Fourier analysis to determine whether any further periodicity is present in the data is inconclusive, in part due to the limited baseline of the data. To determine whether fitting an occultation model to the data is justified, we performed an F-test to compare our best-fitting occultation model with a constant-flux model. The resulting p = 0.013 (i.e. there is only a ∼ 1% chance of our occultation fit arising if there were no underlying occultation) indicates that the occultation model is preferred at the > 95 per cent level.
The remaining correlated noise present in the occultation light curve means that the uncertainty on the occultation depth is likely to have been underestimated. To calculate a realistic uncertainty, we employ the residual permutation or 'prayer bead' method (e.g. Gillon et al. 2007). We took the light curve residuals and shifted each residual to the subsequent observation before adding back the trend, pulsation and occultation models. This process was repeated, shifting each residual two observations, and so on, resulting in a total of n obs − 1 = 202 light curves, each of which have the time-structure of the correlated noise preserved. A separate MCMC analysis was performed on each of these light curves.
We then calculated the median of the best-fitting occultation depths and pulsation periods from the new fits, as well as the uncertainty limits which enclose 68 per cent of the values around the median. The uncertainties on these parameters quoted in Table 2 are those calculated from the permutation of the residuals. In all cases, the median values from the residual permutation analysis are consistent with, and close to the best-fitting values from our single MCMC analysis. As expected, the uncertainties on these parameters derived from the permutation analysis are larger. The original occultation depth we found was 0.00108±0.016, whereas the value from the residual permutation is 0.00109 ± 0.030.

DISCUSSION
To calculate the brightness temperature of the planet, we modelled the star using a synthetic spectrum of an A5 star Pickles (1998) normalised to the integrated flux of a 7430 K (Table 2) black-body emitter, and the planet as a blackbody of temperature TB. We then defined the measured occultation depth as the product of the ratio of the bandpassintegrated planet and star fluxes and the planet-to-star area ratio (e.g. Charbonneau et al. 2005). We calculate a brightness temperature, TB = 3620 +200 −250 K. The uncertainty on TB only accounts for the uncertainty on the measured occultation depth, which is the dominant source of error. Additional, smaller sources of uncertainty on the ratio of planetary and stellar areas, and the stellar effective temperature are not accounted for.
The equilibrium temperature for a zero-albedo planet is given by TP,A=0 = f 1 4 T * ,eff R * 2a , where f = 1 indicates isotropic re-radiation over the whole planet (i.e. the redistribution of heat from the day-side to the nightside is fully efficient). The case where the heat is uniformly distributed on the day-side of the planet, but there is no redistribution of heat to the night-side corresponds to f = 2. A third case, corresponding to f = 8 3 (equivalent to ε = 0 and f = 2 3 in the notations of Cowan & Agol 2011 and López-Morales & Seager 2007 respectively), occurs when the incident radiation is immediately re-radiated and so there is no redistribution of heat even within the day-side. Because the hottest region of the day-side is most visible close to occultation, a deeper occultation is observed than if the day-side had a uniform temperature. The brightness temperature of 3620 K is greater than either the equilibrium temperature for uniform or no redistribution to the night-side (respectively 2750 K and 3271 K, Table 2). The brightness temperature is consistent with the zero-redistribution (f = 8 3 ) case (3515 K), however. This suggests that the heat-transport from the day-side to the night-side is inefficient, and further that the day-side is unlikely to have a uniform temperature. This is in line with observations of other highly irradiated planets, which are also found to have poor heat redistribution efficiencies (e.g. Cowan & Agol 2011).

CONCLUSIONS
We have detected thermal emission from WASP-33b. We measure the occultation depth at 0.91µm to be 0.109±0.030 per cent, which corresponds to a brightness temperature of 3620 +200 −250 K, the hottest such temperature recorded for an exoplanet.
We also detect the non-radial pulsations of the host star in the photometry; and conclude that a multi-periodic solution is required to fit the pulsation signal. Our best-fitting model of the pulsation signal consists of sine terms with periods of 53.5, 76.3, 41.9, and 66.6 min, with amplitudes of 1.5, 0.6, 0.8 and 0.6 mmag respectively. None of these periods corresponds directly to the period of 68 min fitted by H11, but are similar, and our greatest amplitude agrees with their 0.9 mmag. We do however find evidence for the H11 68-min period in our periodogram, and it appears that the 54-min and 68-min periods may be aliases of each other. We do not claim that these periods are definitive, but rather they are the best fit to our data. Our analysis suggests that the pulsations of WASP-33 are complex and multi-periodic in nature. This multi-periodic, ∼ 1 h, pulsation signature lends support to the conclusion of H11 that WASP-33 is a δ Scuti-type variable. Our ability to determine the true periodicity of the pulsations is limited by the relatively short baseline of our data.
Although fitting for the orbital eccentricity returns a significantly non-zero value of e, we argue that because e cos ω is essentially zero, the eccentric solution is improbable. Furthermore, the eccentric solution is incompatible with the stellar analysis of C10, so we therefore adopt a circular model for the orbit.
High signal-to-noise photometry with a longer baseline is required to study better the pulsations, to determine the complete nature of the periodicity and to investigate whether, as H11 suggest, there is evidence for star-planet interactions in the WASP-33 system. Such a characterisation of the pulsation periods would also allow the pulsation signal to be more cleanly subtracted from occultation photometry. This, in turn, would allow a more precise measurement of the occultation depth and the unresolved questions concerning the orbital eccentricity to be answered. Measurements of the occultation depth at different wavelengths to the one presented here are required in order to construct a spectral energy distribution and hence to characterise the planetary atmosphere.