THE 4.5 μm FULL-ORBIT PHASE CURVE OF THE HOT JUPITER HD 209458b

The sky background is estimated in each image by first masking out the target star with a 10 pixel radius. 3σ outliers are then trimmed over three iterations to remove hot pixels and cosmic ray hits from the set of 64 measurements at that pixel position from a given data cube. The filtered background counts are fit with a Gaussian function. This count is then subtracted from each image.

We determine the x- and y-position of HD 209458 in each image using a flux-weighting methodology similar to a centroid that employs a five-pixel-radius circular aperture (e.g., Knutson et al. 2008). We find the latter method results in a more stable position estimate than the center of a two-dimensional Gaussian fit, and a correspondingly smaller scatter in our final photometry, in good agreement with our previous results for Spitzer 4.5 μm phase curve observations of HD 189733b (Knutson et al. 2012) and HAT-P-2 (Lewis et al. 2013). We first calculate our photometry using a fixed circular aperture with radii equal to 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.5, 4.0, 4.5, and 5.0 pixels.

We also explore apertures that scale according to the noise pixel parameter $\tilde{\beta }$, which is defined in Section 2.2.2 (IRAC Image Quality) of the IRAC instrument handbook as

$$\begin{equation} \tilde{\beta } = \frac{(\sum I_{i})^{2}}{\sum \left(I_{i}^{2}\right)} \end{equation} \tag{ 1 }$$

in which I_i is the intensity in a given pixel i. The noise pixel parameter is equal to one over the sharpness parameter S₁, (Muller & Buffington 1974), and proportional to the full width at half-maximum of the stellar point-spread function squared (Mighell 2005); for a more thorough discussion of the noise pixel parameter, please see Appendix A of Lewis et al. (2013). We generate two additional aperture populations using the noise pixel parameter: $\sqrt{{{{\tilde{\beta}}}}} \times$[0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 1.0, 1.05, 1.1, 1.15, 1.2] and $\sqrt{{{{\tilde{\beta }}}}} + [-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,$ 0, 0.1, 0.2, 0.3, 0.4]. We find that a fixed aperture size of 2.7 pixels produces the smallest amount of scatter in the final residuals, which is also consistent with our conclusions for previous 4.5 μm phase curve observations (Knutson et al. 2012; Lewis et al. 2013).

Outliers in the measured flux, x- and y-positions of the target, or the noise pixel parameter $\tilde{\beta }$ are removed by first discarding 6σ outliers and then by using a moving median filter with a width of 10 points and discarding outliers greater than 4σ. After applying this filter, we found that there were two points in our observations (at BJD_UTC ≈2455214.9 and ≈2455217.29) in which the star displayed a sudden, sharp excursion in position typically lasting several minutes. These excursions are most likely the result of micrometoerite hits and we exclude them from our subsequent analysis. The two excluded segments of the light curve contain a total of 984 images, corresponding to 0.13% of the total data set. This raw, filtered data is presented in Table 1.

Previous Spitzer/IRAC observations (e.g., Beerer et al. 2011; Todorov et al. 2012; Knutson et al. 2012; Lewis et al. 2013) note a ramp-like change in the observed flux at 4.5 μm, lasting approximately one hour. Agol et al. (2010) examined a similar effect in the 8 μm bandpass and postulated that the effect might be due to a combination of thermal settling of the telescope and charge-trapping in the array. This effect is mitigated by trimming the first hour of data after the start of the observation and again after the downlink break (at phase ≈0.75). We then fit the trimmed data with the ramp function given by Agol et al. (2010), as part of our global fit described in Section 2.3:

$$\begin{equation} F^{\prime }/F = 1 + a_{1}e^{-t/a_{2}} + a_{3}e^{-t/a_{4}} , \end{equation} \tag{ 2 }$$

in which F' is the flux with the ramp, F is the flux corrected for the ramp, a₁–a₄ are the correction coefficients, and t is the time. A Bayseian Information Criterion (BIC¹²; Schwarz 1978) analysis indicates that a ramp correction is unnecessary in order to fit the trimmed data. We therefore leave this ramp function out of our subsequent analysis, consistent with our previous 4.5 μm phase curve analyses (Knutson et al. 2012; Lewis et al. 2013).

Previous studies of transiting exoplanets find that the largest error source in Spitzer/IRAC 3.6 and 4.5 μm data is due to intra-pixel sensitivity variations, which cause variations in the telescope pointing to manifest as changes in the measured flux from the target star (e.g., Charbonneau et al. 2008; Beaulieu et al. 2010; Ballard et al. 2010; Knutson et al. 2012; Lewis et al. 2013). This effect can be easily seen in our 4.5 μm full-orbit data (Figure 1) in which the raw flux tracks with the x- and y-positions on the detector. We de-trend the data using a pixel-mapping Gaussian weight decorrelation method (Ballard et al. 2010; Knutson et al. 2012; Lewis et al. 2013) in which each photometric measurement is corrected by its 50 nearest neighbors in x-, y-, and $\tilde{\beta }$-space (see Equation (2) in Knutson et al. 2012). We find, consistent with Lewis et al. (2013), that inclusion of the noise pixel parameter $\tilde{\beta }$ terms increases the scatter in our final solutions. We therefore use the following simplified version of this correction instead:

$$\begin{equation} F_{{\rm meas},j} = F_{0,j} \sum _{i=0}^{n}e^{-(x_{i}-x_{j})^{2}/2\sigma _{x,j}^{2}} \times e^{-(y_{i}-y_{j})^{2}/2\sigma _{y,j}^{2}}, \end{equation} \tag{ 3 }$$

in which F_{meas, j} is the measured flux in the j^th image; F_{0, j} is the intrinsic flux; x_j and y_j are the measured x-position and y-position; and σ_{x, j} and σ_{y, j} are the standard deviations of the x and y vectors over the full range in i (0 to n, where n = 50 nearest neighbors). The 50 nearest neighbors to the j^th image are determined via

$$\begin{equation} d_{i,j} = \sqrt{(x_{i}-x_{j})^{2} + (y_{i}-y_{j})^{2}} , \end{equation} \tag{ 4 }$$

in which d_{i, j} is the distance vector between the j^th image and the i^th image; x_i and x_j are the measured x-positions of the i^th and j^th images, respectively; and y_i and y_j are the measured y-positions of the i^th and j^th images, respectively. As a first-order estimation, we set the uncertainties on individual flux values equal to the standard deviation of the normalized residuals (SDNR) of the raw flux divided by the Gaussian weight (see Equation (3)).

Zoom In Zoom Out Reset image size

We employ a Levenburg–Marquardt least-squares minimization via the Interactive Data Language mpfit function (Markwardt 2009) to solve simultaneously for the best-fit transit, secondary eclipse, and phase curve functions while optimizing the corresponding pixel map at each step in our minimization. We first use a Mandel & Agol (2002) model light curve modified for full-orbit observations to calculate primary transit and secondary eclipse light curves using an initial guess for the system parameters, then divide those out of the raw light curve before making an initial pixel map. We next calculate a pixel map and corresponding correction (Equation (3)) for the intrapixel sensitivity variations at each point in the light curve. Thus, the Levenburg–Marquardt performs an iterative fit in which it finds the light curve model that, when the data is divided by the model, gives the best x- and y-position decorrelation (or pixel map) which minimizes the scatter in the corrected data/residuals.

For this analysis, the orbital period (P = 3.5247455 days) and eccentricity (e = 0) are fixed to values from a previous multi-transit study (Torres et al. 2008) while the stellar limb darkening coefficients (ld_{1 − 4} = [0.4614, −0.4277, 0.3362, −0.1074]) are derived from a three-dimensional (3D) model stellar atmosphere (Hayek et al. 2012). The inclination i, ratio of semi-major axis to stellar radius a/R_s, mid-transit time T_c, ratio of planetary to stellar radius R_p/R_s, secondary mid-eclipse times T_c, and secondary eclipse depths are left as free parameters. We also search for longitudinal brightness variations across the face of HD 209458b by fitting the full-orbit light curve with the Cowan & Agol (2008) phase curve function:

$$\begin{eqnarray} F &=& 1 + c_{0} + c_{1} \cos (2 \pi t/P) + c_{2} \sin (2 \pi t/P) \nonumber \\ &&+\,c_{3} \cos (4 \pi t/P) + c_{4} \sin (4 \pi t/P) , \end{eqnarray} \tag{ 5 }$$

in which c₀ is the secondary emission of the planet (here, we use the deeper measured secondary eclipse depth to establish this baseline so that it has a relative flux of unity), c₁ − c₄ are free parameters, t is the time, and P is the planetary orbital period. The BIC indicates that inclusion of this phase curve function is necessary and that the c₃ and c₄ terms do not improve the quality of the fit; therefore, only the c₁ and c₂ terms are used in the final fit. Fitting simultaneously for the phase curve, transit, and secondary eclipse models allows us to accurately account for the effect that the phase curve shape has on our estimates of the various transit and secondary eclipse parameters (Kipping & Tinetti 2010). After finding an initial full-orbit solution to the data set, the photometric uncertainties are then inflated by a factor of 2.2 in order to produce a global best-fit solution with a reduced χ² equal to unity.

Next we explore the global solution-space with a Markov Chain Monte Carlo analysis (MCMC; e.g., Ford 2005). The initial best-fit parameter solutions are randomly perturbed to seed six independent MCMC chains each consisting of 10⁵ links (steps). In each MCMC chain, a new trial solution is drawn from a Gaussian-distributed parameter space based upon the Levenburg–Marquardt best-fit parameters and associated uncertainties. Multiple (six) sufficiently long (10⁵ links) MCMC chains are run with different randomly selected starting values to ensure both that an adequate amount of parameter space is sampled and that all chains converge to the global, and not a local, solution. After running each MCMC chain, we search for the point in which the chain has become well-mixed, as defined by where the χ² value first falls below the median of all the χ² values in the chain, and discard all links up to that "burn-in" point. A Gelman–Rubin test (Gelman & Rubin 1992) indicates a potential scale reduction factor ⩽1.07 for all parameters (Ford 2005), suggesting that all six MCMC chains have converged to the same global solution.

Because MCMCs assume that the noise in the data is Gaussian and uncorrelated from one measurement to the next, they will typically underestimate the true uncertainties for data with a significant component of time-correlated noise (Carter & Winn 2009). To account for this "red" correlated noise, we independently estimate the uncertainties on our fit parameters using the "residual permutation" or "prayer bead" method (e.g., Jenkins et al. 2002; Southworth 2008; Bean et al. 2008; Winn et al. 2008), in which the residuals from the best-fit MCMC solution are circularly permutated and added back onto the best-fit MCMC model to generate a new data set. A Levenburg–Marquardt fit is applied to each of these new "simulated" data sets in order to build up histograms of the best parameter distributions from each permutation. We then compare the resulting uncertainties to their counterparts from the MCMC analysis and take the larger of the two for each fit parameter as our final uncertainties. The prayer bead uncertainties were up to 2.6 times larger than the corresponding MCMC uncertainties. Our final fit parameters and their associated 1σ uncertainties are reported in Table 2 and the final, corrected photometry is presented in Table 1. The brightness temperatures calculated from the final fit parameters are listed in Table 3.

Despite the success of this reduction, some systematics still remain. For example, small oscillations in flux occur before and during the first secondary eclipse, and again near an orbital phase of 0.2. These features are not correlated with any change in x, y, or $\tilde{\beta }$, suggesting that they are not associated with any changes in Spitzer's pointing. In addition, they do not correspond with any background flux changes. These features could be due to stellar brightness variations from HD 209458, as their 0.1% amplitude is similar to the amplitude of flux variations measured for the Sun (Eddy 2009). However, the shape of the observed flux oscillations does not appear to match the sharp rise and slow decline expected for stellar flares (e.g., Figure 2 in Gary et al. 2012), and the timescale is too short for rotational spot modulations. Based on HD 209458's Ca ii H and K emission measurements (log$(R^{\prime }_{HK})=-4.97$; Knutson et al. 2010), HD 209458 appears to have an activity level comparable to that of our own Sun (log$(R^{\prime }_{{HK}})=-4.96$; Noyes et al. 1984); therefore, it might be reasonable to expect similar levels of infrared variability. Knutson et al. (2012) observed similar features in their Spitzer/IRAC HD 189733b data just after the first secondary eclipse at 3.6 μm and possibly before and after the second secondary eclipse at 4.5 μm. This similarity suggests that the features in our light curve might instead be the result of residual uncorrected Spitzer/IRAC systematics. We leave a full exploration of these features for a future study and note that they do not significantly affect the overall shape of the phase curve, the primary transit, or the second secondary eclipse.

The revision of the previous 4.5 μm emission measurement by Knutson et al. (2008) is significant because the high brightness temperatures at 4.5 and 5.8 μm were interpreted as evidence for a thermal inversion (Burrows et al. 2007; Knutson et al. 2008; Madhusudhan & Seager 2009, 2010; Line et al. 2014). Our revision of HD 209458b's 4.5 μm emission is not without precedent. The original studies of both HD 189733b (Charbonneau et al. 2008) and HD 209458b (Knutson et al. 2008) were carried out by continuously cycling between the four IRAC arrays in order to cover all four bands during a single secondary eclipse. This method reduced the effective cadence of these observations to one-fifth that of the now standard staring mode observations, in which the star is observed continuously in a single band. The decision to cycle between detectors resulted in a significantly higher level of pointing jitter and correspondingly large flux variations in the resulting light curves as compared to staring mode observations. The disadvantages of this mode became apparent soon after these data were taken, and no other stars were observed in this manner. Although the original data were reduced using the standard techniques available at the time (namely, fitting with a polynomial function of x- and y-position), it would not be surprising if the large pointing jitter resulted in a biased estimate of the eclipse depth. Subsequent observations of HD 189733b by Knutson et al. (2012) used the now standard staring mode and pixel mapping techniques, and found a 3.6 μm eclipse depth that deviated by approximately 7.5σ from the previous value. It is therefore not surprising that our new staring mode observations of HD 209458b also result in a 4.4σ revision to the published eclipse depth.

Our revised estimate for the 4.5 μm brightness temperature agrees significantly better with published models for this planet (Figure 5). Despite the original interpretation of the Knutson et al. (2008) data as resulting from a thermal inversion, it is worth emphasizing that even one-dimensional radiative-equilibrium models with several tunable free parameters (Burrows et al. 2007; Fortney et al. 2008) had difficulty simultaneously explaining all four IRAC observations. Models in which the relative abundances and pressure-temperature profiles were allowed to vary as free parameters were able to match all four of the original measurements (Madhusudhan & Seager 2009, 2010; Line et al. 2014). These "data-driven" best-fit models find the temperature inversion placed where they have the most leverage over the 4.5 and 5.8 μm bandpasses. GCM simulations including a hot stratosphere generally do not produce elevated brightness temperatures at 4.5 and 5.8 μm relative to the other IRAC bandpasses Showman et al. (2009), since the pressures probed by the IRAC bandpasses largely overlap, and the depth of the inversion is not a free parameter, but is controlled by the optical and infrared opacity sources. As shown in Showman et al. (2009), in models where pressure–temperature profiles are not left as a free variable, much of the flux in the IRAC bandpasses is emitted to space from the bottom edge of, or even below, the inversion, implying that high stratospheric temperatures do not necessarily exert significant leverage on the flux emitted at 4.5 and 5.8 μm in these models. Furthermore, the inversion covers only part of the dayside in the 3D models, which also mutes its effect. However, there could be an absorber not accounted for in the models that gives rise to an inversion at a different level.

Our new analysis is in better agreement with the predictions of the general circulation models for this planet from Showman et al. (2009), including cases both with and without dayside temperature inversions, and also blackbody models (see Figure 5). Although our new measurement does not rule out a dayside temperature inversion, it suggests that an inversion may not be necessary to explain the current 4.5 μm broadband data. New high-precision eclipse measurements in the other Spitzer bandpasses, specifically photometry at 5.8 and 8 μm, and ideally mid-infrared spectroscopic measurements (e.g., Madhusudhan & Seager 2010; Burrows 2014; Line et al. 2014) would help to provide a more definitive answer to this question.

Our derived 4.5 μm light curve indicates an eastward shifted hot spot, similar to the one observed for HD 189733b (Knutson et al. 2007, 2009b, 2012). Such eastward offsets were first predicted for hot Jupiters by Showman & Guillot (2002), and subsequent 3D circulation models confirm that they are a robust feature of the hot-Jupiter circulation regime, at least over a certain range of incident stellar fluxes and atmospheric radiative time constants (e.g., Cooper & Showman 2005; Showman et al. 2008a, 2009, 2013; Menou & Rauscher 2009; Rauscher & Menou 2010, 2012; Heng et al. 2011b, 2011a; Perna et al. 2012; Dobbs-Dixon et al. 2012). In these circulation models, the eastward offset of the hot spot results from advection of the temperature field by a fast, broad eastward equatorial jet stream—so-called equatorial superrotation. In turn, the equatorial superrotation results from interactions with the mean flow of standing, planetary-scale waves induced by the day-night heating contrast (Showman & Polvani 2011). When the radiative time constant is comparable to the characteristic timescale for air parcels to advect eastward over a planetary radius, then significant eastward hot spot offset from the substellar point should occur. Overall, our observations are consistent with this body of theory and models, and suggest that HD 209458b likely exhibits equatorial superrotation at photospheric levels.

Figure 6 compares our best-fit model of HD 209458b's observed full-orbit flux to theoretical phase curves generated from the GCMs presented in Showman et al. (2009). The Showman et al. (2009) 3D dynamical models of HD 209458b are the first with non-gray radiative transfer across the entire wavelength range from the visible through the infrared. In our observations, we find a phase curve maximum that occurs 9.6 ± 1.4 hr before secondary eclipse, corresponding to a hot spot shifted 409 ± 60 eastward of the substellar point. The theoretical phase curves generated from the GCM are shown for no-thermal-inversion models (i.e., models lacking TiO and VO), and models with TiO and VO included at 1× and 3× solar abundance. These molecules cause a thermal inversion due to the extreme opacity of TiO and VO in the visible.¹³ Interestingly, our observations exhibit a peak flux in the light curve which leads secondary eclipse by an amount intermediate to that predicted from the no-inversion and thermal-inversion models. Overall, our observed dayside emission (phase ∼0.25 to 0.75 in Figure 6) matches the predictions of the Showman et al. (2009) models well, especially for models including visible absorbers and a dayside inversion.

Zoom In Zoom Out Reset image size

There exist a variety of factors that can influence the offsets of the hot spot for a given hot Jupiter. Generally speaking, a larger atmospheric opacity would move the photospheres to higher altitude (lower pressure), in which the radiative time constant is shorter, leading to a smaller hot-spot offset; and conversely for smaller atmospheric opacity (Dobbs-Dixon & Lin 2008; Showman et al. 2009; Lewis et al. 2010; Heng et al. 2011a). Such opacity variations could result from the gas metallicity (higher metallicity generally implies greater opacity and vice versa), or potentially from hazes, though there is currently no strong evidence that hazes significantly affect the emission on HD 209458b (Deming et al. 2013). Moreover, for specified opacities, the existence of atmospheric frictional drag could lead to slower wind speeds, decreasing the hot spot offset, potentially even to zero if the drag is sufficiently strong (Rauscher & Menou 2012; Showman et al. 2013; Dobbs-Dixon et al. 2012). Such drag could result from Lorentz forces due to the partial thermal ionization at high temperatures (Perna et al. 2010; Rauscher & Menou 2013; Rogers & Showman 2014). Nevertheless, the fact that our observed offset is significantly nonzero (with 6.9σ confidence) and agrees reasonably well with GCM simulations performed in the absence of strong drag at photospheric levels (Figure 6; Showman et al. 2009) suggests that these magnetohydrodynamic effects do not play a dominant role in controlling the hot spot offset for HD 209458b.

HD 209458b's observed phase curve, which represent the first 4.5 μm measurements of its nightside emission, indicates that HD 209458b has a day-to-night temperature contrast of ΔT_obs = 527 ± 46 K. We find that HD 209458b has a smaller contrast at 4.5 μm (A_obs = 0.352 ± 0.031¹⁴) than other hot Jupiters with higher levels of incident flux (e.g., WASP-12b and WASP-18b; Cowan et al. 2012; Maxted et al. 2013, respectively), consistent with the idea that this temperature contrast is driven by insolation (Perez-Becker & Showman 2013) as the radiative time constant decreases with increasing temperature (Showman & Guillot 2002). However, the only day-night temperature contrast measurements that exist for HD 209458b are at 4.5 μm (presented here) and an upper-limit at 8 μm (Cowan et al. 2007). Thus full-orbit phase curve measurements both at additional wavelengths, such as 3.6 μm, and for other targets are necessary to confirm the hypothesis of Perez-Becker & Showman (2013).

In addition, HD 209458b's phase curve suggests that the nightside is much cooler than predicted by Showman et al. (2009). Yet, while the location of the 1× and 3× solar abundance inversion phase curve minima (Showman et al. 2009) more closely agree with our observations compared to the no inversion model, the GCM overestimates the nightside flux, thereby underestimating the total day-to-night temperature contrast (ΔT_GCM ≈ 300 K versus ΔT_obs = 527 ± 46 K). The Showman et al. (2009) models similarly overpredict HD 189733b's nightside emission compared to Spitzer/IRAC full phase observations by Knutson et al. (2012). This discrepancy is attributed to disequilibrium carbon chemistry not included in the GCM, in particular quenching, which would increase the abundances of CO and CH₄ at higher altitudes so that measurements would be probing a comparably higher optically thick layer with cooler temperatures. In the HD 209458b abundance profiles of Showman et al. (2009) and Moses et al. (2011), the vertical CO abundance profile is driven by chemical equilibrium, resulting in a relatively constant mixing ratio from 10 to 10⁻⁸ bars. Therefore vertical quenching of CO would have a minimal effect on its abundance profile, suggesting that vertical quenching of CO is likely not the culprit for HD 209458b's cooler observed nightside. However, the Moses et al. (2011) abundance profiles of CH₄ suggest that vertical quenching can increase the CH₄ abundance by nearly an order of magnitude from 1 to 10⁻⁵ bars compared to the equilibrium-driven Showman et al. (2009) profiles. This additional CH₄ could potentially help radiate heat and result in an overall cooler nightside than predicted by the GCM. To determine if quenched CH₄ is causing the nightside cooling, we are in the process of measuring HD 209458b's full-orbit phase curve at 3.6 μm, which overlaps the CH₄ ν₃ band. Analyses of these data with future GCM studies, which include both non-equilibrium chemistry and new hot CH₄ line lists (Yurchenko & Tennyson 2014), will better constrain the properties of these two carbon-bearing species and the nightside cooling mechanism.