Thermal Phase Curves of XO-3b: an Eccentric Hot Jupiter at the Deuterium Burning Limit

We report \textit{Spitzer} full-orbit phase observations of the eccentric hot Jupiter XO-3b at 3.6 and 4.5 $\mu$m. Our new eclipse depth measurements of $1770 \pm 180$ ppm at 3.6 $\mu$m and $1610 \pm 70$ ppm at 4.5 $\mu$m show no evidence of the previously reported dayside temperature inversion. We also empirically derive the mass and radius of XO-3b and its host star using Gaia DR3's parallax measurement and find a planetary mass $M_p=11.79 \pm 0.98 ~M_{\rm{Jup}}$ and radius $R_p=1.295 \pm 0.066 ~R_{\rm{Jup}}$. We compare our \textit{Spitzer} observations with multiple atmospheric models to constrain the radiative and advective properties of XO-3b. While the decorrelated 4.5 $\mu$m observations are pristine, the 3.6 $\mu$m phase curve remains polluted with detector systematics due to larger amplitude intrapixel sensitivity variations in this channel. We focus our analysis on the more reliable 4.5 $\mu$m phase curve and fit an energy balance model with solid body rotation to estimate the zonal wind speed and the pressure of the bottom of the mixed layer. Our energy balance model fit suggests an eastward equatorial wind speed of $3.13 ^{+0.26} _{-0.83}$ km/s, an atmospheric mixed layer down to $2.40 ^{+0.92} _{-0.16}$ bar, and Bond albedo of $0.106 ^{+0.008} _{-0.106}$. We assume that the wind speed and mixed layer depth are constant throughout the orbit. We compare our observations with a 1D planet-averaged model predictions at apoapse and periapse and 3D general circulation model (GCM) predictions for XO-3b. We also investigate the inflated radius of XO-3b and find that it would require an unusually large amount of internal heating to explain the observed planetary radius.


INTRODUCTION
As a transiting planet orbits around its host star, the apparent brightness of the planet varies as seen by a distant observer. Infrared phase variations reveal a planet's response to spatial, diurnal and seasonal forcing. Shortperiod planets on circular orbits are subject to strong tidal interaction with their host star and hence are expected to have zero obliquity and to be tidally spun down into synchronous rotation. This means that they don't experience seasons or diurnal forcing, so their atmospheric circulation is driven by the fixed spatial pattern of stellar irradiation and Coriolis forces, resulting in steady-state circulation patterns (e.g. Showman & Guillot 2002). Their thermal phase curves can therefore be translated into a longitudinal thermal map of the planet.
While most hot Jupiters have circular orbits, a few have been found on eccentric orbits. These gas giants are expected to form either in-situ (Bodenheimer et al. 2000) or beyond the snow line in their protoplanetary disk far from their stellar hosts and later migrate inwards via gas disk (Lin et al. 1996), planet-planet scattering (Rasio & Ford 1996), secular interaction (Wu & Lithwick 2011) or Kozai-Lidov migration (Weidenschilling & Marzari 1996;Naoz et al. 2013). In addition to providing insights into gas giants migration mechanisms, hot Jupiters on eccentric orbits are of particular interest for atmospheric studies.
Unlike their circular counterparts, eccentric hot Jupiters experience time-variable heating such that their phase curve reflects a balance between incoming flux, heat transport efficiency (a combination of rotation and winds), and time required to radiate away energy. Therefore, the variable stellar irradiation experienced by eccentric hot Jupiters allows us to break the degeneracy between the heat transport and radiative timescales that limits our studies of typical hot Jupiters on circular orbit (Langton & Laughlin 2008;Iro & Deming 2010;Cowan & Agol 2011;Kataria et al. 2013). In contrast with short-period planets on circular orbits, exoplanets on eccentric orbits present additional challenges when one attempts to retrieve information about their atmosphere from thermal phase observations. In particular, it is difficult to disentangle the flux variation due to the planet's rotation and the change in stellar irradiation. While the dayside of the planet should experience spatial and temporal variability over the course of an orbit, the phase curve should be relatively stable from one or-bit to the next (Showman et al. 2009a;Lewis et al. 2010;Kataria et al. 2013).
In this paper, we present and analyze the full-orbit 3.6 and 4.5 µm phase curves of XO-3b obtained with the Spitzer Space Telescope. In addition, we analyze an unpublished 3.6 µm Spitzer secondary eclipse observation. We use these observations to constrain the radiative and dynamical properties of the atmosphere of XO-3b. The observation and data reduction are presented in Section 2 and the models and methods are described in Sections 3 and 4. Our results are presented in Section 5. Section 6 summarizes the main conclusions from our analysis and presents ideas for future work.

OBSERVATION AND REDUCTION
We observed two continuous full-orbit phase curves of XO-3b (PI H.A. Knutson,PID: 90032): one in each of the 3.6 µm (channel 1) and 4.5 µm (channel 2) bands of the Infrared Array Camera (IRAC; Fazio et al. 2004) on board the Spitzer Space Telescope (Werner et al. 2004). The observations at 3.6 µm and 4.5 µm were acquired on UT 2013 April 12-16 andUT 2013 May 5-8, respectively. Both observations were scheduled to start approximately 5 hours before the start of a secondary eclipse and end approximately 2 hours after the subsequent secondary eclipse. Due to long-term drift of the spacecraft pointing, the telescope was repositioned approximately every 12 hr in order to re-center the target. Hence the observations at each waveband were separated into 9 Astronomical Observation Requests (AORs). We used the sub-array mode with a 2 second frame time (effective exposure time of 1.92 s) for 3.56 days in each waveband which yielded a total of 2400 datacubes in each channel. Every datacube consists of 64 frames of 32×32 pixels (39"×39") resulting in a total of 153,600 images in each waveband.
As explained later, the 3.6 µm phase observations exhibit strong detector systematics during one of the secondary eclipse that biases the eclipse depth measurements. To better constrain the 3.6 µm eclipse depth, we also analyze the eclipse portion of XO-3b's 3.6 µm partial phase curve (PI: P. Machalek, PID 60058). The 3.6 µm phase observations were acquired on UT 2010 March 21-23. Again, the sub-array mode with a 2 second frame time (effective exposure time of 1.92 s) was used and a total 419 datacubes were included in our analysis.

Data Reduction
We use the Spitzer Science Center's basic calibrated data (BCD), which have been dark-subtracted, flat-fielded, linearized and flux calibrated using version S19.2.0 of the IRAC software pipeline. Deming et al. (2011) first noted that the post-cryogenic Spitzer /IRAC data collected in sub-array mode exhibit frame-dependent background flux: they display a settling effect over the 64 frames as well as a sudden increase or decrease in background flux at the 58 th frame. Dang et al. (2018) found that data calibrated using the S19.2.0 pipeline exhibit frame-dependent flux modulation introduced during the sky dark subtraction stage.  The top and bottom panels show the rootmean-squared (RMS) scatter for various photometry schemes performed on 3.6 µm and 4.5 µm data, respectively. The schemes resulting in the smallest RMS scatter are hard-edge circular apertures with a radius of 3.0 pixels for the 3.6 µm observations and a soft-edge aperture of 2.4 pixel for the 4.5 µm observations; these are denoted by a gray circle in each panel.
The 58 th frame error and the flux modulation are both fixed using correction image stacks for each AOR; these were provided by the IRAC team.
We use the Spitzer Phase Curve Analysis Pipeline 1 (SPCA; Dang et al. 2018;Bell et al. 2021), an opensource, modular, and automated pipeline to reduce and analyze our data. After correcting for frame-dependent systematics, we convert the pixel intensity from MJy/str to electron counts and obtain time stamps for each exposure using values from each FITS file. We masked NaN pixels and perform a pixel-by-pixel 4σ sigma clip where the standard deviation, σ, for each pixel is determined along its respective datacube. We choose to mask rather than replace sigma-clipped pixels with average values to minimize the manipulation of the data. We then perform a pixel-level sigma-clipping by comparing each pixel with the pixel located at the same coordinate on all 64 frames of the same datacube and masking all 5σ outliers. Frames containing a sigma-clipped pixel located in a 5×5 pixel box centered on the central pixel of the target are discarded entirely.
We then evaluate the level of background flux for each frame by masking all the pixels within a 7×7 pixel box centered on the target and measuring the median pixel intensity of the remaining unmasked pixels. We then perform background subtraction on each frame. Finally, we estimate the centroid coordinates (x 0 , y 0 ) for each frame using the flux-weighted mean along the x and y axes and measure the Point-Spread-Function (PSF) width (σ x0 , σ y0 ) along the x and y axes.

Aperture Photometry
We perform aperture photometry using soft-edge and hard edge circular apertures as defined in Bell et al. (2021) with radii from 1.25 to 7.25 pixels . We center the aperture on the PSF flux-weighted mean centroid of each frame. To determine the best photometric scheme, we calculate the root-mean-squared (RMS) scatter from a smoothed lightcurve by boxcar averaging with a length of 50. Figure 2 shows the resulting RMS for all our considered aperture choices; we select the scheme with the smallest RMS. We find that the best photometric schemes are hard-edge apertures with a radius of 3.0 and soft-edge apertures with a radius 2.4 pixels for the 3.6 µm and 4.5 µm channels, respectively. The raw photometry and PSF metrics are presented in Figure 3.

Pixel Level Decorrelation Photometry
Pixel Level Decorrelation (PLD) models the systematics as a function of the fractional flux measured by each pixel within a stamp (Deming et al. 2015). SPCA's PLD photometry routine takes a 3 × 3 or 5 × 5 pixel stamp centered on the pixel position (15, 15). The cleaning routine applied to each pixel lightcurve is described in Bell et al. (2021).

PSF Diagnostics
As noted by Lanotte et al. (2014) and Challener et al. (2021), sharp fluctuations in the PSF width can alter the photometry. We search for anomalous PSF behavior by exploring the correlation between the PSF centroids and width shown in Figure 4. Inspecting the PSF centroids for the 3.6 µm data, we find 2 distinct centroid clusters: a large cluster containing most of the data and a smaller cluster which corresponds to the centroid of the last AOR. Unfortunately, it is difficult to constrain the instrumental systematics of the smaller cluster due to  Figure 3. Raw photometry at 3.6 µm (left) and 4.5 µm (right). The top panel shows photometry using the preferred extraction scheme. The second and third panels show the x and y centroid coordinates, respectively. The fourth and fifth panels show PSF width along the x and y axes, respectively. The colored data points represent those used in the analysis while the gray points are discarded from the analysis. The vertical dark colored line represents the time of periastron passage. The vertical dashed lines represent the AOR breaks when the pointing of the spacecraft is readjusted. Note that the data were binned by 64-frame data cube for better visualization. . PSF metrics for the 3.6 µm (left) and the 4.5 µm (right) observations. Each point represent the median of a datacube. The grey dots in the 3.6 µm plot represent data that were discarded from our analysis due to deviant PSF centroids or widths. The darker points represent the in-eclipse and in-transit points data. the sparse data covering the area. We therefore elect to discard the last AOR of 3.6 µm observations. As seen in Figure 4, for both channels, the PSF width follows a parabolic function of centroid. However, with closer inspection of the PSF size plotted against the centroid along the x-direction of the 3.6 µm data, there is a significant deviation from the parabolic trend marked in gray. These deviant points corresponds to the greyed out points in the second to last AOR in Figure 3. Note that this spike in PSF size coincides with a v-shape in the photometry and occurs during a secondary eclipse. The flux decrement is still seen using a larger aperture, which rules out the hypothesis that the inflated PSF causes some of the flux to fall outside the aperture. This detector systematic therefore remains unexplained but is presumably electronic in origin. We incorporate the PSF size into our detector sensitivity model as explained in Section 3, but we were still unable to completely model out the instrumental signal. Upon analyzing the 2010 partial phase curve, we find that the spike feature led to an over-estimation of the eclipse depth. For this reason, we opt to discard these deviant points from our analysis.
Additionally, we elect to discard the first AOR from each phase curves containing 12 datacubes since the target was placed on a different pixel than the rest of the dataset for calibration purposes. After data reduction, a total of 2191 datacubes and 2388 datacubes are kept for analysis for the 3.6 µm and 4.5 µm channels, respectively. The products of our photometry extraction are shown in Figure 3 and 4. We then bin each dataset by datacube (64 frames) to reduce the computational cost of fitting the data with our many different decorrelation models.
3. MODEL SPCA models the photometry as the product of the astrophysical model, A(t), and the detector response, D: F total = A(t) × D. Both models are evaluated simultaneously. We experiment with astrophysical models of varying complexity and with different parametric and non-parametric detector response models, as described below. This experiment results in a statistical analysis to determine which combination of astrophysical and detector model is preferred by the data.

Astrophysical Model
The stellar flux is assumed to be constant except during transit. The shape of the transit is modeled using batman  with quadratic limb darkening. We modeled the astrophysical signal A(t) as the sum of the stellar flux, F (t), and the planetary flux, F p (t): The planetary flux, F p = E(t) × Φ(t), is modeled as a sinusoidal phase variation multiplied by the secondary eclipse, E(t), modeled assuming a uniform disk using batman . We did not account for the light travel time as it is only 45.02 seconds at superior conjunction and does not affect our analysis. The phase variation is modeled as a second order sinusoidal function, Φ, and is expressed following (Lewis et al. 2013a): where θ = f + ω + π/2 is the orbital phase measured from mid-eclipse and f and ω are the true anomaly and the argument of periastron, respectively. We also experiment with a first-order sinusoid with C and D set to 0 to determine the degree of complexity statistically preferred. The phase variation function is scaled such that it is 0 during secondary eclipse and E(t) = F p /F * outside of occultation, where F p /F * is the eclipse depth in terms of stellar flux. This parameterization allows for eclipse depth to be an explicit fit parameter. Despite its simple sinusoidal appearance, this parameterization captures the basic behaviour of more sophisticated simulations: rapid changes in flux near periastron when the planet's orbital phase and temperature both vary quickly, and slower flux evolution near apoastron.

Detector Model
The Spitzer /IRAC 3.6 µm and 4.5 µm channels exhibit significant intrapixel sensitivity variations: for a given astrophysical flux, the electron count varies with the location and spread of the target's PSF on the detector (Charbonneau et al. 2005;Lanotte et al. 2014). Over the years, several decorrelation techniques have been developed to achieve an exquisite level of precision (e.g., Ingalls et al. 2016). We tested several of these methods, namely 2D polynomials, BiLinear Interpolated Sub-pixel Sensitivity (BLISS) mapping, and Pixel Level Decorrelation (PLD). We briefly describe these methods below.

2D Polynomial Model
The 2D polynomial model uses the PSF centroids as regressors and polynomial coefficients are fit parameters (Charbonneau et al. 2005;Cowan et al. 2012). We experiment with second order to fifth order polynomials, including all cross terms.

BLISS Mapping
We also experiment with BLISS mapping, first proposed by Stevenson et al. (2012). In summary, BLISS -- Figure 5. The best-fit models to the 3.6 µm and 4.5 µm observations are presented in the left and right panels, respectively. The vertical dashed lines represent the AOR breaks and the dark vertical lines represent the time of periapse passage. The raw photometry is plotted in the top panels in pale yellow and red while the best-fit signal, F total , is shown in dark yellow and red. The corrected photometry is shown in the second panels in pale yellow and red and the best-fit astrophysical models are shown in darker colors. The third panels are a zoomed-in version of the second panels to more clearly show the phase variations. The last panels shows the residuals. Note that the 3.6 µm and 4.5 µm channels data were fitted independently.

RMS
N binned RMS N binned Figure 6. Red-noise test for our best-fit models as a function of bin size for 3.6 µm (left) and 4.5 µm (right) fit. We binned the data into bins of different size, N binned , and computed the binned residuals RMS. The lighter shaded area is the uncertainty that the MC3 package computed (Cubillos et al. 2017). The grey line represents the expected decrease in RMS if the residuals are purely white noise. We find that there is significant red noise in the residuals of the 3.6 µm fit, while the residuals of the 4.5 µm fit are less correlated. The vertical dashed line represents the number of bins contained in the duration of an eclipse depth.
mapping is a data-driven iterative process to interpolate an intra-pixel sensitivity map of the central pixel, which will in turn be used to decorrelate the data. The area over which the PSF centroids are distributed is divided into subpixels, also called "knots", and each datum is associated with a knot. The data-astrophysical residuals are used to estimate the sensitivity of each knot and the sensitivity of each location of the detector is then estimated by bilinearly interpoating between the surrounding knots. We note that BLISS performs best with continuous uninterrupted Spitzer observations. Since our observation scheme included more AOR breaks than most hot Jupiter phase curves, it may not be the bestsuited dataset for this decorrelation method.

Pixel-Level Decorrelation
In contrast with BLISS mapping and polynomials, Pixel Level Decorrelation (Deming et al. 2015;Luger et al. 2018) does not explicitly depend on the PSF centroids. Rather, this method uses the fractional flux measured by each pixel to model the detector systematics. In principle, astrophysical flux variations would change the intensities of all the pixels in the stamp that encompasses the target's PSF; however, they would not change the fractional fluxes of these pixels in the absence of detector systematics (e.g., spacecraft drift, thermal fluctuation). We experiment with 3 × 3 and 5 × 5 pixel stamps to explore the trade-off between capturing more stellar flux and more background flux. We test fits using a linear PLD and a second-order PLD that does not include cross-terms as they don't improve the quality of the fit (Zhang et al. 2018).

PSF Width
Previous studies have shown that detector models including a function of the PSF width in x and y dramatically improve the photometric residuals (Knutson et al. 2012;Lewis et al. 2013a). Mansfield et al. (2020) and Challener et al. (2021) experimented with linear, quadratic, and cubic dependencies on the PRF width in both the x and y directions and find the linear model to be preferred. Hence, we include the function of PSF width, D PSFw (σ x0 , σ x0 ), as a multiplicative term to the detector models listed above: where d i are the fit parameters. We indeed find that including a PSF shape dependent model significantly reduces the red noise in the residuals when decorrelating data from both Spitzer /IRAC channels.

Step Function
For the 4.5 µm data, we could not get a good fit using any combination of the above detector models: the residuals showed a discontinuity between the first 2 AORs and the following 5 AORs. Often AOR discontinuities can be addressed by also fitting for a linear trend (e.g. Bell et al. 2019), however, an unmodeled linear trend would exhibit a discontinuity at all the AOR breaks, which isn't the case here. To mitigate the problem, we added an additional detector sensitivity model as a multiplicative term: where H(t) is a Heaviside step function, h 1 is the amplitude of the DC offset, which we set as a fit parameter and h 2 is the time of the discontinuity which we fixed to the break between the 2 nd and 3 rd AOR. Our phase curve also includes 2 eclipses observations on either sides of the step to help constrain its magnitude.

PARAMETER ESTIMATION AND MODEL COMPARISON
We use the Affine Invariant Markov Chain Monte Carlo (MCMC) Ensemble Sampler from the emcee package to estimate the parameters and their respective uncertainties . We elect to fix the orbital period, P , the semi-major axis, a, the inclination, i, the eccentricity, e, and the argument of periastron, ω, to the values reported by Wong et al. (2014). We opt to fix the values rather than to impose a Gaussian prior on these parameters to significantly improve the analysis runtime. Although using fixed values instead of distributions can lead to an underestimation of the other model parameter uncertainties, the published uncertainties on the fixed parameters represent less than 0.1% of the reported value and therefore their contribution is negligible.
We initialize the astrophysical parameters to be the values reported by Wong et al. (2014). We begin an initial stage of parameter optimization as described in Bell et al. (2021). We require that transit depths and eclipse depths be between zero and unity. We use the parametrization of Kipping (2013) for the limbdarkening coefficients to ensure that our walkers only explore physically plausible solutions with uniform uninformative sampling.
By default, our pipeline includes a prior rejecting all models with phase variation coefficient that yield negative phase curves (Keating & Cowan 2017). Since an eccentric planet like XO-3b is subject to eccentricity seasons and is expected to have a time-variable atmosphere, we cannot map the planet following Cowan & Agol (2008) and hence cannot apply the more strin- The uncertainties for the fitted parameters to the 3.6 µm data have been inflated by a factor of 1.72 to account for the correlated noise in the residuals (Fig. 6). gent constraint that the implied planetary map is nonnegative (Keating & Cowan 2017;Keating et al. 2019). Nonetheless, our ability to fit the phase variations with an energy balance model in § 5.2 demonstrates that they do not require regions with negative flux. In any case, due to the deep eclipse depths, the fraction of rejected phase curves is less than 0.0001%.
We make the photometric uncertainty, σ F , a fitted parameter and use emcee to estimate the set of parameters that maximizes the log-likelihood: where N dat is the number of data. The badness-of-fit is defined as where F i are brightness measurements and F i,model are the predicted brightness from the astrophysical model described in section 3. We initialize 300 MCMC walkers as a Gaussian ball distributed tightly around our initial guess. We perform an initial burn-in to let the walkers explore a wide region in parameter space during which each walker performs 5000 steps. We then perform a 1000 steps production run while making sure we meet our convergence criteria: 1) the log-likelihood of the best walker does not change over last 1000 steps of the MCMC chain and 2) the distribution of walkers is constant over the last 1000 steps along each parameter. We then obtain a posterior distribution and estimate the 1σ confidence region as the 16 th to 84 th percentile of the posterior distribution of each parameter using all the walkers over the last 1000 MCMC steps.

Model Comparison
As mentioned previously, we experiment with various astrophysical and detector response models. Generally, when the number of fit parameters increase, the fit to the data also improves since the model becomes more flexible. Instead of comparing the badness-of-fit or loglikelihood of the best-fit obtained by each model, we therefore compare the Bayesian Information Criterion fp = 1.77 fp = 1.61 Figure 7. Model comparisons for Spitzer 3.6 µm (top) and 4.5 µm (bottom) phase variations of XO-3b. The best-fit astrophysical model to the Spitzer data obtained using different detector models and phase variation models are shown here. Each column indicates the detector model used. The rows indicate the phase variation model used and whether or not a prior was imposed on the secondary eclipse depth. The ∆BIC from the preferred solution of each fit are indicated in each box and the opacity of the background of each box reflects fit preference (darker is better). The eclipse depth for each fit is denoted by fp, in partsper-mille. Top: The fit using a first-order PLD with a 5 × 5 stamp using a first-order sinusoidal model is preferred. In principle, the shape of the best-fit phase curve is dependent on the models and priors chosen. Fortunately, the shape of the first-order phase variation model is robust against the choice of detector model and prior on the secondary eclipse depth. Bottom: The preferred fit uses a first-order PLD with a 3 × 3 stamp and a second-order sinusoidal model. Note that the shape of the 4.5 µm phase curves is robust against prior and model choices, but the 3.6 µm phase curve shape is not.
(BIC; Schwarz et al. 1978) of the different models: where N par is the number of fit parameters. By definition, a smaller BIC is preferred. A comprehensive model comparison can be found in Figure 7 which shows the shape of the astrophysical phase variation for each model fit and their relative BIC. In principle, the more fit parameters a model has, the more flexible it is and the better goodness-of-fit it achieves. The BIC allows us to justify or rule out having a more complex model with more fit parameters. We note that BLISS does not have any explicit fit parameters, instead we use the number of BLISS knots as number of parameters to estimate its BIC. This could explain why BLISS seems to exhibits systematically worst BIC than the other decorrelation methods (Figure 2), however, we note that the BLISS phase curves also have systematically different amplitudes and shapes than that of the phase curves retrieved using different decorrelation methods.

SPCA Fits
As shown in Figure 7, for our 3.6 µm phase observations, the first-order PLD model using a 5×5 pixels stamp and a first-order sinusoidal phase curve was the preferred model. For our 4.5 µm observations, the preferred solution is a first-order PLD model using a 3×3 pixels stamp and a second-order sinusoidal phase curve. Both preferred fits are shown in detail in Figure 5. We note that the 3.6 µm fit leaves noisier residuals than the 4.5 µm. Figure 6 shows a red-noise test performed on the residuals of the 3.6 µm and 4.5 µm fits. The black line on both plots represent the expected decrease in root-mean-squared (RMS) scatter when uncorrelated data are binned. While the 4.5 µm RMS is in good agreement with this line, the larger than expected RMS of the 3.6 µm channel is indicative of leftover detector or astrophysical variations that our model could not fit.
The best-fit astrophysical parameters for the Spitzer observations are presented in Table 1. We find a ratio of planet to star radius of R p /R * = 0.0866 +0.0014 −0.0012 and 0.0891 ± 0.0008 with the 3.6 µm and 4.5 µm observations, respectively. We measure secondary eclipse depths of 1770 ± 180 ppm and 1610 ± 70 ppm at 3.6 µm and 4.5 µm, respectively. The 4.5 µm eclipse depth is within 1σ of that reported by Wong et al. (2015) and Ingalls et al. (2016). While both phase curves were analyzed independently of each other, their shapes are similar (Figure 8). The 3.6 µm and 4.5 µm phase curves peak 2.11 ± 0.09 and 2.17 ± 0.03 days after transit, respectively. The minimum of the phase curves occur 0.04 ± 0.10 and 0.27 ± 0.04 days before transit. We also note that the 3.6 µm and 4.5 µm best-fit transit time differs by ∆t 0 = 0.0013 ± 0.0005 days. This less than 3 σ discrepancy is about the temporal resolution of our phase curve, hence the difference in the 3.6 and 4.5 µm t 0 isn't meaningful.

Energy Balance Model
To extract radiative and advective properties of the atmosphere of XO-3b, we use the Bell EBM semi-analytical energy balance model (EBM) to fit our detrended phase variations (Bell & Cowan 2018). As the atmospheric temperature does not exceed 2500 K, we did not include the effect of recombination and dissociation of hydrogen. The model treats the advection of atmospheric gas via solid body rotation at angular velocity ω wind , assumes a uniform Bond albedo, A B , and a uniform atmospheric pressure, P 0 , at the bottom of the mixed layer (the part of the atmosphere that responds to diurnal and seasonal forcing). These quantities are treated as fit parameters and remain constant throughout the orbit. Depending on the efficiency of turbulent mixing (Youdin & Mitchell 2010;Bordwell et al. 2018) and varying barotropic largescale flows vertical thickness, P 0 could be deeper in the atmosphere than the pressure at which incoming optical light is absorbed.
Given the significant correlation in the 3.6 µm residuals and that different wavelengths probe different photospheres (Dobbs-Dixon & Cowan 2017), we opt to evaluate each phase curve separately. The energy balance model is fit to the phase curves using the MCMC package emcee (Foreman-Mackey et al. 2013). Again, we decide to fix the orbital period, P , the semi-major axis, a, the inclination, i, the eccentricity, e, and the argument of periastron, ω, with values reported in Table 2. We use the stellar effective temperature from Gaia DR2 (Gaia Collaboration et al. 2018) as well as the updated radius and mass of XO-3 based on Gaia DR2 parallaxes (Stassun et al. 2017, Stassun et al. in prep).
The Bell EBM is agnostic about the underlying rotation of XO-3b, since its atmosphere is not expected to remain stationary with respect to the deeper regions. Nonetheless, in order to convert the inertial-frame atmospheric angular frequency, ω atm , to a zonal wind velocity, we must assume a rotational frequency for the interior of the planet. The interiors of short-period eccentric planets are expected to be pseudo-synchronously rotating. Roughly speaking, this means that the planet is momentarily tidally locked near periapse passage, when the tidal forces are strongest. We adopt the prescription of Hut (1981) for the pseudo-synchronous rotation frequency, ω ps 0.8ω max , where the maximum orbital ---- In both panels, the grey transparent curve represent the limiting case with a short advective timescale τ adv τ rad P and the pink transparent line represents the limiting case of a short radiative timescale, τ rad = 0. angular velocity at periastron is (Cowan & Agol 2011): The equatorial wind velocity is therefore v wind (ω atm − ω ps )R p .
We transform the Spitzer phase curve into an orbital apparent brightness temperature profile and fit it with our energy balance model. Physically motivated uniform priors are imposed to the fit parameters (see Table 2). Our initial attempts to fit the 3.6 µm Spitzer phase-dependent temperatures with the Bell EBM could not reproduce the high brightness temperatures at all orbital phases. Due to the eccentricity of XO-3b's orbit, the planet is expected to experience tidal heating. We therefore add an internal energy source flux term, E int , to equation (1) of Bell & Cowan (2018) and fit for this extra parameter. We experiment with and without this term and find that models allowing for an internal energy source are preferred. The best-fit EBM models to each channel are presented in Table 2 and Figure 8. For comparison, we also show the temperature curve of a limiting case with a very short advective timescale τ adv τ rad P , such that the incident energy is uniformly redistributed instantaneously (the advectiondominated phase curve). We also show the special case where τ rad = 0, i.e., the radiation-dominated phase curve. Neither of these limiting cases account for the presence of an internal energy source.

4.5 µm EBM Fit
The energy balance model is able to reproduce the timing and amplitude of the 4.5 µm phase curve peak and the minimum is only 1σ from the Spitzer observations. The best-fit model has v wind of 3.13 +0.26 −0.83 km/s: approximately the speed of sound. The model suggests that the mixed layer extends down to P 0 = 2.40 +0.09 −0.16 bar and a Bond albedo A b = 0.106 +0.008 −0.106 . From these, we estimate an advective timescale of ∼ 1 hour and an average radiative timescale of ∼ 30 hours. Unlike the 3.6 µm phase curve, we do not need internal heating to explain the 4.5 µm phase curve.

Strong Detector Systematics
Spitzer /IRAC's 3.6 µm channel is known to be less stable than the 4.5 µm channel: stronger detector systematics generally plague the 3.6 µm observations (e.g. Zhang et al. 2018). As discussed in Section 2.3, the 3.6 µm observations exhibit a sharp PSF fluctuation coinciding with one of the secondary eclipses. We experiment with and without excluding the anomalous observations and find an eclipse depth 2130 ± 110 ppm when the aberrant observations are included and 1770 ± 100 ppm when discarded. The deeper eclipse depth estimate is likely a result of the large decrease in flux caused by the sharp PSF width fluctuation, hence we elect to omit the anomalous portion of the observations for the analysis. Consequently, without a reliable second eclipse, it is significantly more difficult to distinguish astrophysical trends from detector systematics. Furthermore, Figure  6 indicates that the residuals from our best-fit model to XO-3b's 3.6 µm phase curve are significantly correlated; on the eclipse duration timescale the 3.6 µm residual RMS is 1.72 times larger than expected if the residuals were uncorrelated. Hence, we elect to inflate our SPCA fit uncertainties by 1.72 for the 3.6 µm fit and the 3.6 µm observations should be interpreted with caution.
- Figure 9. SPCA analysis of the secondary eclipse obtained as part of the 2010 partial phase curve of XO-3b at 3.6 µm with Spitzer. The detrended data are connected with a pale orange line. The dark line represents the best-fit astrophysical model and the orange circles represent the binned calibrated photometry in 20 bins. The bottom panel shows the residuals of the SPCA fit.

3.6 µm Secondary Eclipse Inconsistencies
While the many repeat observations of the 4.5 µm eclipse of XO-3b enable a precise and robust measurements (Wong et al. 2014), the planet has only been observed in 3.6 µm with Spitzer two other times: one secondary eclipse in 2009 (Machalek et al. 2010, PID 525) and an unpublished partial phase curve obtained in 2010 (PI: P. Machalek, PID 60058). Our 1770±180 ppm eclipse depth is 5σ discrepant with the 3.6 µm eclipse depth of 1010 ± 40 ppm reported by Machalek et al. (2010) taken during the cryogenic era Spitzer data. Such significant discrepancies between cryogenic and warm Spitzer eclipse depths are not unheard of (Hansen et al. 2014), and we note that there is visible correlated noise in Machalek et al. (2010)'s 3.6 µm residuals, hence their eclipse uncertainty is likely underestimated.
We investigate the 3.6 µm eclipse depths inconsistencies by analyzing the partial phase curve. Unfortunately, the 2010 observations are difficult to detrend because they only cover half an orbit and have large AOR breaks. Nonetheless, we were able to fit the secondary eclipse portion of the 2010 time series, shown in Figure 9. We find an eclipse depth of 1520 ± 130 ppm, within 2σ our fit to the 2013 full-orbit phase curve. Machalek et al. (2010)'s data were taken with a twochannel mode, i.e., two 2 s exposures at 3.6 µm for every 12 s exposure at 5.8 µm in order to avoid saturating in the shorter wavelength. As a result, the early eclipse has approximately 30% the efficiency of the continuous observation mode used for the later phase curves. Hence, we elect to discard Machalek et al. (2010)'s and the spoiled second eclipse in our full phase curve. When fitting the 2013 phase curve we experiment with using a Gaussian prior centered on the depth we obtained using the 2010 data. We find that the eclipse depth posteriors are consistent with or without the prior.

Tentative 3.6 µm EBM Fit
While the average 4.5 µm Spitzer temperature curve is consistent with the expected T eq of XO-3b, the 3.6 µm temperature curve is higher and does not intersect with the 4.5 µm temperature curve at any point in the orbit. In fact, the 3.6 µm brightness temperature curve is comparable to the planet's irradiation temperature, which means one of two scenarios: 1) a large excess flux or 2) leftover detector systematics. We favour the second scenario-detector systematics-but explore the implications of the model fit below for completeness.
The best-fit model gives an unexpectedly high v wind of 12.1 +3.0 −4.4 km/s, an order of magnitude faster than the equatorial wind speed obtained for the 4.5 µm fit. Such a large equatorial wind speeds are in general not physically possible. As highlighted in Koll & Komacek (2018) although wind speeds approaching or even slightly exceeding the speed of sound (∼2 km/s) are possible in hot Jupiters, the development of shocks and shear instabilities in the atmosphere will naturally limit the maximal wind speeds at the atmospheric pres-sures being probed by our observations of XO-3b. The fit suggests a mixed layer down to P 0 = 0.15 +0.21 −0.11 bar, a short advective timescale of 0.2 hours and an average radiative timescale of ∼ 6 hours. An internal energy source flux of 7.21 +0.11 −0.17 × 10 5 W/m 2 is required to fit the 3.6 µm brightness temperature. This internal flux is 1.3 times the average incident stellar flux. If this is interpreted as tidal heating, it would require a tidal quality factor of Q ∼ 6 × 10 3 . Given the above issue with the 3.6 µm data, these results should be taken lightly. The model predicts that incoming shortwave radiation from the star is deposited entirely at pressures P < 10 bars above the isothermal region. An approximate infrared photosphere is indicated by the red swath from 0.01 to 1 bar.

Vertical Thermal Structure
We present in Figure 10 a planet-averaged 1D model at apoapse and periapse to look at the deposition of stellar energy as a function of pressure (Marley et al. 2002;Fortney et al. 2008). The model uses the code, physics, and chemistry described in Fortney et al. (2008). The specific entropy of the deep adiabat is relatively uncertain, here a value of the intrinsic flux (parameterized as T int ) was set to 300 K a periapse, which sets the adiabat for the self-consistent model in radiativeconvective equilibrium. The temperatures in the deep part of the atmosphere are expected to be horizontally uniform. Therefore, the T int value for apoapse was iterated until the converged apoapse model fell upon the same deep adiabat, yielding T int =520 K. The planet's temperature should be homogenized at depth, but, T int can't be thought of as constant in a simplified 1D modeling framework because the T-P profiles won't lie on the µm eclipse depth rules out the thermal inversion reported by Machalek et al. (2010): water vapour opacity dictates that the 3.6 µm photons originate from deeper in the atmosphere than those at 4.5 µm.
same adiabat due to limitations of a 1D model. We note that a more sophisticated solution with time-stepping atmospheric structure code has been developed to investigate the continuous atmospheric response to the variable incident flux that eccentric planets experience (Mayorga et al. 2021) that uses a different approach by fixing T int . The 1D atmospheric model shown in Figure 10 predicts that incoming shortwave light from the star has all been absorbed by P < 10 bars. This is slightly deeper than the bottom of the mixed layer of P 0 = 2.40 +0.92 −0.16 bar inferred from our EBM fit to the 4.5 µm data. We note that mixed layer here is a relevant model quantity in terms of explaining heat transport at observable levels, but it might not be accurate to extrapolate its meaning to the full depth of the circulation. In reality, because of the barotropic nature of the flow, winds will be strong throughout the observable atmosphere. A hotter interior would result in an adiabat at lower pressures and the isothermal region would take less room in pressure space. The infrared photosphere, on the other hand, is expected to be at the 0.01-1 bar altitude due to the atmosphere's greater IR opacity. In reality, for strong narrow absorption lines, photons could absorbed even down to 10 −4 bars. We only expect winds to become significant at higher altitude than the deposition layer, where horizontal temperature gradients are greater; we therefore expect the mixed layer to lie above the shortwave deposition depth.
By comparing the simulated temperature-pressure profiles with condensation and molecular transition curves in Figure 10, we expect that CO is the main carbon carrier on the planet throughout its orbit. Clouds of TiO/VO, FeSiO 3 and MgSiO 3 might be expected to form at or above the IR photosphere when the planet is near periapse, but would be too deep to affect the emergent spectrum at apoapse. Hence the emergent spectrum of the planet near periapse-possibly including the secondary eclipse-could be affected by the presence of clouds above the notional clear-sky photosphere. There is no evidence of such clouds in the planet's eclipse spectrum. The Spitzer eclipse depth measurements are shown in Figure 11: a deeper 3.6 µm eclipse depth disfavors the dayside thermal inversion reported by Machalek et al. (2010). Unfortunately, our analysis is inconclusive as to temperature inversions for the rest of the orbit due to systematics spoiling our 3.6 µm phase curve. Further investigations with JWST at different orbital phases could provide a better understanding of XO-3b's atmospheric thermal structure and its response to the changing incident stellar flux.

3D General Circulation Model
We present three-dimensional atmospheric circulation models of XO-3b using the SPARC/MITgcm (Showman et al. 2009b). The SPARC/MITgcm couples the MITgcm. a three-dimensional (3D) general circulation model (GCM) (GCM; Adcroft et al. 2004) with a two-stream adaptation of a multi-stream radiative transfer code Marley & McKay (1999). The MITgcm solves the primitive equations using the finite-volume method over a cubed sphere grid. The radiative transfer code solves the two-stream radiative transfer equations, and employs the correlated-k method to solve for upward/downward fluxes and heating/cooling rates through the atmosphere (e.g., Goody et al. 1989;Marley & McKay 1999). The correlated-k method retains most of the accuracy of full line-by-line calculations, while drastically increasing computational efficiency. The SPARC/MITgcm has been applied to a range of exoplanets and brown dwarfs (e.g., Showman et al. 2009b;Kataria et al. 2014;Lewis et al. 2014;Parmentier et al. 2018;Steinrueck et al. 2019).
In this model we adopt a horizontal resolution of 32×64 in latitude and longitude, and a vertical resolution of 40 layers evenly spaced in log pressure from 200 bars at the bottom boundary to 200 microbars at the top. Given that XO-3b is on an eccentric orbit, we assume the planet is "pseudo-synchronously" rotat----- Figure 12. General circulation model (GCM) thermal phase curve predictions for the 3.6 µm (left) and 4.5 µm (right) Spitzer band passes. The coloured swath represent the 1σ uncertainties of the best-fit Spitzer phase curves. The grey swath represents the magnitude of the orbit-to-orbit variability seen in the GCM. The large discrepancy between the 3.6 µm Spitzer phase curve and GCM predictions is likely due to issues with the observations. At 4.5 µm, our model correctly predicts the cooling rate of the planet between eclipse and transit, but underestimates its heating rate between transit and eclipse.The dashed black lines represents the EBM fit to each Spitzer phase curve.
ing, i.e., that the planet's tidal interactions with the star force a single side of the planet to approximately face the star every periapse passage. We estimate the planet's rotation rate following the Hut (1981) formulation for binary stars, T rot = 1.852 × 10 5 seconds (or approximately 2 days). We assume a Solar atmospheric composition without TiO/VO (whose opacities are used to produce a thermal inversion). Given the high gravity and eccentricity of the planet, the GCM was run for 63 Earth days. Despite this short run time, this amounts to approximately 21 orbits of XO-3b, sufficient time for the model to have converged. Unlike hot Jupiters on circular orbits, eccentric hot Jupiters allow us to investigate the atmospheric response to time-varying incident flux. We use our 3D simulations to inspect wind patterns and temperature gradients at apoapse and periapse at pressures of 10 and 100 mbar (Figure 13). At periapse, the GCM predicts a large temperature contrast between the dayside and nightside of the planet with large zonal (east-west) and meridional (north-south) flows from the substellar point to the limbs and antistellar point. Conversely, at apoapse, the atmosphere is comparatively quiescent, with low temperature contrasts from dayside to nightside, and unorganized flow. This behavior is in broad agreement with previous GCMs of highly eccentric exoplanets, including HD 80606b (Lewis et al. 2017) and HAT-P-2b .
Comparing GCM predictions to our Spitzer phase curves in Figure 12, we find that the planetary 3.6 µm flux is greatly underestimated throughout the orbit. As noted in section 5.3, this discrepancy is likely due to issues with the 3.6 µm observations. However, it is sur-prising that the 3.6 µm eclipse depth is also underestimated. Assuming the absence of the formation of a strong thermal inversion in the dayside atmosphere of XO-3b, the relative flux from the planet at 3.6 µm vs. 4.5 µm is a strong function of the pressures and hence atmospheric temperatures being probed by each channel. As the GCM assumes instantaneous equilibrium chemistry in the atmosphere it cannot capture possible disequilibrium processes that may affect abundances of key species such as CH 4 and CO that are strong absorbers in the 3.6 and 4.5 µm Spitzer bandpasses respectively. Visscher (2012) highlights that for eccentric hot Jupiters like XO-3b orbit induced thermal quenching can produce a significant reduction in the abundance of CH 4 in the planet's atmosphere throughout its orbit. Such a scenario would naturally allow the 3.6 µm channel to probe deeper into XO-3b's atmosphere resulting in a deeper than expected secondary eclipse depth in that channel.
The numerical model is able to correctly predict the amplitude of the 4.5 µm phase curve, although the peak and trough occurs slightly later than the Spitzer's data. The model seems to adequately predict the cooling timescale of the planet at 4.5 µm, but underestimates the heating timescale where the discrepancy is more apparent after transit. Although predictions from the GCM match well with the flux measured from XO-b from eclipse through periastron and into transit, especially at 4.5 µm, the GCM predicts a significantly shallower increase in the planetary flux between the transit and eclipse events. As highlighted in other studies of eccentric hot Jupiters such as HAT-P-2b (e.g, Lewis et al. 2013bLewis et al. , 2014 and HD80606b (e.g., de Wit et al. 2016b), Figure 13. Temperature (colorscale) and winds (arrows) at the ≈ 10 mbar and ≈ 120 mbar level of our 1× solar model at snapshots corresponding to apoapse (left) and periapse (right). The substellar longitude is indicated by the solid vertical line. At periapse, the model exhibits a large temperature contrast between the dayside and nightside of the planet. At apoapse, the temperature gradient is attenuate and the predominant wind zonal and meridional flows are suppressed.
the assumption of a "pseudo-synchronous" rotation rate for hot Jupiters on eccentric orbits can result in inconsistencies between model predictions and observations. Near periastron passage, the thermal structure of the planet is dominated by the intense transient heating that results in the theoretical flux from the planet to be fairly insensitive to the assumed rotation rate. However away from periaston passage the assumed rotation rate plays a stronger role in shaping the phase dependent flux from the planet (see discussion in Lewis et al. 2014 in the context of HAT-P-2b and de Wit et al. 2016b in the context of HD80606b). A rotation rate that is slower than the assumed pseudo-synchronous rotation rate would result in the cooler hemisphere of XO-3b being projected toward an earth observer for more the time between the transit and secondary eclipse event that would mimic a slower than expected heating rate for the planet. The energy balance model allows P 0 , A b , E int and V wind to be free-parameters, constant across the planet and throughout the orbit. In contrast, these quantities are spatially and temporally variable in the GCM and the local pressure levels of absorption and re-emission, Bond Albedo and winds are computed self-consistently assuming equilibrium chemistry. Additionally, the EBM is compared with each Spitzer phase curve separately while the simulated GCM phase curves are derived from the same simulation. The internal heat, E int , is a free parameter that is explored in the EBM but not in the GCM. Given the atmospheric temperature expected for XO-3b and assumption of chemical equilibrium, the 3.6 micron photosphere will generally be located at deeper pressures ( 400 mbar) compared to the 4.5 µm photosphere ( 100 mbar). Therefore, increasing internal heat in the GCM could serve to increase the temperature at depth and provide a better prediction to the 4.5 µm phase curve. The discrepancy between the 3.6 µm phase curve and the GCM prediction could also be due instrumental issues with the 3.6 µm channel. However, the GCM also under-predicts the robust 3.6 µm eclipse depth and instrumental effects are unlikely to be the cause for this difference. 5.5. Possible Inflated Radius of XO-3b Figure 14. Theoretical planet radius vs age for XO-3b with typical heating (blue) and with additional heating (red) compared to XO-3b's observed radius. The shaded swatchs represent model uncertainties, dominated by the compositional uncertainties.
Given the significant difference between the radius and mass of XO-3 reported by Stassun et al. (2017) and previously reported parameters, we re-determined these parameters by updating the stellar T eff = 6759 ± 79 K from the latest PASTEL spectroscopic catalog (Soubiran et al. 2020) and the parallax to the Gaia EDR3 value (Gaia Collaboration et al. 2021), then applied all of the other empirical parameters and calculations described in Stassun et al. (2017) with a correction described in Stassun & Torres (2018). We find a stellar mass and radius of M * = 1.21 ± 0.15 M and R * = 1.407 ± 0.038R and a planetary mass and radius of M p = 11.79 ± 0.98 M J and R p = 1.295 ± 0.066R J .
To contextualize these new constraints, we compare them with planetary interior structure models based on Thorngren & Fortney (2018). These are 1-D evolution models that solve the equations of hydrostatic equilibrium, mass conservation, and the relevant equations of state. The most important free parameters are the bulk metallicity Z p and the anomalous heating, which we parametrize as a fraction of the incident flux. Using heating values fitted from the observed hot Jupiter population (Thorngren & Fortney 2018) and a distribution of bulk metallicities inferred from the warm Jupiter population (Thorngren et al. 2016), we predict a range of expected radii for the observed mass and flux. Figure 14 shows this radius plotted against age, and suggests that the observed radius of XO-3b is about 2σ larger than expected. Interestingly the expected radius of XO-3b is consistent with the observed radius of the planet if we use a high 20% insolation interior heating. Such a relationship between the internal heating and planetary radius has been previously proposed for KELT-1b, a 27-M J brown dwarf companion (von Essen et al. 2021). Given the observational uncertainties, this means that we are either measuring the radius at 2 sigma above the true value or the planet is inflated beyond the level expected for its time-averaged flux (Thorngren et al. 2019).
If we interpret the radius of the planet as being linked to internal heating, then there are many possible candidate sources of heat. First, the planet is expected to experience tidal heating and it is possible we are catching the planet in a few Myr window during which it is rapidly circularizing (Mardling & Lin 2002;Ibgui & Burrows 2009;Millholland 2019). Since the observed mass of the planet is less than 1σ below the deuterium limit, it is possible that the true mass of the planet is slightly above the deuterium limit and it could take gigayears to finish burning (Spiegel et al. 2011a;Phillips et al. 2020). Even if somewhat below this limit, it is likely that at least some of the deuterium in the planet has or will be burned (Spiegel et al. 2011b;Bodenheimer et al. 2013); the main issue therefore is whether the heating is sufficient to explain the radius at this age, but that is a difficult question outside the scope of this work. Another possibility is that the planet could be experiencing a Cassini State 2 with high obliquity that would increase the tidal dissipation of the planet (Fabrycky et al. 2007;Adams et al. 2019). Thermal tides, advection of potential temperatures, and Ohmic dissipation are proposed mechanisms for the general radius inflation problem that may also play an important role on XO-3b (e.g. Socrates 2013;Tremblin et al. 2017;Thorngren & Fortney 2018). Ultimately, unusually hot interior would likely be explained by a combination of the usual hot Jupiter inflation effect operating in conjunction with one or more of these other, less common heat sources.

SUMMARY AND CONCLUSION
We presented the analysis of new Spitzer /IRAC observations of the curious XO-3 system harbouring a massive M p = 11.79 ± 0.98 M jup inflated hot Jupiter with R p = 1.295 ± 0.066R jup on a 3.2 day orbit with an orbital eccentricity of e = 0.2769. The full-orbit 3.6 µm and 4.5 µm Spitzer /IRAC phase curves of XO-3b yield a secondary eclipse depths of 1770 ± 180 ppm and 1610 ± 70 ppm at 3.6 µm and 4.5 µm, respectively. From the secondary eclipse portion of the 3.6 µm partial phase curve of XO-3b obtained in 2010 (PI: P. Machalek, PID 60058), we retrieve an eclipse depth of 1520 ± 130 ppm which agrees with our more recent phase curve observations, but is 5σ discrepant with Machalek et al. (2010)'s results. Our observations therefore suggest no evidence for the thermal inversion on the dayside at secondary eclipse proposed by Machalek et al. (2010). The discrepancy is likely due to the less efficient observing mode used by Machalek et al. (2010) and resulting systematics.
Unfortunately, detector systematics are difficult to decorrelate from our 3.6 µm phase curve. Nonetheless, we compare our reliable 4.5 µm phase curve observations to multiple atmospheric models to constraint the radiative and advective properties of XO-3b. We use an energy balance model, assuming the Hut (1981) prescription for pseudosynchronous rotation rate, to fit the more reliable 4.5 µm observations and find a Bond albedo of A b = 0.106 +0.008 −0.106 best-fits our data. We also estimates an average equatorial wind speed v wind of 3.13 +0.26 −0.83 km/s, in agreement with the ∼ 2.5 km/s equatorial wind speeds predicted near periastron by a general circulation model. Our energy balance model fit suggest that the mixed layer of the atmosphere on a planet-averaged extends down to P 0 = 2.40 +0.92 −0.16 bars which is consistent with our 1D radiative transfer model that predicts shortwave light absorbed at deeper pressures. We also compare our phase curves with predictions from a GCM and find good agreement at 4.5 µm and large discrepancies at 3.6 µm. While the disagreement could be due to detector systematics spoiling the 3.6 µm phase curve, it's unlikely to be culprit for the difference in Spitzer -measured and GCM-predicted 3.6 µm eclipse depths. Planetary evolution models suggest that XO-3b is unusually large for its mass. Interestingly, additional heating equivalent to 20% insolation could explain its observed radius. If our results are interpreted as internal heat, the cryptic source of heating could be deuterium burning or tidal dissipation due to the orbital eccentricity or the high planetary obliquity.
Better characterization of stellar properties resulting in stringent constraints on the planet's mass would allow us to determine if the radius of XO-3b is really unusual. Further investigations with the James Webb Space Telescope would enable a search for clouds and could better constrain the presence of a temperature inversion at orbital phases other than at superior conjunction. Phase curve observations at other wavelengths can also better constraint the planetary flux and hence cryptic heating. Along with, HD 80606b, a giant planet with an orbital eccentricity of 0.93, gas giants with moderate orbital eccentricity, such as XO-3b and HAT-P-2b, offer a unique opportunity to characterize the gas giants at different stages of planet migration and help constrain planetary evolution theories.