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

Lisa Dang Department of Physics, McGill University, 3600 University St,
Montréal, QC H3A 2T8, Canada  Université de Montréal, Institut de Recherche
sur les Exoplanètes, 1375 Ave.Thérèse-Lavoie-Roux, Montréal, QC H2V 0B3
Taylor J. Bell Bay Area Environmental Research Institute, NASA Ames Research
Center Moffett Field CA 94035  Nicolas B. Cowan Department of Physics, McGill
University, 3600 University St, Montréal, QC H3A 2T8, Canada  Université de
Montréal, Institut de Recherche sur les Exoplanètes, 1375 Ave.Thérèse-Lavoie-
Roux, Montréal, QC H2V 0B3  Department of Earth and Planetary Sciences, McGill
University, 3450 rue University, Montréal, QC H3A 2A7, Canada  Space Science
Institute, 4765 Walnut St, Suite B, Boulder, CO 80301  Daniel Thorngren
Université de Montréal, Institut de Recherche sur les Exoplanètes, 1375
Ave.Thérèse-Lavoie-Roux, Montréal, QC H2V 0B3  Tiffany Kataria  NASA Jet
Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA  Heather
A. Knutson Division of Geological & Planetary Sciences, California Institute
of Technology, 1200 E California Blvd MC 150-21, Pasadena, CA 91125 USA
Nikole K. Lewis  Department of Astronomy and Carl Sagan Institute, Cornell
University, 122 Sciences Drive, Ithaca, NY 14853, USA  Keivan G. Stassun
Vanderbilt University, Department of Physics & Astronomy, 6301 Stevenson
Center Lane, Nashville, TN, 37235, USA  Jonathan J. Fortney  Department of
Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064,
USA  Eric Agol Astronomy Department, University of Washington, Box 351580,
Seattle, WA 98195 USA  Gregory P. Laughlin  Department of Astronomy, Yale
University, 52 Hillhouse Avenue, New Haven, CT 06511-8101, USA  Adam Burrows
Department of Astrophysical Sciences, Princeton University Princeton, NJ 08544
USA  Karen A. Collins  Vanderbilt University, Department of Physics &
Astronomy, 6301 Stevenson Center Lane, Nashville, TN, 37235, USA  Drake Deming
Department of Astronomy, University of Maryland, College Park, MD 20742-2421
USA  Diana Jovmir  Department of Physics, McGill University, 3600 University
St, Montréal, QC H3A 2T8, Canada  Jonathan Langton  Physics Department,
Principia College, 1 Maybeck Place, Elsah, Illinois 62028, USA  Sara Rastegar
Department of Earth & Planetary Sciences, Northwestern University, 2145
Sheridan Road, Evanston, IL 60208, USA  Adam P. Showman  Lunar and Planetary
Lab, University of Arizona Tucson, AZ 85721-0092 USA

(Received Jul 28, 2021; Revised Oct 25, 2021; Accepted Oct 31, 2021)

###### Abstract

We report 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 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.

extra-solar planets, atmosphere, photometry

††journal: AJ††facilities: Spitzer Space Telescope (IRAC)††software: astropy
(Astropy Collaboration et al., 2013), emcee (Foreman-Mackey et al., 2013),
batman (Kreidberg, 2015), SPCA (Dang et al., 2018; Bell et al., 2021),
Bell_EBM (Bell & Cowan, 2018)

tablenum SIXtablenum

##  1 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. Short-
period 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 orbit to the next (Showman et
al., 2009a; Lewis et al., 2010; Kataria et al., 2013).

Thermal phase variations have now been published for more than a dozen gas
giants on circular orbits with the Spitzer Space Telescope: CoRoT-2b (Dang et
al., 2018, PID 11073); HAT-P-7b (Wong et al., 2016, PID 60021); HD 149026b
(Zhang et al., 2018, PID 60021); HD 189733b (Knutson et al., 2012, PID 60021);
HD 209458b (Zellem et al., 2014, PID 60021); KELT-1b (Beatty et al., 2019, PID
11095); KELT-9b (Mansfield et al., 2020, PID 14059); KELT-1b (Beatty et al.,
2019, PID 11095); KELT-16b (Bell et al., 2021, PID 14059); MASCARA-1b (Bell et
al., 2021, PID 14059); Qatar-1b (Keating et al., 2020, PID 13038); WASP-12b
(Cowan et al. 2012, PID 70060; Bell et al. 2019, PID 90186); WASP-14b (Wong et
al., 2015, PID 80073); WASP-18b (Maxted et al., 2013, PID 60185); WASP-19b
(Wong et al., 2016, PID 80073); WASP-33b (Zhang et al., 2018, PID 80073);
WASP-43b (Stevenson et al., 2017, PID 11001); WASP-76b (May et al., 2021, PID
13038); and WASP-103b (Kreidberg et al., 2018, PID 11099). The many published
Spitzer thermal phase curves have also enabled various comparative studies
(Adams & Laughlin, 2018; Keating et al., 2019; Bell et al., 2021). Thermal
phase curves have also been observed with the Hubble Space Telescope (e.g.
Stevenson et al., 2014; Kreidberg et al., 2018; Arcangeli et al., 2019, 2021)
and more recently with the Transiting Exoplanet Survey Satellite (e.g. Wong et
al., 2020; Daylan et al., 2021; von Essen et al., 2021). In contrast, only 3
exoplanets with an eccentricity greater than 0.15 have published phase curves:
HAT-P-2b (Lewis et al., 2013a, PID 90032), GJ 436b (Lanotte et al., 2014, PID
30129) and HD 80606b (de Wit et al., 2016a, PID 60102).

Figure 1: Top-view of XO-3b’s orbit with parameters from Wong et al. (2014).
The gray dots represent the planet’s position at 1 hour intervals and the
planet revolves counter-clockwise. The shaded area represents our line of
sight (the Earth is off the bottom of the plot).

XO-3b (Johns-Krull et al., 2008) is an eccentric hot Jupiter ($e=0.2769$)
orbiting a F5V star (Figure 1) and is a tantalizing target for follow-up
observations. With a mass of $M_{p}=11.79\pm 0.98M_{\rm{Jup}}$, XO-3b provides
an important link between giant exoplanets and low-mass brown dwarfs. In
addition, its unusually large radius measurements of $R_{p}=1.295\pm
0.066R_{\rm{jup}}$ is difficult to explain with traditional evolution models
for hot Jupiters with an age of $2.82^{+0.58}_{-0.82}$ Gyr (e.g. Liu et al.,
2008; Winn et al., 2008). Finally, the incident stellar flux on XO-3b at
periapse is 3.3 times that at apoapse, which could cause large variations in
atmospheric temperature, wind speed, chemistry, and clouds.

Naturally, the puzzling system of XO-3 has been the subject of various
observational studies. Planets on eccentric orbits are often misaligned with
their stellar spin. Hébrard et al. (2008) measured a sky-projected spin-orbit
misalignment for the star XO-3 of $\lambda=70\pm 15^{\circ}$ using SOPHIE
observations. This quantity has been revised by Winn et al. (2009) to
$\lambda=37.3\pm 3.7^{\circ}$. Turner et al. (2017) obtained ground-based
transit observations and suggest that the anomalously large transit depth they
measure in B-band could be indicative of scattering in its atmosphere.
Machalek et al. (2010) measured eclipses of XO-3b in the four Spitzer/IRAC
wavebands to infer the planet’s vertical temperature profile. They found that
the dayside of the planet exhibited a temperature inversion (temperature
increasing with height over a limited pressure range). Wong et al. (2014) and
Ingalls et al. (2016) analyzed 12 secondary eclipses of XO-3b at 4.5 $\mu$m.
These measurements favor a greater eclipse depth than reported by Machalek et
al. (2010), and hence strengthen the claimed temperature inversion. With the
baseline eclipse observations of XO-3b at $4.5~{}\mu$m extended to 3 years,
Wong et al. (2014) place an upper limit on the periastron precession rate of
$2.9\times 10^{-3}$ deg/day.

In this paper, we present and analyze the full-orbit 3.6 and 4.5 $\mu$m phase
curves of XO-3b obtained with the Spitzer Space Telescope. In addition, we
analyze an unpublished 3.6 $\mu$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.

##  2 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 $\mu$m (channel 1) and 4.5 $\mu$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 $\mu$m
and 4.5 $\mu$m were acquired on UT 2013 April 12-16 and UT 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$\times$32 pixels (39”$\times$39”)
resulting in a total of 153,600 images in each waveband.

As explained later, the 3.6 $\mu$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 $\mu$m eclipse depth, we also
analyze the eclipse portion of XO-3b’s 3.6 $\mu$m partial phase curve (PI: P.
Machalek, PID 60058). The 3.6 $\mu$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.

Figure 2: The top and bottom panels show the root-mean-squared (RMS) scatter
for various photometry schemes performed on 3.6 $\mu$m and 4.5 $\mu$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 $\mu$m observations
and a soft-edge aperture of 2.4 pixel for the 4.5 $\mu$m observations; these
are denoted by a gray circle in each panel.

Figure 3: Raw photometry at 3.6 $\mu$m (left) and 4.5 $\mu$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.

Figure 4: PSF metrics for the 3.6 $\mu$m (left) and the 4.5 $\mu$m (right)
observations. Each point represent the median of a datacube. The grey dots in
the 3.6 $\mu$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.

###  2.1 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 58th 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 58th 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 Pipeline111Details about how to
install and use SPCA can be found at https://spca.readthedocs.io (SPCA; Dang
et al., 2018; Bell et al., 2021), an open-source, 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$ sigma clip where the
standard deviation, $\sigma$, 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$\sigma$ outliers. Frames containing a sigma-clipped pixel
located in a 5$\times$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$\times$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
$(\sigma_{x_{0}},\sigma_{y_{0}})$ along the $x$ and $y$ axes.

###  2.2 Extracting the Photometry

####  2.2.1 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 $\mu$m and 4.5 $\mu$m channels,
respectively. The raw photometry and PSF metrics are presented in Figure 3.

####  2.2.2 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\times 3$ or $5\times 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).

###  2.3 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 $\mu$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 the sparse data covering the area. We therefore
elect to discard the last AOR of 3.6 $\mu$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 $\mu$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
$\mu$m and 4.5 $\mu$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_{\rm total}=A(t)\times 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.

###  3.1 Astrophysical Model

The stellar flux is assumed to be constant except during transit. The shape of
the transit is modeled using batman (Kreidberg, 2015) with quadratic limb
darkening. We modeled the astrophysical signal $A(t)$ as the sum of the
stellar flux, $F_{\star}(t)$, and the planetary flux, $F_{p}(t)$:

$A(t)=F_{\star}(t)+F_{p}(t).$ (1)

The planetary flux, $F_{p}=E(t)\times\Phi(t)$, is modeled as a sinusoidal
phase variation multiplied by the secondary eclipse, $E(t)$, modeled assuming
a uniform disk using batman (Kreidberg, 2015). 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, $\Phi$, and is expressed following (Lewis et al., 2013a):

$\begin{split}\Phi(\theta)=1&+A[\cos(\theta)-1]+B\sin(\theta)\\\
&+C[\cos(2\theta)-1]+D\sin(2\theta),\end{split}$ (2)

where $\theta=f+\omega+\pi/2$ is the orbital phase measured from mid-eclipse
and $f$ and $\omega$ 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.

Figure 5: The best-fit models to the 3.6 $\mu$m and 4.5 $\mu$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_{\rm 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 $\mu$m and 4.5 $\mu$m channels data were fitted independently.

Figure 6: Red-noise test for our best-fit models as a function of bin size for
3.6 $\mu$m (left) and 4.5 $\mu$m (right) fit. We binned the data into bins of
different size, $N_{\rm 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 $\mu$m fit, while the residuals of the 4.5 $\mu$m
fit are less correlated. The vertical dashed line represents the number of
bins contained in the duration of an eclipse depth.

###  3.2 Detector Model

The Spitzer/IRAC 3.6 $\mu$m and 4.5 $\mu$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.

####  3.2.1 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.

####  3.2.2 BLISS Mapping

We also experiment with BLISS mapping, first proposed by Stevenson et al.
(2012). In summary, BLISS 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 best-suited
dataset for this decorrelation method.

####  3.2.3 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\times 3$ and
$5\times 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).

####  3.2.4 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_{\rm PSFw}(\sigma_{x_{0}},\sigma_{x_{0}})$, as a multiplicative term to the
detector models listed above:

$D_{\rm
PSFw}(\sigma_{x_{0}},\sigma_{y_{0}})=d_{0}+d_{1}\sigma_{x_{0}}+d_{2}\sigma_{y_{0}}$
(3)

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.

####  3.2.5 Step Function

For the 4.5 $\mu$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:

$D_{step}=h_{1}[H(t-h_{2})]+1,$ (4)

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 2nd and 3rd AOR. Our
phase curve also includes 2 eclipses observations on either sides of the step
to help constrain its magnitude.

Table 1: Best-fit Model Parameter Estimates 333The uncertainties for the fitted parameters to the 3.6$~{}\mu$m data have been inflated by a factor of 1.72 to account for the correlated noise in the residuals (Fig. 6). Name | Symbol | Prior | Reference | 3.6 $\mu$m | 4.5 $\mu$m
---|---|---|---|---|---
Fitted |  |  |  |  |
Time of transit (BMJD) | $t_{0}$ | – | – | $56418.5434^{+0.0004}_{-0.0005}$ | $56418.5421\pm 0.0002$
Radius of planet | $R_{p}/R_{*}$ | [0,1] | – | $0.0866^{+0.0014}_{-0.0012}$ | $0.0891\pm 0.0008$
Eclipse depth (ppm) | $F_{p}/F_{*}$ | [0,1] | – | $1770\pm 180$ | $1610\pm 70$
Limb Darkening coeff. | $q_{1}$ | [0,1] | – | $0.05^{+0.06}_{-0.05}$ | $0.05^{+0.03}_{-0.02}$
Limb Darkening coeff. | $q_{2}$ | [0,1] | – | $0.2^{+0.6}_{-0.4}$ | $0.8^{+0.3}_{-0.4}$
Phase var. coeff (order 1) | $A$ | – | – | $0.21\pm 0.04$ | $0.34\pm 0.02$
Phase var. coeff (order 1) | $B$ | – | – | $-0.01\pm 0.04$ | $-0.08^{+0.02}_{-0.03}$
Phase var. coeff (order 2) | $C$ | – | – | – | $0.05\pm 0.02$
Phase var. coeff (order 2) | $D$ | – | – | – | $0.07^{+0.01}_{-0.02}$
Fixed |  |  |  |  |
Period (days) | $P$ | 3.19153285 $\pm$ 0.000000058 | Wong et al. (2014) | – | –
Semi-major axis | $a/R*$ | 7.052 ${}^{+0.076}_{-0.097}$ | Wong et al. (2014) | – | –
Inclination | $i$ | 84.11 $\pm$ 0.16 | Wong et al. (2014) | – | –
Eccentricity | $e$ | 0.2769 ${}^{+0.0017}_{-0.0016}$ | Wong et al. (2014) | – | –
Longitude of periapse | $\omega$ | 347.2 ${}^{+1.7}_{-1.6}$ | Wong et al. (2014) | – | –
Derived |  |  |  |  |
Periastron passage after transit (days) | – | 2.5645 | – | – | –
Phase curve peak after transit (days) | – | – | – | 2.11 $\pm$ 0.09 | 2.17 $\pm$ 0.03
Phase curve trough after transit (days) | – | – | – | -0.04 $\pm$ 0.10 | -0.27 $\pm$ 0.04
Phase curve peak after eclipse (hrs) | – | – | – | -0.53 $\pm$ 1.98 | 0.95 $\pm$ 0.82

##  4 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 (Foreman-Mackey et al., 2013). We elect to fix the orbital
period, $P$, the semi-major axis, $a$, the inclination, $i$, the eccentricity,
$e$, and the argument of periastron, $\omega$, 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 limb-darkening 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 stringent
constraint that the implied planetary map is non-negative (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, $\sigma_{F}$, a fitted parameter and use
emcee to estimate the set of parameters that maximizes the log-likelihood:

$\ln(L)=-\frac{1}{2}\chi^{2}-N_{\rm dat}\ln\sigma_{F}-\frac{N_{\rm
dat}}{2}\ln(2\pi),$ (5)

where $N_{\rm dat}$ is the number of data. The badness-of-fit is defined as

$\chi^{2}=\frac{\sum_{i}^{N_{\rm dat}}[F_{i}-F_{i,\rm
model}]^{2}}{\sigma_{F}^{2}},$ (6)

where $F_{i}$ are brightness measurements and $F_{i,\rm 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$\sigma$ confidence region as the
16th to 84th percentile of the posterior distribution of each parameter using
all the walkers over the last 1000 MCMC steps.

###  4.1 Model Comparison

Figure 7: Model comparisons for Spitzer $3.6~{}\mu$m (top) and $4.5~{}\mu 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 $\Delta$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 $f_{p}$, in parts-per-mille. Top: The fit using a first-order
PLD with a $5\times 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\times 3$ stamp and a second-order sinusoidal model. Note that the
shape of the 4.5 $\mu$m phase curves is robust against prior and model
choices, but the 3.6 $\mu$m phase curve shape is not.

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 log-likelihood of the best-fit obtained by
each model, we therefore compare the Bayesian Information Criterion (BIC;
Schwarz et al., 1978) of the different models:

$BIC=N_{\rm par}\ln N_{\rm dat}-2\ln L,$ (7)

where $N_{\rm 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.

##  5 Results

###  5.1 SPCA Fits

As shown in Figure 7, for our 3.6 $\mu$m phase observations, the first-order
PLD model using a 5$\times$5 pixels stamp and a first-order sinusoidal phase
curve was the preferred model. For our 4.5 $\mu$m observations, the preferred
solution is a first-order PLD model using a 3$\times$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 $\mu$m fit leaves noisier residuals than the
4.5 $\mu$m. Figure 6 shows a red-noise test performed on the residuals of the
3.6 $\mu$m and 4.5 $\mu$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 $\mu$m RMS is in good agreement with this line, the
larger than expected RMS of the 3.6 $\mu$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\pm 0.0008$ with the 3.6
$\mu$m and 4.5 $\mu$m observations, respectively. We measure secondary eclipse
depths of $1770\pm 180$ ppm and $1610\pm 70$ ppm at $3.6~{}\mu$m and
$4.5~{}\mu$m, respectively. The 4.5 $\mu$m eclipse depth is within 1$\sigma$
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 $\mu$m and 4.5 $\mu$m phase curves peak $2.11\pm
0.09$ and $2.17\pm 0.03$ days after transit, respectively. The minimum of the
phase curves occur $0.04\pm 0.10$ and $0.27\pm 0.04$ days before transit. We
also note that the 3.6 $\mu$m and 4.5 $\mu$m best-fit transit time differs by
$\Delta t_{0}=0.0013\pm 0.0005$ days. This less than 3 $\sigma$ discrepancy is
about the temporal resolution of our phase curve, hence the difference in the
3.6 and 4.5 $\mu$m $t_{0}$ isn’t meaningful.

###  5.2 Energy Balance Model

Figure 8: The yellow (left) and red (right) line and swath represent the brightness temperature and uncertainties of XO-3b based on the Spitzer 3.6 $\mu$m and 4.5 $\mu$m observations, respectively. The dark yellow and red dashed lines represent the energy balance model fits to the temperature curves. The light grey lines represent the time of secondary eclipse and the dark grey line indicates the transit time. The vertical brown line represents the periastron passage. In both panels, the grey transparent curve represent the limiting case with a short advective timescale $\tau_{\rm adv}\ll\tau_{\rm rad}\ll P$ and the pink transparent line represents the limiting case of a short radiative timescale, $\tau_{\rm rad}=0$. Table 2: Energy Balance Model Fit Parameter | Prior | Reference | 3.6 $\mu$m | 4.5 $\mu$m
---|---|---|---|---
Fitted |  |  |  |
$P_{0}$ (bar) | [$10^{-3}$, 20] | – | $0.15^{+0.21}_{-0.11}$ | $2.40^{+0.92}_{-0.16}$
Equatorial Wind Speed, $V_{\rm{wind}}$, (km/s) | – | – | $12.1^{+3.0}_{-4.4}$ | $3.13^{+0.26}_{-0.83}$
Internal Energy, $E_{\rm{int}}$, ($W/m^{2}$) | [0, $+\infty$] | – | $(7.21^{+0.11}_{-0.17})\times 10^{5}$ | $(1.6^{+11.8}_{-27.0})\times 10^{2}$
Albedo, $A_{B}$ | [0,1] | – | $0.0$ | $0.106^{+0.008}_{-0.106}$
Fixed |  |  |  |
Period, $P$, (days) | 3.19153285 $\pm$ 0.000000058 | Wong et al. (2014) | – | –
Eccentricity, $e$ | 0.2769 ${}^{+0.0017}_{-0.0016}$ | Wong et al. (2014) | – | –
Semi-major axis, $a$, (AU) | 0.04589 ${}^{+0.0063}_{-0.0049}$ | Wong et al. (2014) | – | –
Planet Radius, $R_{p}$, ($R_{\rm{Jup}}$) | 1.295 $\pm$ 0.066 | this work | – | –
Planet Mass, $M_{p}$, ($M_{\rm{Jup}}$) | 11.79 $\pm$ 0.98 | this work | – | –
Inclination, $i$, (deg) | 84.11 $\pm$ 0.16 | Wong et al. (2014) | – | –
Longitude of periapse, $\omega$ | 347.2 ${}^{+1.7}_{-1.6}$ | Wong et al. (2014) | – | –
Pseudo-Synchronous Rotation Period, $P_{\rm{ps}}$ (days) | 2.170 | Hut (1981) prescription | – | –
Stellar Effective Temperature, $T_{\rm{eff}}$, (K) | $6759\pm 79$ | Soubiran et al. 2020 | – | –
Stellar Radius, $R_{*}$, ($R_{\odot}$) | $1.407\pm 0.038$ | this work | – | –
Stellar Mass, $M_{*}$, ($M_{\odot}$) | 1.21 $\pm$ 0.15 | this work | – | –
Derived |  |  |  |
Average Radiative Timescale, $\tau_{\rm{rad,ave}}$, (hrs) | – | – |  $\sim$ 6 |  $\sim$ 30
Advective Timescale, $\tau_{adv}$, (hrs) | – | – |  $\sim$ 0.1 |  $\sim$ 1

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 $\omega_{\rm
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 large-scale 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 $\mu$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, $\omega$, 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, $\omega_{\rm 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,
$\omega_{ps}\simeq 0.8\omega_{\rm max}$, where the maximum orbital angular
velocity at periastron is (Cowan & Agol, 2011):

$\omega_{\rm max}=\frac{2\pi}{P}\frac{(1+e)^{1/2}}{(1-e)^{3/2}}.$ (8)

The equatorial wind velocity is therefore

$v_{\rm wind}\mathbf{\simeq}(\omega_{\rm atm}-\omega_{ps})R_{p}.$ (9)

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 $\mu$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_{\rm{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 $\tau_{\rm adv}\ll\tau_{\rm rad}\ll P$, such that
the incident energy is uniformly redistributed instantaneously (the advection-
dominated phase curve). We also show the special case where $\tau_{\rm
rad}=0$, i.e., the radiation-dominated phase curve. Neither of these limiting
cases account for the presence of an internal energy source.

####  5.2.1 $4.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ EBM Fit

The energy balance model is able to reproduce the timing and amplitude of the
4.5 $\mu$m phase curve peak and the minimum is only 1$\sigma$ from the Spitzer
observations. The best-fit model has $v_{\rm 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 $\sim 1$ hour and an average radiative timescale of $\sim 30$
hours. Unlike the 3.6 $\mu$m phase curve, we do not need internal heating to
explain the 4.5 $\mu$m phase curve.

###  5.3 Spitzer/IRAC $3.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$: A
Cautionary Tale

####  5.3.1 Strong Detector Systematics

Spitzer/IRAC’s 3.6 $\mu$m channel is known to be less stable than the 4.5
$\mu$m channel: stronger detector systematics generally plague the 3.6 $\mu$m
observations (e.g. Zhang et al., 2018). As discussed in Section 2.3, the 3.6
$\mu$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\pm 110$ ppm when the aberrant
observations are included and $1770\pm 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 $\mu$m
phase curve are significantly correlated; on the eclipse duration timescale
the 3.6 $\mu$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 $\mu$m fit and the 3.6 $\mu$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 $\mu 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.

####  5.3.2 $3.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Secondary Eclipse
Inconsistencies

While the many repeat observations of the $4.5~{}\mu$m eclipse of XO-3b enable
a precise and robust measurements (Wong et al., 2014), the planet has only
been observed in 3.6 $\mu$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$\pm$180
ppm eclipse depth is $5\sigma$ discrepant with the 3.6 $\mu m$ eclipse depth
of $1010\pm 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
$\mu$m residuals, hence their eclipse uncertainty is likely underestimated.

We investigate the 3.6 $\mu$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\pm 130$ ppm,
within $2\sigma$ our fit to the 2013 full-orbit phase curve.

Machalek et al. (2010)’s data were taken with a two-channel mode, i.e., two 2
s exposures at 3.6 $\mu$m for every 12 s exposure at 5.8 $\mu$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.

####  5.3.3 Tentative $3.6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ EBM Fit

While the average 4.5 $\mu$m Spitzer temperature curve is consistent with the
expected $T_{eq}$ of XO-3b, the 3.6 $\mu$m temperature curve is higher and
does not intersect with the 4.5 $\mu$m temperature curve at any point in the
orbit. In fact, the 3.6 $\mu$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_{\rm 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 $\mu$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 ($\sim$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 pressures 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 $\sim 6$ hours. An internal energy source
flux of $7.21^{+0.11}_{-0.17}\times 10^{5}$ W/m2 is required to fit the 3.6
$\mu$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\sim 6\times 10^{3}$. Given the above
issue with the 3.6 $\mu$m data, these results should be taken lightly.

Figure 10: 1D atmospheric model at apoapse and periapse. 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.  Figure 11:
Published Spitzer/IRAC secondary eclipse depths are show along with our new
eclipse depths. Gray dashed lines show notional eclipse spectra if the star
and planet radiate as blackbodies with $T_{\rm eff,*}=6885$ K and planetary
dayside temperature of $T_{\rm day}=[2000,2342,2600]$ K. The deeper 3.6 $\mu$m
eclipse depth in comparison to the 4.5 $\mu$m eclipse depth rules out the
thermal inversion reported by Machalek et al. (2010): water vapour opacity
dictates that the 3.6 $\mu$m photons originate from deeper in the atmosphere
than those at 4.5 $\mu$m.

###  5.4 Theoretical Models of XO-3b’s Atmosphere

####  5.4.1 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_{\rm int}$) was set to 300 K a periapse, which sets the
adiabat for the self-consistent model in radiative-convective equilibrium. The
temperatures in the deep part of the atmosphere are expected to be
horizontally uniform. Therefore, the $T_{\rm int}$ value for apoapse was
iterated until the converged apoapse model fell upon the same deep adiabat,
yielding $T_{\rm int}$=520 K. The planet’s temperature should be homogenized
at depth, but, $T_{\rm{int}}$ can’t be thought of as constant in a simplified
1D modeling framework because the T-P profiles won’t lie on the 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_{\rm 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~{}\mu$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, FeSiO3 and
MgSiO3 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 $\mu$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 $\mu$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.

Figure 12: General circulation model (GCM) thermal phase curve predictions for
the 3.6$~{}\mu$m (left) and 4.5$~{}\mu$m (right) Spitzer band passes. The
coloured swath represent the 1$\sigma$ 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$~{}\mu$m
Spitzer phase curve and GCM predictions is likely due to issues with the
observations. At 4.5$~{}\mu$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.

Figure 13: Temperature (colorscale) and winds (arrows) at the $\approx$ 10
mbar and $\approx$ 120 mbar level of our 1$\times$ 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.

####  5.4.2 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$\times$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” rotating, 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_{\rm rot}=1.852\times 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 (Lewis et al., 2014).

Comparing GCM predictions to our Spitzer phase curves in Figure 12, we find
that the planetary 3.6 $\mu$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 $\mu$m observations. However, it is surprising that the 3.6 $\mu$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 $\mu$m vs. 4.5 $\mu$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 $\mathrm{CH_{4}}$ and CO that are strong absorbers in the
3.6 and 4.5 $\mu$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 $\mathrm{CH_{4}}$ in the
planet’s atmosphere throughout its orbit. Such a scenario would naturally
allow the 3.6 $\mu$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
$\mu$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 $\mu$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 $\mu$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., 2013b, 2014) and HD80606b (e.g.,
de Wit et al., 2016b), 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_{\rm 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 $\mu$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 $\mu$m phase curve. The discrepancy between the 3.6 $\mu$m phase curve
and the GCM prediction could also be due instrumental issues with the 3.6
$\mu$m channel. However, the GCM also under-predicts the robust 3.6 $\mu$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_{\rm eff}=6759\pm 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\pm 0.15$ $M_{\rm{\odot}}$ and
$R_{*}=1.407\pm 0.038R_{\rm{\odot}}$ and a planetary mass and radius of
$M_{p}=11.79\pm 0.98$ $M_{\rm J}$ and $R_{p}=1.295\pm 0.066R_{\rm 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$\sigma$ 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_{\rm 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$\sigma$ 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.

##  6 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\pm 0.98$ $M_{\rm{jup}}$ inflated hot
Jupiter with $R_{p}=1.295\pm 0.066R_{\rm{jup}}$ on a 3.2 day orbit with an
orbital eccentricity of $e=0.2769$. The full-orbit 3.6 $\mu$m and 4.5 $\mu$m
Spitzer/IRAC phase curves of XO-3b yield a secondary eclipse depths of 1770
$\pm$ 180 ppm and 1610 $\pm$ 70 ppm at 3.6 $\mu$m and 4.5 $\mu$m,
respectively. From the secondary eclipse portion of the 3.6 $\mu$m partial
phase curve of XO-3b obtained in 2010 (PI: P. Machalek, PID 60058), we
retrieve an eclipse depth of 1520 $\pm$ 130 ppm which agrees with our more
recent phase curve observations, but is 5$\sigma$ 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
$\mu$m phase curve. Nonetheless, we compare our reliable 4.5 $\mu$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 $\mu$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_{\rm wind}$ of $3.13^{+0.26}_{-0.83}$ km/s,
in agreement with the $\sim$ 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 $\mu$m and large discrepancies at 3.6 $\mu$m. While
the disagreement could be due to detector systematics spoiling the 3.6 $\mu$m
phase curve, it’s unlikely to be culprit for the difference in Spitzer-
measured and GCM-predicted 3.6 $\mu$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.

L.D. acknowledges support in part through the Technologies for Exo-Planetary
Science (TEPS) PhD Fellowship, and the Natural Sciences and Engineering
Research Council of Canada (NSERC)’s Postgraduate Scholarships-Doctoral
Fellowship. This work is based on archival data obtained with the Spitzer
Space Telescope, which is operated by the Jet Propulsion Laboratory,
California Institute of Technology under a contract with NASA. T.J.B.
acknowledges support from the McGill Space Institute Graduate Fellowship, the
Natural Sciences and Engineering Research Council of Canada’s Postgraduate
Scholarships-Doctoral Fellowship, and from the Fonds de recherche du Québec –
Nature et technologies through the Centre de recherche en astrophysique du
Québec.

##  Appendix A Tidal Heating

Time-dependent tidal distortion of a body leads to internal heating (Peale &
Cassen, 1978; Peale et al., 1979; Wisdom, 2004). Hence, eccentricity and
obliquity tides have been proposed as the missing energy source to explain the
anomalously large radius of some hot Jupiters. We estimate the rate of energy
dissipation using the formalism described in Levrard et al. (2007) which takes
into account the effect of synchronous rotation:

$\frac{dE_{\rm{tides}}}{dt}=2K\left[N_{a}(e)-\frac{N^{2}(e)}{\Omega(e)}\frac{2x^{2}}{1+x^{2}}\right]$
(A1)

where $\Omega(e)$, $N(e)$, and $N_{a}(e)$ are functions of the orbital
eccentricity and are defined as:

$\Omega(e)=\frac{1+3e^{2}+\frac{3}{8}e^{4}}{\left(1-e^{2}\right)^{9/2}},$ (A2)
$N(e)=\frac{1+\frac{15}{2}e^{2}+\frac{45}{8}e^{4}+\frac{5}{16}e^{6}}{\left(1-e^{2}\right)^{6}},$
(A3)
$N_{a}(e)=\frac{\left(1+\frac{31}{2}e^{2}+\frac{255}{8}e^{4}+\frac{185}{16}e^{6}+\frac{25}{64}e^{8}\right)}{(1-e^{2})^{15/2}}.$
(A4)

and where $x=\cos\epsilon$, $\epsilon$ is the planet’s obliquity (the angle
between the equatorial and orbital planes) and

$K=\frac{3}{2}\frac{k_{2}}{Q_{n}}\left(\frac{GM_{p}^{2}}{R_{p}}\right)\left(\frac{M_{*}}{M_{p}}\right)^{2}\left(\frac{R_{p}}{a}\right)^{6}n$
(A5)

where $k_{2}$ is the planet’s potential Love number of degree 2, $Q_{n}$ is
the planet’s annual tidal quality factor, $G$ is the gravitational constant,
and $n$ is the planet’s mean motion which is approximately
$\sqrt{GM_{*}/a^{3}}$. Assuming that XO-3b has a zero obliquity and Jupiter’s
tidal Love number $k_{2}=0.565\pm 0.006$ (Durante et al., 2020), we find that
a tidal quality factor of $Q\sim 6\times 10^{3}$ is required to explain the
excess flux of $1.3\times$ insolation inferred for the 3.6 $\mu$m phase curve.

## References

  * Adams & Laughlin (2018) Adams, A. D., & Laughlin, G. 2018, AJ, 156, 28, doi: 10.3847/1538-3881/aac437
  * Adams et al. (2019) Adams, A. D., Millholland, S., & Laughlin, G. P. 2019, AJ, 158, 108, doi: 10.3847/1538-3881/ab2b35
  * Adcroft et al. (2004) Adcroft, A., Campin, J.-M., Hill, C., & Marshall, J. 2004, Monthly Weather Review, 132, 2845, doi: 10.1175/MWR2823.1
  * Arcangeli et al. (2021) Arcangeli, J., Désert, J. M., Parmentier, V., Tsai, S. M., & Stevenson, K. B. 2021, A&A, 646, A94, doi: 10.1051/0004-6361/202038865
  * Arcangeli et al. (2019) Arcangeli, J., Désert, J.-M., Parmentier, V., et al. 2019, A&A, 625, A136, doi: 10.1051/0004-6361/201834891
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Beatty et al. (2019) Beatty, T. G., Marley, M. S., Gaudi, B. S., et al. 2019, AJ, 158, 166, doi: 10.3847/1538-3881/ab33fc
  * Bell & Cowan (2018) Bell, T. J., & Cowan, N. B. 2018, ApJ, 857, L20, doi: 10.3847/2041-8213/aabcc8
  * Bell et al. (2019) Bell, T. J., Zhang, M., Cubillos, P. E., et al. 2019, MNRAS, 489, 1995, doi: 10.1093/mnras/stz2018
  * Bell et al. (2021) Bell, T. J., Dang, L., Cowan, N. B., et al. 2021, MNRAS, 504, 3316, doi: 10.1093/mnras/stab1027
  * Bodenheimer et al. (2013) Bodenheimer, P., D’Angelo, G., Lissauer, J. J., Fortney, J. J., & Saumon, D. 2013, ApJ, 770, 120, doi: 10.1088/0004-637X/770/2/120
  * Bodenheimer et al. (2000) Bodenheimer, P., Hubickyj, O., & Lissauer, J. J. 2000, Icarus, 143, 2, doi: 10.1006/icar.1999.6246
  * Bordwell et al. (2018) Bordwell, B., Brown, B. P., & Oishi, J. S. 2018, ApJ, 854, 8, doi: 10.3847/1538-4357/aaa551
  * Challener et al. (2021) Challener, R. C., Harrington, J., Jenkins, J., et al. 2021, The Planetary Science Journal, 2, 9, doi: 10.3847/PSJ/abc954
  * Charbonneau et al. (2005) Charbonneau, D., Allen, L. E., Megeath, S. T., et al. 2005, ApJ, 626, 523, doi: 10.1086/429991
  * Cowan & Agol (2008) Cowan, N. B., & Agol, E. 2008, ApJ, 678, L129, doi: 10.1086/588553
  * Cowan & Agol (2011) —. 2011, ApJ, 726, 82, doi: 10.1088/0004-637X/726/2/82
  * Cowan et al. (2012) Cowan, N. B., Machalek, P., Croll, B., et al. 2012, ApJ, 747, 82, doi: 10.1088/0004-637X/747/1/82
  * Cubillos et al. (2017) Cubillos, P., Harrington, J., Loredo, T. J., et al. 2017, AJ, 153, 3, doi: 10.3847/1538-3881/153/1/3
  * Dang et al. (2018) Dang, L., Cowan, N. B., Schwartz, J. C., et al. 2018, Nature Astronomy, 2, 220, doi: 10.1038/s41550-017-0351-6
  * Daylan et al. (2021) Daylan, T., Günther, M. N., Mikal-Evans, T., et al. 2021, AJ, 161, 131, doi: 10.3847/1538-3881/abd8d2
  * de Wit et al. (2016a) de Wit, J., Lewis, N. K., Langton, J., et al. 2016a, ApJ, 820, L33, doi: 10.3847/2041-8205/820/2/L33
  * de Wit et al. (2016b) —. 2016b, ApJ, 820, L33, doi: 10.3847/2041-8205/820/2/L33
  * Deming et al. (2011) Deming, D., Knutson, H., Agol, E., et al. 2011, ApJ, 726, 95, doi: 10.1088/0004-637X/726/2/95
  * Deming et al. (2015) Deming, D., Knutson, H., Kammer, J., et al. 2015, ApJ, 805, 132, doi: 10.1088/0004-637X/805/2/132
  * Dobbs-Dixon & Cowan (2017) Dobbs-Dixon, I., & Cowan, N. B. 2017, ApJ, 851, L26, doi: 10.3847/2041-8213/aa9bec
  * Durante et al. (2020) Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572, doi: 10.1029/2019GL086572
  * Fabrycky et al. (2007) Fabrycky, D. C., Johnson, E. T., & Goodman, J. 2007, ApJ, 665, 754, doi: 10.1086/519075
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10, doi: 10.1086/422843
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306, doi: 10.1086/670067
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419, doi: 10.1086/528370
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gaia Collaboration et al. (2021) —. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Goody et al. (1989) Goody, R., West, R., Chen, L., & Crisp, D. 1989, J. Quant. Spec. Radiat. Transf., 42, 539, doi: 10.1016/0022-4073(89)90044-7
  * Hansen et al. (2014) Hansen, C. J., Schwartz, J. C., & Cowan, N. B. 2014, MNRAS, 444, 3632, doi: 10.1093/mnras/stu1699
  * Hébrard et al. (2008) Hébrard, G., Bouchy, F., Pont, F., et al. 2008, A&A, 488, 763, doi: 10.1051/0004-6361:200810056
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Ibgui & Burrows (2009) Ibgui, L., & Burrows, A. 2009, ApJ, 700, 1921, doi: 10.1088/0004-637X/700/2/1921
  * Ingalls et al. (2016) Ingalls, J. G., Krick, J. E., Carey, S. J., et al. 2016, AJ, 152, 44, doi: 10.3847/0004-6256/152/2/44
  * Iro & Deming (2010) Iro, N., & Deming, L. D. 2010, ApJ, 712, 218, doi: 10.1088/0004-637X/712/1/218
  * Johns-Krull et al. (2008) Johns-Krull, C. M., McCullough, P. R., Burke, C. J., et al. 2008, ApJ, 677, 657, doi: 10.1086/528950
  * Kataria et al. (2014) Kataria, T., Showman, A. P., Fortney, J. J., Marley, M. S., & Freedman, R. S. 2014, ApJ, 785, 92, doi: 10.1088/0004-637X/785/2/92
  * Kataria et al. (2013) Kataria, T., Showman, A. P., Lewis, N. K., et al. 2013, ApJ, 767, 76, doi: 10.1088/0004-637X/767/1/76
  * Keating & Cowan (2017) Keating, D., & Cowan, N. B. 2017, ApJ, 849, L5, doi: 10.3847/2041-8213/aa8b6b
  * Keating et al. (2019) Keating, D., Cowan, N. B., & Dang, L. 2019, Nature Astronomy, 3, 1092, doi: 10.1038/s41550-019-0859-z
  * Keating et al. (2020) Keating, D., Stevenson, K. B., Cowan, N. B., et al. 2020, AJ, 159, 225, doi: 10.3847/1538-3881/ab83f4
  * Kipping (2013) Kipping, D. M. 2013, MNRAS, 435, 2152, doi: 10.1093/mnras/stt1435
  * Knutson et al. (2012) Knutson, H. A., Lewis, N., Fortney, J. J., et al. 2012, ApJ, 754, 22, doi: 10.1088/0004-637X/754/1/22
  * Koll & Komacek (2018) Koll, D. D. B., & Komacek, T. D. 2018, ApJ, 853, 133, doi: 10.3847/1538-4357/aaa3de
  * Kreidberg (2015) Kreidberg, L. 2015, Publications of the Astronomical Society of the Pacific, 127, 1161, doi: 10.1086/683602
  * Kreidberg et al. (2018) Kreidberg, L., Line, M. R., Parmentier, V., et al. 2018, AJ, 156, 17, doi: 10.3847/1538-3881/aac3df
  * Langton & Laughlin (2008) Langton, J., & Laughlin, G. 2008, ApJ, 674, 1106, doi: 10.1086/523957
  * Lanotte et al. (2014) Lanotte, A. A., Gillon, M., Demory, B. O., et al. 2014, A&A, 572, A73, doi: 10.1051/0004-6361/201424373
  * Levrard et al. (2007) Levrard, B., Correia, A. C. M., Chabrier, G., et al. 2007, A&A, 462, L5, doi: 10.1051/0004-6361:20066487
  * Lewis et al. (2017) Lewis, N. K., Parmentier, V., Kataria, T., et al. 2017, arXiv e-prints, arXiv:1706.00466.  https://arxiv.org/abs/1706.00466
  * Lewis et al. (2014) Lewis, N. K., Showman, A. P., Fortney, J. J., Knutson, H. A., & Marley, M. S. 2014, ApJ, 795, 150, doi: 10.1088/0004-637X/795/2/150
  * Lewis et al. (2010) Lewis, N. K., Showman, A. P., Fortney, J. J., et al. 2010, ApJ, 720, 344, doi: 10.1088/0004-637X/720/1/344
  * Lewis et al. (2013a) Lewis, N. K., Knutson, H. A., Showman, A. P., et al. 2013a, ApJ, 766, 95, doi: 10.1088/0004-637X/766/2/95
  * Lewis et al. (2013b) —. 2013b, ApJ, 766, 95, doi: 10.1088/0004-637X/766/2/95
  * Lin et al. (1996) Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606, doi: 10.1038/380606a0
  * Liu et al. (2008) Liu, X., Burrows, A., & Ibgui, L. 2008, ApJ, 687, 1191, doi: 10.1086/592189
  * Luger et al. (2018) Luger, R., Kruse, E., Foreman-Mackey, D., Agol, E., & Saunders, N. 2018, AJ, 156, 99, doi: 10.3847/1538-3881/aad230
  * Machalek et al. (2010) Machalek, P., Greene, T., McCullough, P. R., et al. 2010, ApJ, 711, 111, doi: 10.1088/0004-637X/711/1/111
  * Mansfield et al. (2020) Mansfield, M., Bean, J. L., Stevenson, K. B., et al. 2020, ApJ, 888, L15, doi: 10.3847/2041-8213/ab5b09
  * Mardling & Lin (2002) Mardling, R. A., & Lin, D. N. C. 2002, ApJ, 573, 829, doi: 10.1086/340752
  * Marley & McKay (1999) Marley, M. S., & McKay, C. P. 1999, Icarus, 138, 268, doi: 10.1006/icar.1998.6071
  * Marley et al. (2002) Marley, M. S., Seager, S., Saumon, D., et al. 2002, ApJ, 568, 335, doi: 10.1086/338800
  * Maxted et al. (2013) Maxted, P. F. L., Anderson, D. R., Doyle, A. P., et al. 2013, MNRAS, 428, 2645, doi: 10.1093/mnras/sts231
  * May et al. (2021) May, E. M., Komacek, T. D., Stevenson, K. B., et al. 2021, arXiv e-prints, arXiv:2107.03349.  https://arxiv.org/abs/2107.03349
  * Mayorga et al. (2021) Mayorga, L. C., Robinson, T. D., Marley, M. S., May, E. M., & Stevenson, K. B. 2021, arXiv e-prints, arXiv:2105.08009.  https://arxiv.org/abs/2105.08009
  * Millholland (2019) Millholland, S. 2019, ApJ, 886, 72, doi: 10.3847/1538-4357/ab4c3f
  * Naoz et al. (2013) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, MNRAS, 431, 2155, doi: 10.1093/mnras/stt302
  * Parmentier et al. (2018) Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, A&A, 617, A110, doi: 10.1051/0004-6361/201833059
  * Peale & Cassen (1978) Peale, S. J., & Cassen, P. 1978, Icarus, 36, 245, doi: 10.1016/0019-1035(78)90109-4
  * Peale et al. (1979) Peale, S. J., Cassen, P., & Reynolds, R. T. 1979, Science, 203, 892, doi: 10.1126/science.203.4383.892
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38, doi: 10.1051/0004-6361/201937381
  * Rasio & Ford (1996) Rasio, F. A., & Ford, E. B. 1996, Science, 274, 954, doi: 10.1126/science.274.5289.954
  * Schwarz et al. (1978) Schwarz, G., et al. 1978, The annals of statistics, 6, 461
  * Showman et al. (2009a) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009a, ApJ, 699, 564, doi: 10.1088/0004-637X/699/1/564
  * Showman et al. (2009b) —. 2009b, ApJ, 699, 564, doi: 10.1088/0004-637X/699/1/564
  * Showman & Guillot (2002) Showman, A. P., & Guillot, T. 2002, A&A, 385, 166, doi: 10.1051/0004-6361:20020101
  * Socrates (2013) Socrates, A. 2013, arXiv e-prints, arXiv:1304.4121.  https://arxiv.org/abs/1304.4121
  * Soubiran et al. (2020) Soubiran, C., Le Campion, J. F., Brouillet, N., & Chemin, L. 2020, VizieR Online Data Catalog, B/pastel
  * Spiegel et al. (2011a) Spiegel, D. S., Burrows, A., & Milsom, J. A. 2011a, ApJ, 727, 57, doi: 10.1088/0004-637X/727/1/57
  * Spiegel et al. (2011b) —. 2011b, ApJ, 727, 57, doi: 10.1088/0004-637X/727/1/57
  * Stassun et al. (2017) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136, doi: 10.3847/1538-3881/aa5df3
  * Stassun & Torres (2018) Stassun, K. G., & Torres, G. 2018, ApJ, 862, 61, doi: 10.3847/1538-4357/aacafc
  * Steinrueck et al. (2019) Steinrueck, M. E., Parmentier, V., Showman, A. P., Lothringer, J. D., & Lupu, R. E. 2019, ApJ, 880, 14, doi: 10.3847/1538-4357/ab2598
  * Stevenson et al. (2012) Stevenson, K. B., Harrington, J., Fortney, J. J., et al. 2012, ApJ, 754, 136, doi: 10.1088/0004-637X/754/2/136
  * Stevenson et al. (2014) Stevenson, K. B., Désert, J.-M., Line, M. R., et al. 2014, Science, 346, 838, doi: 10.1126/science.1256758
  * Stevenson et al. (2017) Stevenson, K. B., Line, M. R., Bean, J. L., et al. 2017, AJ, 153, 68, doi: 10.3847/1538-3881/153/2/68
  * Thorngren et al. (2019) Thorngren, D., Gao, P., & Fortney, J. J. 2019, ApJ, 884, L6, doi: 10.3847/2041-8213/ab43d0
  * Thorngren & Fortney (2018) Thorngren, D. P., & Fortney, J. J. 2018, AJ, 155, 214, doi: 10.3847/1538-3881/aaba13
  * Thorngren et al. (2016) Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64, doi: 10.3847/0004-637X/831/1/64
  * Tremblin et al. (2017) Tremblin, P., Chabrier, G., Mayne, N. J., et al. 2017, ApJ, 841, 30, doi: 10.3847/1538-4357/aa6e57
  * Turner et al. (2017) Turner, J. D., Leiter, R. M., Biddle, L. I., et al. 2017, MNRAS, 472, 3871, doi: 10.1093/mnras/stx2221
  * Visscher (2012) Visscher, C. 2012, ApJ, 757, 5, doi: 10.1088/0004-637X/757/1/5
  * von Essen et al. (2021) von Essen, C., Mallonn, M., Piette, A., et al. 2021, A&A, 648, A71, doi: 10.1051/0004-6361/202038524
  * Weidenschilling & Marzari (1996) Weidenschilling, S. J., & Marzari, F. 1996, Nature, 384, 619, doi: 10.1038/384619a0
  * Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1, doi: 10.1086/422992
  * Winn et al. (2008) Winn, J. N., Holman, M. J., Torres, G., et al. 2008, ApJ, 683, 1076, doi: 10.1086/589737
  * Winn et al. (2009) Winn, J. N., Johnson, J. A., Fabrycky, D., et al. 2009, ApJ, 700, 302, doi: 10.1088/0004-637X/700/1/302
  * Wisdom (2004) Wisdom, J. 2004, AJ, 128, 484, doi: 10.1086/421360
  * Wong et al. (2014) Wong, I., Knutson, H. A., Cowan, N. B., et al. 2014, ApJ, 794, 134, doi: 10.1088/0004-637X/794/2/134
  * Wong et al. (2015) Wong, I., Knutson, H. A., Lewis, N. K., et al. 2015, ApJ, 811, 122, doi: 10.1088/0004-637X/811/2/122
  * Wong et al. (2016) Wong, I., Knutson, H. A., Kataria, T., et al. 2016, ApJ, 823, 122, doi: 10.3847/0004-637X/823/2/122
  * Wong et al. (2020) Wong, I., Benneke, B., Shporer, A., et al. 2020, AJ, 159, 104, doi: 10.3847/1538-3881/ab6d6e
  * Wu & Lithwick (2011) Wu, Y., & Lithwick, Y. 2011, ApJ, 735, 109, doi: 10.1088/0004-637X/735/2/109
  * Youdin & Mitchell (2010) Youdin, A. N., & Mitchell, J. L. 2010, ApJ, 721, 1113, doi: 10.1088/0004-637X/721/2/1113
  * Zellem et al. (2014) Zellem, R. T., Lewis, N. K., Knutson, H. A., et al. 2014, ApJ, 790, 53, doi: 10.1088/0004-637X/790/1/53
  * Zhang et al. (2018) Zhang, M., Knutson, H. A., Kataria, T., et al. 2018, AJ, 155, 83, doi: 10.3847/1538-3881/aaa458