Constraints on the Atmospheric Circulation and Variability of the Eccentric Hot Jupiter XO3b
Abstract
We report secondary eclipse photometry of the hot Jupiter XO3b in the 4.5 m band taken with the Infrared Array Camera (IRAC) on the Spitzer Space Telescope. We measure individual eclipse depths and center of eclipse times for a total of twelve secondary eclipses. We fit these data simultaneously with two transits observed in the same band in order to obtain a global bestfit secondary eclipse depth of and a center of eclipse phase of . We assess the relative magnitude of variations in the dayside brightness of the planet by measuring the size of the residuals during ingress and egress from fitting the combined eclipse light curve with a uniform disk model and place an upper limit of 0.05%. The new secondary eclipse observations extend the total baseline from one and a half years to nearly three years, allowing us to place an upper limit on the periastron precession rate of degrees/day – the tightest constraint to date on the periastron precession rate of a hot Jupiter. We use the new transit observations to calculate improved estimates for the system properties, including an updated orbital ephemeris. We also use the large number of secondary eclipses to obtain the most stringent limits to date on the orbittoorbit variability of an eccentric hot Jupiter and demonstrate the consistency of multipleepoch Spitzer observations.
1. Introduction
The study of exoplanets has matured greatly in the past decade. The Kepler survey has enriched the body of known exoplanets with thousands of newly discovered objects (Petigura et al. 2013). Meanwhile, the improved capabilities of both space and groundbased telescopes have enabled the collection of highprecision photometry for the brightest of these planets. The majority of the observations carried out to date have focused on the class of gas giant planets known as hot Jupiters, which typically have orbital periods of a few days and atmospheric temperatures ranging between 1000 and 3000 K. These planets are some of the most favorable targets for detailed study, as they have relatively frequent and deep eclipses (e.g., Winn et al. 2008).
By observing the secondary eclipse, which occurs when the planet passes behind its host star, we can deduce the shape of the planet’s infrared emission spectrum and thereby study the properties of its dayside atmosphere (e.g., Charbonneau et al. 2005; Deming et al. 2005). The strong thermal emission of hot Jupiters makes them ideal targets for this type of analysis, and measurements of the secondary eclipse in multiple wavelength bands are commonly used to construct lowresolution emission spectra at infrared wavelengths (e.g., Charbonneau et al. 2008; Knutson et al. 2008). By comparing these observations with theoretical atmospheric models, we can constrain the planet’s atmospheric pressuretemperature profile, chemistry, albedo, and global circulation patterns (see Madhusudhan et al. 2014 for a recent review). Infrared secondary eclipses have been measured for more than fifty exoplanets to date, of which the majority were obtained using the Spitzer Space Telescope (e.g., Knutson et al. 2010; Hansen et al. 2014). Although the telescope’s cryogen was exhausted in 2009, Spitzer’s Infrared Array Camera (IRAC; Fazio et al. 2004) continues to operate in the 3.6 and 4.5 m bands.
Hot Jupiters on eccentric orbits are of particular interest for atmospheric studies. In contrast to hot Jupiters that lie on circular orbits and are therefore likely tidallylocked, these eccentric bodies experience diurnal forcing due to their nonsynchronous rotation as well as significant variations in the incident stellar flux throughout an orbit, leading to timevariable atmospheric forcing (Langton and Laughlin 2008; Iro and Deming 2010; Cowan and Agol 2011; Visscher 2012; Kataria et al. 2013). The response of the planet’s atmosphere reflects a balance between the timevarying stellar irradiation, the gravity wave propagation timescale, and the radiative time scale (e.g., PerezBecker and Showman 2013). These mechanisms can lead to thermal gradients in the planet’s atmosphere, as well as temporal variations in those gradients. Unlike the circular case, phase curve observations of eccentric exoplanets cannot uniquely distinguish between flux variations due to the planet’s rotation and variations due to the changing irradiation experienced by the planet. Observations of the secondary eclipse ingress and egress break this degeneracy by providing a nearinstantaneous picture of the planet’s dayside temperature distribution. Specifically, the apparent offset of the center of eclipse phase due to a zonallyadvected hot spot (Williams et al. 2006; Agol et al. 2010), along with deviations in the eclipse light curve morphology from a uniform disk model during ingress and egress (Rauscher et al. 2007; de Wit et al. 2012; Majeau et al. 2012), can be used to construct a twodimensional map of the planet’s dayside brightness distribution and constrain important properties of the atmosphere. The only successful secondary eclipse mapping observation to date is of the circular hot Jupiter HD 189733b and was carried out independently by de Wit et al. (2012) and Majeau et al. (2012) using the 8 m band on Spitzer.
Recent general circulation models for eccentric hot Jupiters have shown that while the dayside atmospheric brightness of these planets may exhibit some spatial and temporal variability(Langton and Laughlin 2008), the phase curves of these planets should be relatively constant from one orbit to the next (Showman et al. 2009; Lewis et al. 2010; Cowan and Agol 2011; Kataria et al. 2013; Heng and Showman 2014). The presence or absence of longterm variability has important implications for atmospheric studies of these planets. It is common practice to combine data at different wavelengths from separate epochs years apart and assume that these data offer a consistent picture of the dayside emission spectrum. By comparing many individual secondary eclipse observations over a long time baseline, we can directly constrain the magnitude of orbittoorbit variations in the planet’s atmospheric circulation.
In this paper, we analyze twelve secondary eclipse and two transit observations of the eccentric hot Jupiter XO3b obtained in Spitzer’s 4.5 m micron band. This planet has a mass of (JohnsKrull et al. 2008; Hirano et al. 2011) and orbits its host star (spectral type F5V, , K, and ; Winn et al. 2008; Torres et al. 2012) with a period of 3.19 days and an orbital semimajor axis of AU (Winn et al. 2008). With an orbital eccentricity of (Knutson et al. 2014), the stellar flux received by this planet varies by a factor of 3.2 between periapse and apoapse. Machalek et al. (2010) reported secondary eclipse depths for this planet in the 3.6, 4.5, 5.8, and 8.0 micron Spitzer bands, but these data were taken in stellar mode and as a result have a lower cadence and precision than the nowstandard staring mode observations. Here, we present twelve additional secondary eclipse and two transit observations with a total baseline of almost three years.
The paper is organized as follows: The observations and methodology for extracting photometry are described in Section 2. In Section 3, we present the bestfit eclipse parameters for individual secondary eclipse and transit observations as well as the combined eclipse light curve. We then use our results to obtain an updated orbital ephemeris and discuss the implications of our results for the planet’s atmospheric dynamics in Section 4.
2. Observations and photometry
A total of thirteen secondary eclipse and two transit observations of XO3b were carried out using the Infrared Array Camera (IRAC) on the Spitzer Space Telescope in the 4.5 m channel. For all observations, we utilized the IRAC subarray mode, which produced 32 32 pixel (3939) images with 2.0 s integration times. The observations include ten individual secondary eclipse observations, each of which contains 14912 images over 8.4 hr; additionally, two secondary eclipses and one transit are contained within a fullorbit observation, which was obtained on UT 2013 May 58 and has a total duration of 86.7 hr, corresponding to 153600 images (PID: 90032). The long duration of observation and limited onboard memory necessitated multiple breaks for downlinks. One of these downlinks occurred during the second of the two secondary eclipses in the fullorbit observation, which led to difficulties during our analysis. Therefore, this eclipse observation is omitted from the analysis presented here. For the other secondary eclipse and transit contained in this fullorbit observation, we extracted 14912 images from the fullorbit observation spanning the duration of each event to use in our analysis. This number was chosen to match the lengths of the other ten individual eclipse observations; we obtain consistent results when including more or fewer images. The remaining secondary eclipse and transit are contained within a 66 hr observation obtained on UT 2012 April 58 (PID: 60058); as with the fullorbit observation, we extracted 14912 images spanning the secondary eclipse and transit for use in our analysis. In Table 1, we list the start times of the twelve secondary eclipse and two transit observations that are included in our analysis.
Observation Start Time  Median Photometric Aperture  Type  
(BJD2455000)  (pixels)  
Eclipses 
1  294.36810  2.40  fixed 
2  1242.25466  3.37  scale  1.75  
3  1248.66323  2.78  shift  1.2  
4  1251.83395  3.45  scale  1.75  
5  1255.03213  2.98  shift  1.1  
6  1264.60591  2.71  scale  1.70  
7  1270.99427  3.53  scale  1.65  
8  1405.03292  2.41  scale  1.55  
9  1417.74069  3.28  scale  1.80  
10  1430.56660  2.83  scale  1.55  
11  1433.75742  2.82  shift  0.8  
12  1436.94105  3.26  scale  1.95  
Transits 

1  292.22203  2.20  fixed  
2  1418.83774  3.52  scale  1.90 
Data from previous observations using the IRAC instrument displayed a shortduration, ramplike behavior at the start of each observation likely due to the settling of the telescope at a new pointing position (e.g., Knutson et al. 2012). To address this, each individual eclipse observation was preceded by a 30 minute “peak up” observation, which allows the telescope pointing to stabilize before making a final position adjustment to place the star in the middle of the central pixel in the subarray. We exclude these initial peakup observations from our subsequent analysis. We find that our postpeakup light curves do not display any obvious ramplike behavior. Nevertheless, we experimented with trimming the first 10 or 30 minutes of data before fitting to our model eclipse light curve and obtained consistent results with no significant improvement in the residual scatter. In the analysis presented here, we utilize the full eclipse observations.
We extract photometry using methods similar to those used in previous secondary eclipse analyses (e.g., Todorov et al. 2013; O’Rourke et al. 2014). The basic calibrated data (BCD) files are darksubtracted, flatfielded, linearized, and fluxcalibrated using version S19.1.0 of the IRAC pipeline. The exported data from a secondary eclipse or transit observation comprise a set of 233 FITS files, each with 64 images and a UTCbased Barycentric Julian Date (BJD) time stamp designating the start of the first image. The midexposure time stamp for each individual image in a FITS file is calculated assuming uniform spacing and using the difference between the AINTBEG and ATIMEEND headers, which indicate the start and end of each 64image series. We then transform each BJD time stamp into the Barycentric Julian Date based on the Terrestrial Time standard (BJD), using the conversion at the time of our observations (Eastman et al. 2010). The continuous BJD standard is preferred because leap seconds are occasionally added to the BJD standard.
To estimate the sky background, we create a histogram of all pixel values in each image and fit a Gaussian function. We avoid contamination from the wings of the star’s pointspread function (PSF) by excluding pixels within a radius of 15 pixels from the center of the image, as well as the 13th16th rows and the 14th and 15th columns, where the stellar PSF extends close to the edge of the array. In addition, we exclude the top (32nd) row of pixels, since they have values that are consistently lower than those from the rest of the array. Before binning the remaining pixel values and fitting a Gaussian, we iteratively trim outlier values that are more than from the median value. After subtracting the bestfit sky background from the subarray images, we correct for transient “hot pixels” in each set of 64 images by comparing the intensity of each pixel to its median value. Pixel intensities varying by more than from the median value are replaced by the median value. The average percentage of pixels replaced across all pixels contained within an eclipse/transit observation is less than .
The position of the star in each image is determined using fluxweighted centroiding (see, for example, Knutson et al. 2008; Charbonneau et al. 2008): Defining a parameter , we iteratively calculate the fluxweighted centroid within a circular region of radius pixels centered on the estimated position of the star. We optimize separately for each data set to minimize the scatter in the residuals from our eventual bestfit light curve (Section 3). For each image, we also estimate the width of the star’s PSF by computing the noise pixel parameter (Mighell 2005), which is defined in Section 2.2.2 of the Spitzer/IRAC instrument handbook as
(1) 
where is the intensity detected in the th pixel. We define the parameter as the radius of the circular aperture used to calculate and optimize the value for each individual eclipse/transit to minimize the scatter in the residual of the bestfit light curve from each individual fit to an eclipse or transit (Section 3.4).
We perform aperture photometry to calculate the flux in each image, using either a fixed radius ranging from 2.0 to 5.0 in pixels, or a timevarying radius equal to the squareroot of the noise pixel parameter , with either a constant scaling factor or a constant shift (see Lewis et al. 2013 for a full discussion of the noisepixelbased aperture). To determine the optimal choice of aperture photometry for each eclipse or transit, we subtract our bestfit model light curve from each photometric time series and compute the RMS scatter in the resultant residuals, binned in fiveminute intervals, and pick the version of the photometry that gives the minimum scatter. The binning is chosen so as to assess scatter in the data over time scales comparable to the duration of ingress and egress (21.8 minutes), which is more likely to affect the results of the light curve fits. We found that for 11 of the 12 secondary eclipses and one of the two transits, choosing the timevarying aperture photometry yielded the smallest residual scatter. The medians of the timevarying radius for the photometric aperture and types of aperture (fixed, scale, or shift) used for all data sets are listed in Table 1. For the case of timevarying aperture, the scaling factor or number of pixels shifted is also listed.
It is important to note that by choosing the optimal aperture through minimizing the resulting RMS scatter, we may be partially fitting away any deviation during ingress and egress from our eclipse model due to nonuniform dayside brightness. In particular, the pixelmapping method we use to estimate and remove the instrumental noise (Section 3.3) eliminates as much structure as possible from the raw light curve, with no distinction made between astrophysical signals and noise due to intrapixel sensitivity variations. However, there are many flux measurements at a given pixel location before and after ingress and egress, which statistically dominate the pixel map calculation. Moreover, we obtain consistent results using a polynomial noise decorrelation function in the and position of the star centroid, a method that is not capable of removing such smallscale structure in the light curve during ingress and egress. Ultimately, our lack of a priori knowledge of the precise center of eclipse time is much more important in potentially fitting away any astrophysical signal during ingress and egress (see Section 4.2).
Before fitting to the model, we iteratively filter out points in our light curves with uncorrected measured fluxes that vary by more than 3 from the median values in the adjacent 64 frames (i.e., the length of one FITS file) in the time series. Choosing a larger or smaller interval for computing the median values does not significantly affect the number of excised points. For each of the twelve eclipses, the percentage of data points removed is less than 0.5.
3. Data analysis
3.1. Orbital parameters
The routines for fitting the secondary eclipse and transit light curves require precise values for the planet’s orbital parameters. In particular, we need the planet’s orbital period , inclination , eccentricity , scaled semimajor axis , and pericenter longitude , as well as the ratio between the planetary and stellar radii . When optimizing our choice of photometric aperture, we take values from Winn et al. (2008) for all orbital parameters except eccentricity and pericenter longitude, for which we take values from Knutson et al. (2014) ( and ). In the final version of our fits of individual eclipses and transits, we use the updated values for , , and in Table 3 computed from our global fit (see Section 3.4).
3.2. Secondary eclipse/transit model
We calculate the transit and eclipse light curves using the formalism of Mandel and Agol (2002) and include free parameters for the eclipse/transit depth and the center of eclipse/transit time . In addition, we approximate the planet’s phase curve in the region of the secondary eclipse or transit as a linear function of time where the slope (“phase constant”) is a free parameter in our fits. Because this phase curve slope is absent during the secondary eclipse, we keep the observed flux between second and third contacts constant at and scale the amplitudes of ingress and egress appropriately to match them to the outofeclipse phase curve. For fits to individual secondary eclipses we allow the depths and times to vary as free parameters along with the phase constant, while fixing , , and . In the case of transits, the depth is equal to the square of the radius of the planet relative to that of the host star, , which we allow to vary as a free parameter while keeping and fixed. In our global fits, which include both the secondary eclipse and the transit light curves, we allow all of these parameters to vary freely in the fit. The host star XO3 has an effective temperature K and a specific gravity (Torres et al. 2012). When fitting to the transit light curves, we use a fourparameter nonlinear limbdarkening law (Sing 2010) with parameter values calculated for a star with K and (, , , and ).
3.3. Correction for intrapixel sensitivity variations
The fluxes measured within our photometric aperture display a strong correlation with variations in the position of the star on the array. This effect is due to the nonuniform sensitivity of an individual pixel across its area (Charbonneau et al. 2005) and must be accounted for when fitting the fluxes with the eclipse/transit model. In our analysis, we correct for this systematic effect in two different ways.
First, we adopt an approach similar to that used in previous studies (e.g., Todorov et al. 2013) and include a second order polynomial function of the and positions of the star as part of the eclipse fit:
(2) 
We evaluate the utility of individual terms in this expression using the Bayesian information criterion (BIC):
(3) 
where is the number of free parameters in the model (including the parameters in the eclipse model), is the number of data points, and is the standard metric for goodness of fit. Using this metric, we find that the constant and linear terms in the decorrelation function suffice for eclipses 1, 3, 7, and 10. The BIC is minimized by including the term for eclipses 5 and 12, by including the term for eclipses 2, 5, 6, 8, and 9, and by including all the secondorder terms for eclipse 11. It is important to note that the intrinsic error in the light curve, which formally enters into the calculation of , is unknown, and our values are based solely on the spread in the residuals from the fits. However, the large number of data points available ensures we can obtain a reasonably accurate estimate of the noise in each data set, with a correspondingly small uncertainty in our calculation of the BIC. We obtain consistent results for the eclipse depths and times using any decorrelation polynomial of at least first order.
Eclipse/Transit  Eclipse/Transit  Eclipse/Transit  Phase  
Center Time  Phase  Depth  Constant  
(BJD2455000)  ()  ( d)  
Eclipses 
1  294.5736  0.67067  0.1722  
2  1242.4563  0.66987  0.1732  3.0  
3  0.67022  0.1653  2.5  
4  0.66968  0.1561  17.2  
5  1255.2215  0.66957  0.1467  7.9  
6  0.67042  0.1614  5.6  
7  0.66987  0.1587  23.0  
8  1405.2229  0.66937  0.1533  13.4  
9  0.67024  0.1585  2.9  
10  0.66948  0.1578  15.2  
11  1433.9484  0.66991  0.1565  17.0  
12  0.67040  0.1468  11.0  
Transits 
()  
1  292.4320  2.08  0.7830  5.3  
2  1419.0441  1.10  0.7760  7.8 
The second approach we use to remove the intrapixel sensitivity effect is to create a map of the pixel response. Our method is similar to the one described in Ballard et al. (2010). We approximate the star in an image as a point source at the location on the array given by the measured centroid position and model the sensitivity of that location by comparing other images with measured centroid positions near . The effective pixel sensitivity at a given position is calculated using
(4) 
where is the flux measured in the th image, is the intrinsic flux, and are the measured and positions of the star centroid, and and are the standard deviations of the and over the full range in , reflecting the relative spread in position spanned by the points included in the sum. In essence, we are implementing an adaptive smoothing technique where the spatial resolution is smaller in more denselysampled regions. The index in Eq. (4) sums over the nearest neighbors with distance defined as . We chose this number of neighbors to be large enough to adequately map the pixel response and also low enough to remain relatively computationally efficient (Lewis et al. 2013).
The advantage of using the pixelmapping approach is that it uses the measured fluxes themselves to account for the intrapixel sensitivity effect and does not require additional parameters to be fitted, as is the case with the polynomial approximation. This is especially helpful in our global fit to all twelve secondary eclipses and two transits, for which the use of decorrelation functions incurs a prohibitively large computational overhead. Comparing the fitted eclipse parameters for individual eclipses obtained using both approaches, we see that the values and relative uncertainties are consistent. We therefore use the pixelmapping technique in the final version of our analysis.
3.4. Parameter fits
We fit the data for each of the twelve secondary eclipses and two transits, with intrapixel sensitivity correction calculated via pixelmapping (Section 3.3), to our eclipse/transit model (Section 3.2) using a LevenbergMarquardt leastsquares algorithm. The bestfit eclipse/transit parameters for the individual data sets are listed in Table 2; the center of eclipse/transit phase values were calculated using an updated orbital ephemeris derived from fitting all available transit times (Section 4.1). We adjust the center of eclipse phases to correct for the 42.1 s delay (relative to the midtransit time) due to the light travel time across the system (Loeb 2005). The individual secondary eclipse and transit light curves with intrapixel sensitivity effects removed are shown in Figures 1 and 2, respectively, along with the corresponding residuals from the bestfit light curves. The standard deviations of the bestfit residuals for individual secondary eclipse and transit observations are higher than the predicted photon noise limit by a factor of 1.11 to 1.15.
We use two different methods to estimate the uncertainties in our bestfit parameters. First, we estimate the contribution of timecorrelated noise to the uncertainty with the socalled “prayerbead” (PB) method (Gillon et al. 2009): After extracting the residuals from the bestfit solution, we group them into segments of length 14 and cyclically permute their order segment by segment 1000 times, each time adding the new residual series back to the bestfit solution and recomputing the parameters using the leastsquares algorithm. For each of the three parameters, we create a histogram of all 1000 computed values and obtain the uncertainties from the 68 confidence limits. The choice of segment length here is set to be , where is the number of data points (after outliers have been removed); we obtain consistent uncertainties when using different segment lengths. Second, we use a Markov chain Monte Carlo (MCMC) routine with steps to fit the parameters, setting the initial state to be the bestfit solution from the leastsquares analysis and the uncertainty on individual points to be the standard deviation of the bestfit residuals. We discard an initial burnin on each chain of length equal to 10 of the chain length, which we found ensured the removal of any initial transient behavior in a chain, regardless of the choice of initial state. The distribution of values for all parameters are roughly Gaussian and do not display any correlations with one another. As in the PB method, we set the uncertainties in the eclipse/transit parameters to be the 68 confidence limits. For both uncertainty calculation methods, the median of the resulting parameter value distribution lies within 0.1 of the bestfit value.
For each parameter, we choose the larger of the two errors and report it in Table 2. The MCMC method generally yields larger uncertainties for the parameters of individual eclipses than the PB method; the calculated MCMC errors in phase constant are larger for 9 out of 12 eclipses, while the MCMC errors in center of eclipse time and eclipse depth are larger for 8 out of 12 eclipses. The MCMC method yields larger uncertainties for the transit parameters in all cases except for the phase constant of transit 1. Overall, the size of the PB errors for individual eclipse parameters range between 0.80 and 1.15 times that of the corresponding MCMC errors. In the case where the light curves have a significant component of red (i.e., timecorrelated) noise, we would expect the PB uncertainties to be systematically larger than the corresponding MCMC errors (Carter and Winn 2009). However, the characteristic time scale of variations in the residuals of individual eclipse/transit observations is typically on the order of an hour. As a result, cyclic permutation of residuals from an individual eclipse/transit observation, which has a length of around four hours, may not sufficiently sample the timecorrelated noise signal, leading to random variations in the size of the PB errors of up to 20% relative to the MCMC values.
Parameter^{a} 
Value  68 Confidence Limits 
Eclipse depth, ()  0.1580  +0.0033, 0.0039 
Center of eclipse phase,  0.67004  +0.00015, 0.00010 
Phase constant, ( d)  6.0  +1.3, 1.6 
Planettostar radius ratio,  0.08825  +0.00037, 0.00037 
Inclination, (deg)  84.11  +0.16, 0.16 
Scaled semimajor axis,  7.052  +0.076, 0.097 


^{a}The asterisks indicate that these parameter values are computed from the global fit of all twelve secondary eclipses and two transits.
The global bestfit parameters are determined by fitting all twelve secondary eclipses and two transits simultaneously. For each individual secondary eclipse data set, we convert the time series into a phase series and then define a global center of eclipse phase . In the cumulative fit, the planettostar radius ratio is set as an additional global fit parameter that determines the depth of transit as well as the duration of ingress and egress. We also fit for global values of orbital inclination and scaled orbital semimajor axis . The computed global eclipse parameter values are listed in Table 3. The PB errors are up to 50% larger than the MCMC errors with 200,000 steps for all the parameters except for the planettostar radius ratio and inclination, for which the PB errors are 30% and 5% smaller, respectively. As with the individual bestfit eclipse parameters, the larger of the two errors is reported for each parameter. The combined light curve with intrapixel sensitivity effects removed is shown in Figure 3. We estimate the red noise by calculating the standard deviation of the residuals from the bestfit solution for various bin sizes, which is shown in Figure 4 along with the inverse squareroot dependence of white noise on bin size for comparison. We find that on timescales relevant for the eclipses (e.g., the length of ingress/egress – 21.8 minutes), the red noise contributes a relative scatter of roughly 0.01%.
4. Discussion
Fitting all twelve secondary eclipse and two transit observations simultaneously reduces the relative scatter in the binned residuals of the cumulative light curve (Figure 3) to below the 0.05 level, which enables us to better discern any deviations from a spatiallyuniform dayside brightness distribution at the time of secondary eclipse. At the same time, the large number of visits allows us to assess the orbittoorbit variability of the planet’s atmosphere and to calculate improved estimates of the orbital parameters and ephemeris. We discuss each of these aspects separately below.
4.1. Orbital parameters and ephemeris
In our global fit, we allowed the orbital inclination and scaled semimajor axis to vary as free parameters. These two quantities are important for secondary eclipse mapping because they are the primary determinants of the eclipse/transit duration as well as the length of ingress and egress, and small variations in their values can produce residuals during ingress and egress that mimic the expected signals from a planet with a dayside brightness distribution that is not spatially uniform. Conversely, setting the inclination and scaled semimajor axis as free parameters may cause our optimization routine to partially “fit away” any deviation in the light curve due to a genuine spatiallyuniform dayside brightness.
We examine the effect of allowing inclination and scaled semimajor axis values to vary freely by utilizing a MCMC routine with 200,000 steps to calculate the bestfit values and uncertainties of the global eclipse parameters as well as and , for cases where we fit (a) the twelve secondary eclipses alone and (b) all the secondary eclipses and transits together, as was done in the global fit in Section 3.4. For the first case, we obtain degrees and ; for the second case, we obtain degrees and . We see that these values are consistent at the level; the marginalized posterior probability distribution for inclination and scaled semimajor axis are shown in Figure 5 with the and contours marked.
It is evident that the inclusion of the transits in the cumulative fit places a much stronger constraint on the inclination and scaled semimajor axis, which underlines the advantage of including the two transits in our cumulative fit in Section 3.4. Meanwhile, the consistency between the distributions of inclination and scaled semimajor axis for the cases with and without the transits demonstrates that the twelve eclipses alone do appreciably constrain the value of these two parameters. Indeed, the values of the global eclipse depth and center of eclipse phase computed from the MCMC fits with and without transits included are consistent to well within .
We calculate an updated ephemeris for the XO3 system using the transit times listed in Table 2, combined with all previously published values (JohnsKrull et al. 2008; Hébrard et al. 2008; Winn et al. 2008, 2009; Hirano et al. 2011). We define the zeroth epoch as that of the first Spitzer transit observation (transit 1 in this paper) and carry out a linear fit to all of the measured transit times. Our new observations extend the previous baseline by almost a factor of two, and as a result we derive a new, more precise estimate of the planet’s orbital period and zeroth epoch midtransit time:
(5) 
The observed minus calculated transit times using these updated ephemeris values are plotted in Figure 6
We obtain a second, independent estimate of the orbital period by fitting the secondary eclipse times calculated in Machalek et al. (2010) as well as this paper and arrive at a bestfit value of days. This period differs from the bestfit transit period by 2.1. The observed minus calculated secondary eclipse times are plotted in Figure 7, where the zeroth epoch is defined here as the seventh eclipse analyzed in this work.
We can derive limits on the periastron precession rate from the two different estimates of the orbital period. We use the formalism in Pál et al. (2010) to convert the updated estimates for orbital eccentricity and pericenter longitude into a zeroth epoch center of eclipse time. Next, we introduce a term and calculate the predicted center of eclipse times at each of the other epochs of eclipse observations, which we then fit with a line to obtain a new value of the eclipse period. We limit the value of so that the resulting estimates of the eclipse period do not differ from the transit period calculated above by more than 3 and obtain the following constraint on the periastron precession rate:
(6) 
For comparison, the expected periastron precession rate from general relativity and tides for the XO3 system are and deg./day, respectively (Jordán and Bakos 2008).
The most stringent limit we can place on the periastron precession rate is still roughly a factor of ten larger than the largest expected precession rate from theory. In order to reduce the derived limits to values comparable with the theoretical periastron precession rates, we would need to reduce the uncertainty in the eclipse period by roughly a factor of ten, which can be achieved by obtaining more secondary eclipse measurements of XO3b over a sufficiently long time baseline. To assess this possibility, we envisioned a future campaign of eleven secondary eclipse observations with Spitzer, spaced apart in a manner identical to that of the eleven most recent observations analyzed in this work (PID: 90032). Assuming that the timing of the future secondary eclipse observations match the predictions of our previouslycalculated eclipse ephemeris and have uncertainties comparable with those in our data, we shifted the simulated future observations forward in time and calculated the new predicted uncertainty in orbital period. From this, we showed that a total baseline of about 15 years (versus the current baseline of roughly three years) is needed to detect the presence of periastron precession due to general relativity or tides.
Parameter  Value  Units 

RV Model Parameters  
3.19153247  days  
2456419.04365  
0.2769  
347.2  degrees  
1478  m s  
203  m s  
972  m s  
202  m s  
185  m s  
402  m s  
0.023  m sday  
jitter  45.6  m s 
RV Derived Parameters  
0.27005  
0.0612 
The measured center of eclipse times and updated transit ephemeris can be combined with the radial velocity measurements analyzed in Knutson et al. (2014) to arrive at an updated estimate of the orbital eccentricity and pericenter longitude. Using the methodology of that paper, we obtain and degrees. These values are consistent with the eccentricity and pericenter longitude estimates reported in Knutson et al. (2014) ( and degrees). The updated values reduce the uncertainty in the eccentricity value by about a factor of two. The full results of our radial velocity fits are shown in Table 4 and Figure 8. In addition to a new estimate of the orbital period () and center of transit time (), we obtain values for the semiamplitude of the planet’s radial velocity (), the relative radial velocity zero points for data collected by each of the different spectrographs from which radial velocity measurements of the system have been obtained (), the slope () of the bestfit radial velocity acceleration, as well as the radial velocity jitter. A complete description of the methodology used in our radial velocity fits can be found in Knutson et al. (2014).
4.2. Atmospheric circulation
Examining the combined light curve and associated residual series in Figure 3, we do not discern any anomalous signal in the residuals during ingress and egress above the level of the noise (). Therefore, we place an upper limit of 0.05% on the relative magnitude of any deviation from a spatiallyuniform surface brightness. We generate model secondary eclipse light curves using a onedimensional semianalytic model developed in Cowan and Agol (2011) and compare the observed upper limit on dayside brightness variations with the modelpredicted signals in the residuals during ingress and egress.
The atmospheric model takes as inputs the radiative timescale and the rotational frequency of the planet’s atmosphere in units of the periastron orbital angular frequency. We generate simulated light curves for models with (a) (radiative equilibrium; wind velocities are not important in this case and is set to 0), (b) days and , which entails superrotating winds and a hotspot that is offset from the substellar point, and (c) days and , which entails a zonallyuniform brightness. Figure 9 shows the predicted residuals during ingress and egress from the simulated light curves. In order to determine the expected residuals during ingress and egress when fitting the simulated light curves with our uniform disk brightness eclipse model, we consider the case where we fix the center of eclipse time to the model value as well as the case where we let the center of eclipse time vary as a free parameter. In both cases, we allow the eclipse depth to vary and use a thirdorder polynomial to model the outofeclipse flux. The resulting residuals are shown in Figure 10.
The magnitude of the residual signal during ingress and egress is reduced in both cases. In the case where the eclipse time is fixed, the leastsquares algorithm adjusts the eclipse depth to partially compensate for the relatively large deviations from spatiallyuniform brightness, resulting in a slight reduction in the magnitude of the expected residuals and nonzero residuals during eclipse. In the case where the eclipse time is allowed to vary, there is an additional reduction that most strongly affects the residuals from the model with large antisymmetric residual signals during ingress and egress (i.e., days and ). As discussed in de Wit et al. (2012), antisymmetric residuals during ingress and egress can be “fit away” by adjusting the center of eclipse phase. In our cumulative fit of the twelve secondary eclipses, we have imprecise knowledge of the orbital parameters and must allow the eclipse phase to vary as a free parameter. As a result, the maximum expected relative magnitude of the residual signal during ingress and egress is , much less than the noise level. Therefore, we cannot place meaningful constraints on the radiative timescale and wind speeds of XO3b from these data.
When fitting the secondary eclipse to the model light curve, we included the phase constant as a free parameter to describe variation in the observed outofeclipse flux due to the shape of the planet’s phase curve in the region of the secondary eclipse. The global bestfit value, , is distinct from a flat phase curve slope at the level. It is generally difficult to estimate the slope of the phase curve from individual secondary eclipse observations, as this slope is degenerate with the change in flux due to the longterm position drift of the telescope. This is evident from the large scatter in phase constant values around the global bestfit value (see Figure 11). However, by using all twelve secondary eclipse observations in the combined fit, we are able to break this degeneracy and obtain a unique estimate of the phase curve slope in the region of the secondary eclipse.
The observed increase in the planet’s brightness is likely due to an increase in the atmospheric temperature due to increasing stellar irradiation experienced by the planet during secondary eclipse. Another possible contributing factor may be changes in the apparent atmospheric brightness due to planetary rotation, which may be caused by an offset hotspot situated to the west of the substellar point during secondary eclipse. Such an offset hotspot can be generated by either westward winds or the presence of clouds east of the substellar point. Recent atmospheric circulation models do not readily produce westward winds that would lead to a westward offset of the hotspot at infrared wavelengths (Kataria et al. 2013). While nonuniform clouds have been inferred from visible observations of Kepler7b (Demory et al. 2013), the significantly higher temperature and surface gravity of XO3b place the equilibrium cloud decks at pressure levels greater than those probed by the 4.5micron bandpass. Yet another possibility is nonuniform atmospheric chemical composition across the planet. In particular, variations in the abundance of CO in the atmosphere of hot Jupiters has been shown to induce changes in the 4.5micron flux from one side to the other, though the effect on the variation of emergent flux with phase is largely dampened by longitudinal homogenization of chemical abundances through circulation (Agúndez et al. 2014). Therefore, while all of these mechanisms may contribute to the brightening of XO3b at the time of secondary eclipse, it is expected that the overall warming of the planet due to decreasing planetstar distance is the dominant factor. A fullorbit phase curve analysis of XO3b is contained in a parallel study (Cowan et al. 2014, in prep.).
We can use the detection of brightening during eclipse as an independent constraint on the atmospheric properties. In particular, long radiative timescales tend to keep the planet from warming up as it approaches periastron. The increasing flux at superior conjunction, which occurs shortly before periastron, therefore indicates a relatively short radiative timescale. We generate simulated phase curves using the model in Cowan and Agol (2011) for a range of and values. The positive global phase constant rules out radiative timescales longer than one day. This is consistent with the results of a phase curve analysis of another eccentric hot Jupiter, HATP2b, which constrain the radiative timescale to hours (Lewis et al. 2013).
4.3. Orbittoorbit variability
We next consider whether or not the planet displays any evidence for orbittoorbit variability in its dayside brightness. Recent work by Parmentier et al. (2013) has suggested that measurable variations in the eclipse depth of a planet might result from orbittoorbit variations in the abundance of condensable species such as TiO or silicates. In Figure 11, the bestfit eclipse parameters and uncertainties for the twelve individual eclipses are compared with the global bestfit eclipse parameters and uncertainties. The individual eclipse parameter values are largely consistent with the global values, and there is no apparent timecorrelated variation in their values across the twelve eclipses. We find an average deviation in the individual eclipse depths for XO3b at 4.5 m on the order of 5%, which is roughly consistent with the measurement uncertainty. In Machalek et al. (2010), the eclipse depth at 4.5 m was calculated to be , which is consistent with our global eclipse depth at the 2.1 level. The agreement of individual eclipse depths measured at many different epochs suggests that any variability in the eclipse depth that is due to variability in XO3b’s global circulation patterns is likely small (%). Recent threedimensional (3D) atmospheric circulation models that only consider radiative, dynamical, and equilibrium chemical processes find fairly low levels (%) of orbittoorbit variability in the atmospheres of circular (Showman et al. 2009) and eccentric (Lewis et al. 2010; Kataria et al. 2013) exoplanets alike. There are only two planets with comparable observational limits on their orbittoorbit variability: HD 189733b, a hot Jupiter on a circular orbit, and GJ 436b, which is a warm Neptune with a relatively modest orbital eccentricity (; Knutson et al. 2014). While the effects of condensation and/or turbulent mixing in the atmosphere of XO3b may lead to some orbittoorbit variation in the eclipse depths, the relatively low level of variability evident in our data suggest that these processes do not appear to significantly alter the atmosphere’s thermal structure on orbittoorbit (and longer) timescales at the pressure levels probed by the 4.5 m bandpass.
These same data can also be used to evaluate the reliability of the reported uncertainties in secondary eclipse depths measured by the Spitzer Space Telescope. Hansen et al. (2014) argue that errors in the eclipse parameters calculated from Spitzer data may be significantly underestimated, a statement which has major implications for the ability to use secondary eclipse photometry to deduce characteristics of exoplanet atmospheres. In particular, Hansen et al. (2014) argue that individual secondary eclipse measurements do not fully account for detector systematics, which may lead to poor reproducibility of individual eclipse depth measurements and result in apparently inconsistent values for the eclipse depth from one study to the next. The 2.1 discrepancy between our global bestfit eclipse depth and the depth measured from the single 4.5 m observation in Machalek et al. (2010) is consistent with this purported tendency to underestimate the uncertainties of singleeclipse depths. However, the notable selfconsistency of the individual depth measurements for the twelve secondary eclipses analyzed in our work (Figure 11) indicates that our technique of extracting photometry from the observations does adequately account for instrumental effects. Therefore, we expect that the reported uncertainties in the parameter values from both individual eclipse fits and the combined fit are an accurate representation of the measurement uncertainties. This is the first time that such an extensive data set has been collected in the shortwavelength (3.6 and 4.5 micron) Spitzer channels, which are characterized by larger systematic effects due to the intrapixel sensitivity variations than the longer wavelength 5.8 and 8.0 micron channels. Previous multiplevisit studies of HD 189733b and GJ 436b were both at 8.0 microns and were thus not sufficient to definitively assess the validity of the arguments in Hansen et al. (2014).
5. Conclusion
In this paper, we analyzed twelve secondary eclipse observations of the hot Jupiter XO3b in the 4.5 m Spitzer band. After correcting for the intrapixel sensitivity effect, we fit each photometric time series with a uniform disk model light curve and measured the bestfit eclipse depths and center of eclipse times. We included two transits observed in the same band and fit all of the data simultaneously to obtain a global bestfit secondary eclipse depth of and a center of eclipse phase of , as well as updated values for orbital inclination ( degrees), planettostar radius ratio (), and scaled orbital semimajor axis (). We combined the two transits analyzed in this work with all previouslypublished values and derived a more precise estimate of XO3b’s orbital period ( days). By comparing the transit period with one derived from the secondary eclipse times, we were able to constrain the orbital periastron precession of the planet to between and degrees per day. In addition, we incorporated the measured center of eclipse times and updated transit ephemeris in a radial velocity analysis to arrive at updated orbital eccentricity and pericenter longitude values of and degrees, respectively.
The bestfit eclipse depths and center of eclipse times for individual secondary eclipse observations were found to be consistent with the corresponding global values. The lack of any discernible timecorrelated variation in the eclipse depth values is indicative of a low level () of orbittoorbit variability in the planet’s atmospheric brightness and is consistent with recent threedimensional atmospheric circulation models (Showman et al. 2009; Lewis et al. 2010). Furthermore, the selfconsistency of the individual bestfit eclipse parameter values derived from data collected over a span of more than three years demonstrates the reliability of multipleepoch Spitzer data and our ability to adequately account for instrumental effects.
Our cumulative fit reduced the relative scatter in the binned residuals from the bestfit solution to less than 0.05%, and we did not observe any signals in the residual series during ingress and egress above the noise level. We therefore conclude that any deviations from a spatiallyuniform surface brightness on the planet’s dayside must have a relative magnitude of less than 0.05%. While the maximum expected residual signals during ingress and egress from plausible atmospheric models are smaller than the noise level achieved in our observations, our 4 detection of brightening at the time of eclipse enabled us to constrain the atmospheric radiative time scale of XO3b to day. More stringent constraints on XO3b’s orbital eccentricity from radial velocity fits could eventually yield predictions for the center of eclipse times that are precise enough to circumvent the problem encountered in our analysis, where the magnitude of the expected residual signal during ingress and egress was reduced by about an order of magnitude when we allowed the secondary eclipse time to vary in the fits. Meanwhile, an analysis of the available fullorbit phase curve of XO3b could be combined with the results presented here to construct a more detailed picture of the planet’s atmospheric properties. Lastly, a similarlysized set of 4.5 m Spitzer secondary eclipse observations exists for another eccentric hot Jupiter, HATP2b, which will provide a second, independent look at the issues of atmospheric circulation and variability in the eccentric orbit regime.
Footnotes
 affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA
 affiliation: iwong@caltech.edu
 affiliation: Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA
 affiliation: Center for Interdisciplinary Exploration and Astrophysics (CIERA), Department of Earth & Planetary Sciences, Department of Physics & Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
 affiliation: Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
 affiliation: Sagan Fellow
 affiliation: Department of Astronomy, University of Washington, Seattle, WA 98195, USA
 affiliation: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA
 affiliation: Department of Astronomy, University of Maryland, College Park, MD 20742, USA
 affiliation: Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95604, USA
 affiliation: Institute for Astronomy, University of Hawaii, Honolulu, HI 96822, USA
 affiliation: Department of Physics, Principia College, Elsah, IL 62028, USA
 affiliation: Department of Astronomy and Astrophysics, University of California at Santa Cruz, Santa Cruz, CA 95604, USA
 affiliation: Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA
References
 E. Agol, N. B. Cowan, H. A. Knuston, D. Deming, J. H. Steffen, G. W. Henry, and D. Charbonneau. The climate of HD 189733b from fourteen transits and eclipses measured by Spitzer. ApJ, 721:1861, 2010.
 M. Agúndez, V. Parmentier, O. Venot, F. Hersant, and F. Selsis. Pseudo 2D chemical model of hotJupiter atmospheres: application to HD 209458b and HD 189733b. A&A, 564:A73, 2014.
 S. Ballard, D. Charbonneau, D. Deming, H. A. Knuston, J. L. Christiansen, Holman M. J., D. Fabrycky, S. Seager, and M. F. A’Hearn. A search for a subEarth sized companion to GJ 436 and a novel method to calibrate warm Spitzer IRAC observations. PASP, 122:1341, 2010.
 J. A. Carter and J. N. Winn. Parameter estimation from timeseries data with correlated errors: A waveletbased method and its application to transit light curves. ApJ, 704:51, 2009.
 D. Charbonneau, L. E. Allen, S. T. Megeath, G. Torres, R. Alonso, T. M. Brown, R. L. Gilliland, D. W. Latham, G. Mandushev, F. T. O’Donovan, and A. Sozzetti. Detection of thermal emission from an extrasolar planet. ApJ, 626:523, 2005.
 D. Charbonneau, H. A. Knutson, T. Barman, L. E. Allen, M. Mayor, S. T. Megeath, D. Queloz, and S. Udry. Broadband emission spectrum of HD 189733b. ApJ, 686:1341, 2008.
 N. B. Cowan and E. Agol. A model for thermal phase variations of circular and eccentric exoplanets. ApJ, 726:82, 2011.
 J. de Wit, M. Gillon, B.O. Demory, and S. Seager. Towards consistent mapping of distant worlds: Secondaryeclipse scanning of the exoplanet HD 189733b. A&A, 548:A128, 2012.
 D. Deming, S. Seager, L. J. Richardson, and J. Harrington. Infrared radiation from an extrasolar planet. Nature, 434:740, 2005.
 B.O. Demory, J. de Wit, N. Lewis, J. Fortney, A. Zsom, S. Seager, H. Knutson, K. Heng, N. Madhusudhan, M. Gillon, T. Barclay, J.M. Desert, V. Parmentier, and N. B. Cowan. Inference of inhomogeneous clouds in an exoplanet atmosphere. ApJL, 776:L25, 2013.
 J. Eastman, R. Siverd, and B. S. Gaudi. Achieving better than 1 minute accuracy in the Heliocentric and Barycentric Julian Dates. PASP, 122:935, 2010.
 G. G. Fazio, J. L. Hora, L. E. Allen, M. L. N. Ashby, P. Barmby, L. K. Deutsch, J.S. Huang, S. Kleiner, M. Marengo, S. T. Megeath, G. J. Melnick, M. A. Pahre, B. M. Patten, J. Polizotti, H. A. Smith, R. S. Taylor, Z. Wang, S. P. Willner, W. F. Hoffmann, J. L. Pipher, W. J. Forrest, C. W. McMurty, C. R. McCreight, M. E. McKelvey, R. E. McMurray, D. G. Koch, S. H. Moseley, R. G. Arendt, J. E. Mentzell, C. T. Marx, P. Losch, P. Mayman, W. Eichhorn, D. Krebs, M. Jhabvala, D. Y. Gezari, D. J. Fixsen, J. Flores, K. Shakoorzadeh, R. Jungo, C. Hakun, L. Workman, G. Karpati, R. Kichak, R. Whitley, S. Mann, E. V. Tollestrup, P. Eisenhardt, D. Stern, V. Gorjian, B. Bhattacharya, S. Carey, B. O. Nelson, W. J. Glaccum, M. Lacy, P. J. Lowrance, S. Laine, W. T. Reach, J. A¿ Stauffer, J. A. Surace, G. Wilson, E. L. Wright, A. Hoffman, G. Domingo, and M. Cohen. The Infrared Array Camera (IRAC) for the Spitzer Space Telescope. ApJS, 154:10, 2004.
 M. Gillon, B. Smalley, L. Hebb, D. R. Anderson, A. H. M. J. Triaud, C. Hellier, P. F. L. Maxted, D. Queloz, and D. M. Wilson. Improved parameters for the transiting hot Jupiters WASP4b and WASP5b. A&A, 267:259, 2009.
 C. J. Hansen, J. C. Schwartz, and N. B. Cowan. Broadband eclipse spectra of exoplanets are featureless. MNRAS, accepted.
 G. Hébrard, F. Bouchy, F. Pont, B. Loeillet, M. Rabus, X. Bonfils, C. Moutou, I. Boisse, X. Delfosse, M. Desort, A. Eggenberger, D. Ehrenreich, T. Forveille, A.M. Lagrange, C. Lovis, M. Mayor, F. Pepe, C. Perrier, D. Queloz, N. C. Santos, D. Ségransan, S. Udry, and A. VidalMadjar. Misaligned spinorbit in the XO3 planetary system? A&A, 488:763, 2008.
 K. Heng and A. P. Showman. Atmospheric dynamics of exoplanets. Annu. Rev. Earth Planet. Sci., submitted 2014.
 R. Hirano, N. Narita, B. Sato, J. N. Winn, W. Aoki, M. Tamura, A. Taruya, and Y. Suto. Further observations of the tilted planet XO3: A new determination of spinâorbit misalignment, and limits on differential rotation. PASJ, 63:L57, 2011.
 N. Iro and L. D. Deming. A timedependent radiative model for the atmosphere of the eccentric exoplanets. ApJ, 712:218, 2010.
 C. M. JohnsKrull, P. R. McCullough, C. J. Burke, J. A. Valenti, K. A. Janes, J. N. Heasley, L. Prato, R. Bissinger, M. Fleener, C. N. Foote, E. GarciaMelendo, B. L. Gary, P. J. Howell, F. Mallia, G. Masi, and T. Vanmunster. XO3b: A massive planet in an eccentric orbit transiting an F5 V star. ApJ, 677:657, 2008.
 A. Jordán and G. Á. Bakos. Observability of the general relativistic precession of periastra in exoplanets. ApJ, 685:543, 2008.
 T. Kataria, A. P. Showman, N. K. Lewis, J. J. Fortney, M. S. Marley, and R. S. Freedman. Threedimensional atmospheric circulation of hot Jupiters on highly eccentric orbits. ApJ, 767:76, 2013.
 H. A. Knutson, D. Charbonneau, L. E. Allen, A. Burrows, and S. T. Megeath. The 3.68.0 m broadband emission spectrum of HD 209458b: Evidence for an atmospheric temperature inversion. ApJ, 673:526, 2008.
 H. A. Knutson, A. W. Howard, and H. Isaacson. A correlation between stellar activity and hot Jupiter emission spectra. ApJ, 720:1569, 2010.
 H. A. Knutson, N. K. Lewis, J. J. Fortney, A. Burrows, A. P. Showman, N. B. Cowan, E. Agol, S. Aigrain, D. Charbonneau, D. Deming, J.M. Désert, G. W. Henry, J. Langton, and G. Laughlin. 3.6 and 4.5 m phase curves and evidence for nonequilibrium chemistry in the atmosphere of extrasolar planet HD 189733b. ApJ, 754:22, 2012.
 H. A. Knutson, B. J. Fulton, B. T. Montet, M. Kao, H. Ngo, A. W. Howard, J. R. Crepp, S. Hinkley, G. A. Bakos, K. Batygin, J. A. Johnson, T. D. Morton, and P. S. Muirhead. Friends of Hot Jupiters I: A radial velocity search of massive longperiod companions to closein gas giant planets. ApJ, 785:126, 2014.
 J. Langton and G. Laughlin. Hydrodynamic simulations of unevenly irradiated Jovian planets. ApJ, 674:1106, 2008.
 N. K. Lewis, A. P. Showman, J. J. Fortney, M. S. Marley, R. S. Freedman, and K. Lodders. Atmospheric circulation of eccentric hot Neptune GJ436b. ApJ, 720:344, 2010.
 N. K. Lewis, H. A. Knuston, A. P. Showman, N. B. Cowan, G. Laughlin, A. Burrows, D. Deming, J. R. Crepp, K. J. Mighell, E. Agol, G. A. Bakos, D. Charbonneau, J.M. Désert, D. A. Fischer, J. J. Fortney, J. D. Hartman, S. Hinkley, A. W. Howard, J. A. Johnson, M. Kao, J. Langton, and G. W. Marcy. Orbital phase variations of the eccentric giant planet HATP2b. ApJ, 766:95, 2013.
 A. Loeb. A dynamical method of measuring the masses of stars with transiting planets. ApJ, 623:L45, 2005.
 P. Machalek, T. Greene, P. R. McCullough, A. Burrows, C. J. Burke, J. L. Hora, C. M. JohnsKrull, and D. L. Deming. Thermal emission and tidal heating of the heavy and eccentric planet XO3b. ApJ, 711:111, 2010.
 N. Madhusudhan, H. Knutson, J. J. Fortney, and T. Barman. Exoplanetary atmospheres. In Protostars and Planets VI, 2014.
 C. Majeau, E. Agol, and N. B. Cowan. A twodimensional infrared map of the extrasolar planet HD 189733b. ApJ, 747:L20, 2012.
 K. Mandel and E. Agol. Analytic lightcurves for planetary transit searches. ApJ, 580:L171, 2002.
 K. J. Mighell. Stellar photometry and astrometry with discrete point spread functions. MNRAS, 361:861, 2005.
 J. G. O’Rourke, H. A. Knutson, M. Zhao, J. J. Fortney, A. Burrows, E. Agol, D. Deming, J.M. Désert, A. W. Howard, N. K. Lewis, A. P. Showman, and K. O. Todorov. Warm Spitzer and Palomar nearIR secondary eclipse photometry of two Hot Jupiters: WASP48b and HATP23b. ApJ, 781:109, 2014.
 A. Pál, G. Á. Bakos, G. Torres, R. W. Noyes, D. A. Fischer, J. A. Johnson, G. W. Henry, R. P. Butler, G. W. Marcy, A. W. Howard, B. Sipöcz, D. W. Latham, and G. A. Esquerdo. Refined stellar, orbital and planetary parameters for the eccentric HATP2 planetary system. MNRAS, 401:2665, 2010.
 V. Parmentier, A. P. Showman, and Y. Lian. 3D mixing in hot Jupiter atmospheres. I. Application to the day/night cold trap in HD 209458b. A&A, 558:A91, 2013.
 D. PerezBecker and A. P. Showman. Atmospheric heat redistribution on hot Jupiters. ApJ, 776:134, 2013.
 E. A. Petigura, G. W. Marcy, and A. W. Howard. A plateau in the planet population below twice the size of Earth. ApJ, 770:69, 2013.
 E. Rauscher, K. Menou, J. Y.K. Cho, S. Seager, and B. M. S. Hansen. Hot Jupiter variability in eclipse depth. ApJ, 662:L115, 2007.
 A. P. Showman, J. J. Fortney, Y. Lian, M. S. Marley, R. S. Freedman, H. A. Knuston, and D. Charbonneau. Atmospheric circulation of hot Jupiters: Coupled radiativedynamical general circulation model simulations of HD 189733b and HD 209458b. ApJ, 699:564, 2009.
 D. K. Sing. Stellar limbdarkening coefficients for CoRot and Kepler. A&A, 510:A21, 2010.
 K. O. Todorov, D. Deming, H. A. Knutson, A. Burrows, J. J. Fortney, N. K. Lewis, N. B. Cowan, E. Agol, J.M. Desert, P. B. Sada, D. Charbonneau, G. Laughlin, J. Langton, and A. P. Showman. Warm Spitzer photometry of three Hot Jupiters: HATP3b, HATP4b and HATP12b. ApJ, 770:102, 2013.
 G. Torres, D. A. Fischer, A. Sozzetti, L. A. Buchhave, J. N. Winn, M. J. Holman, and J. A. Carter. Improved spectroscopic parameters for transiting planet hosts. ApJ, 757:161, 2012.
 C. Visscher. Chemical timescales in the atmospheres of highly eccentric exoplanets. ApJ, 757:5, 2012.
 P. K. G. Williams, D. Charbonneau, C. S. Cooper, A. P. Showman, and J. J. Fortney. Resolving the surfaces of extrasolar planets with secondary eclipse light curves. ApJ, 649:1020, 2006.
 J. N. Winn, M. J. Holman, G. Torres, P. McCullough, C. JohnsKrull, D. W. Latham, A. Shporer, T. Mazeh, E. GarciaMelendo, C. Foote, G. Esquerdo, and M. Everett. The transit light curve project. IX. Evidence for a smaller radius of the exoplanet XO3b. ApJ, 683:1076, 2008.
 J. N. Winn, J. A. Johnson, D. Fabrycky, A. W. Howard, G. W. Marcy, N. Narita, I. J. Crossfield, Y. Suto, E. L. Turner, G. Esquerdo, and M. J. Holman. On the spinorbit misalignment of the XO3 exoplanetary system. ApJ, 700:302, 2009.