# Global Climate and Atmospheric Composition of the Ultra-hot Jupiter WASP-103b from HST and Spitzer Phase Curve Observations

Laura Kreidberg, Michael R. Line, Vivien Parmentier, Kevin B. Stevenson, Tom Louden, Mickäel Bonnefoy, Jacqueline K. Faherty, Gregory W. Henry, Michael H. Williamson, Keivan Stassun, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Laura Kreidberg, Michael R. Line, Vivien Parmentier, Kevin B. Stevenson, Tom Louden, et al. Global Climate and Atmospheric Composition of the Ultra-hot Jupiter WASP-103b from HST and Spitzer Phase Curve Observations. The Astronomical Journal, 2018, 156, 10.3847/1538-3881/aac3df. insu-03693548

HAL Id: insu-03693548

https://insu.hal.science/insu-03693548v1

Submitted on 22 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Global Climate and Atmospheric Composition of the Ultra-hot Jupiter WASP-103b from *HST* and *Spitzer* Phase Curve Observations

Laura Kreidberg<sup>1,2</sup>, Michael R. Line<sup>3</sup>, Vivien Parmentier<sup>4</sup>, Kevin B. Stevenson<sup>5</sup>, Tom Louden<sup>6</sup>, Mickäel Bonnefoy<sup>7</sup>, Jacqueline K. Faherty<sup>8</sup>, Gregory W. Henry<sup>9</sup>, Michael H. Williamson<sup>9</sup>, Keivan Stassun<sup>10</sup>, Thomas G. Beatty<sup>11,12</sup>, Jacob L. Bean<sup>13</sup>

<sup>1</sup> Haryard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

<sup>2</sup> Harvard Society of Fellows, 78 Mount Auburn Street, Cambridge, MA 02138, USA

<sup>3</sup> School of Earth & Space Exploration, Arizona State University, Tempe AZ 85287, USA

<sup>4</sup> Aix Marseille Univ, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France

<sup>5</sup> Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

<sup>6</sup> Department of Physics, University of Warwick, Coventry, CV4 7AL, UK

<sup>7</sup> Univ. Grenoble Alpes, IPAG, F-38000 Grenoble, France. CNRS, IPAG, F-38000 Grenoble, France

<sup>8</sup> American Museum of Natural History, Department of Astrophysics, Central Park West at 79th Street, New York, NY 10034, USA

<sup>9</sup> Center of Excellence in Information Systems, Tennessee State University, Nashville, TN 37209, USA

<sup>10</sup> Vanderbilt University, Dept. of Physics and Astronomy, 6301 Stevenson Center Ln, Nashville TN, 37235, USA

<sup>11</sup> Department of Astronomy & Astrophysics, The Pennsylvania State University, 525 Davey Lab, University Park, PA 16802, USA

<sup>12</sup> Center for Exoplanets and Habitable Worlds, The Pennsylvania State University, 525 Davey Lab, University Park, PA 16802, USA

<sup>13</sup> Department of Astronomy & Astrophysics, University of Chicago, 5640 S. Ellis Avenue, Chicago, IL 60637, USA

<sup>14</sup> Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA

<sup>15</sup> Department of Planetary Sciences and Lunar and Planetary Laboratory, University of Arizona, Tucson, Arizona 85721, USA

<sup>16</sup> Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

*Received 2018 February 15; revised 2018 May 2; accepted 2018 May 8; published 2018 June 18*

## Abstract

We present thermal phase curve measurements for the hot Jupiter WASP-103b observed with *Hubble*/WFC3 and *Spitzer*/IRAC. The phase curves have large amplitudes and negligible hotspot offsets, indicative of poor heat redistribution to the nightside. We fit the phase variation with a range of climate maps and find that a spherical harmonics model generally provides the best fit. The phase-resolved spectra are consistent with blackbodies in the WFC3 bandpass, with brightness temperatures ranging from 1880 ± 40 K on the nightside to 2930 ± 40 K on the dayside. The dayside spectrum has a significantly higher brightness temperature in the *Spitzer* bands, likely due to CO emission and a thermal inversion. The inversion is not present on the nightside. We retrieved the atmospheric composition and found that it is moderately metal-enriched ([M/H] = 23±59/13 × solar) and the carbon-to-oxygen ratio is below 0.9 at 3σ confidence. In contrast to cooler hot Jupiters, we do not detect spectral features from water, which we attribute to partial H₂O dissociation. We compare the phase curves to 3D general circulation models and find that magnetic drag effects are needed to match the data. We also compare the WASP-103b spectra to brown dwarfs and young, directly imaged companions. We find that these objects have significantly larger water features, indicating that surface gravity and irradiation environment play an important role in shaping the spectra of hot Jupiters. These results highlight the 3D structure of exoplanet atmospheres and illustrate the importance of phase curve observations for understanding their complex chemistry and physics.

**Key words:** planets and satellites: atmospheres – planets and satellites: gaseous planets – planets and satellites: individual (WASP-103b)

## 1. Introduction

Planets are round, rotating, and irradiated on one hemisphere at a time—all of which contribute to rich spatial structure in their climate and atmospheric composition. Short-period, tidally locked planets are an extreme example, with one hot, continuously illuminated side. This asymmetry is expected to produce large gradients in temperature, chemistry, and cloud coverage with longitude (Showman et al. 2009; Kataria et al. 2016; Parmentier et al. 2016), and provides an opportunity to learn about atmospheric dynamics in a very different regime from the planets of our solar system.

Exoplanets are so distant that they are generally not spatially resolved from their host stars, but it is still possible to reveal inhomogeneities in their atmospheres by observing the total system flux. One approach is to measure a phase curve, which consists of continuous monitoring of the planet-to-star flux ratio over a complete orbital revolution of the planet. This

observation is sensitive to different longitudes at each orbital phase of the planet. The first phase curve of an exoplanet was observed with *Spitzer* for the hot Jupiter *ν* Andromedae b by Harrington et al. (2006), followed by additional *Spitzer* observations for about a dozen more systems (cataloged in Parmentier & Crossfield 2017). These observations revealed large day-night temperature contrasts (in excess of 300 K), and eastward-shifted peak brightness due to heat circulation, as predicted by 3D models (Showman & Guillot 2002). These infrared measurements were complemented by optical phase curves from *Kepler* that show evidence for reflected light from patchy and possibly variable dayside clouds with a range of compositions (Borucki et al. 2009; Demory et al. 2013; Hu et al. 2015; Armstrong et al. 2016; Parmentier et al. 2016). A *spectroscopic* phase curve was observed for WASP-43b with *Hubble*/Wide Field Camera 3 (HST/WFC3) in the near-infrared, which provided the first phase-resolved measurements
of an exoplanet's water abundance and thermal structure (Stevenson et al. 2014c, 2017).

In this paper, we present spectroscopic phase curve observations of the hot Jupiter WASP-103b, measured with HST/WFC3 and Spitzer/IRAC. This planet is an ideal target for phase curve observations, with an orbital period of just 22 hr and an equilibrium temperature of 2500 K. WASP-103b is slightly larger than Jupiter, with a mass and radius of $1.49 \pm 0.09 M_{\text{Jup}}$ and $1.53_{-0.07}^{+0.05} R_{\text{Jup}}$, respectively. The host star is a main-sequence F8 dwarf with an effective temperature of 6110 ± 160 K (Gillon et al. 2014). Previous observations of WASP-103b's atmosphere revealed a blackbody-like dayside emission spectrum, with possible evidence for a K<sub>S</sub>-band emission feature (Cartier et al. 2017; Delrez et al. 2018). The optical transmission spectrum shows evidence for sodium and potassium absorption features that are consistent with expectations for a cloud-free atmosphere (Lendl et al. 2017).

WASP-103b is an archetype of the class of ultra-hot Jupiters with orbital periods of about one day and dayside temperatures typically >2000 K. These very hot planets were initially predicted to have inverted temperature pressure profiles due to strong optical absorption by TiO/VO in the upper atmospheres (Hubeny et al. 2003; Fortney et al. 2008); however, observations of their emission spectra have been inconclusive with regard to their thermal structure and composition. In the near-infrared, where water is the dominant absorber, some spectra show water absorption features, some show emission features, and some are consistent with blackbody models (Madhusudhan et al. 2011; Crossfield et al. 2012; Stevenson et al. 2014a; Haynes et al. 2015; Evans et al. 2016; Beatty et al. 2017a, 2017b; Sheppard et al. 2017; Arcangeli et al. 2018; Mansfield et al. 2018). A variety of explanations have been proposed for these results, including low metallicity or high carbon-to-oxygen compositions, dayside clouds, and finely tuned isothermal temperature pressure profiles. Recently, Arcangeli et al. (2018) and Lothringer et al. (2018) showed that water dissociation and H<sup>-</sup> opacity on the hot dayside play an important role in the atmospheres of these ultra-hot planets and may be responsible for some of the blackbody-like near-IR spectra. In this work, we contextualize these results by investigating the global thermal structure and composition of the ultra-hot Jupiter WASP-103b.

The structure of the paper is as follows. In Section 2, we describe the observations and data reduction. Section 3 details the models fit to the phase curves. In Section 4, we discuss results, including the phase curve amplitudes and hotspot offsets, the phase-resolved spectra, estimates of the planet's climate, and the transmission spectrum. In Sections 6 and 7, we compare the observations to general circulation model (GCM) predictions and spectra from similar temperature stars and directly imaged companions. We offer our conclusions in Section 8.

## 2. Observations and Data Reduction

We observed two full-orbit phase curves of WASP-103b with HST/WFC3 and one each with Spitzer/IRAC at 3.6 and 4.5 µm (from HST Program 14050 and Spitzer Program 11099, PI: L. Kreidberg). We also reduced two HST/WFC3 secondary eclipse observations of WASP-103b from HST Program 13660 (PI: M. Zhao).

### 2.1. HST/WFC3

The HST phase curve observations consisted of two visits on 2015 February 26–27 and August 2–3. Each visit was 15 orbits in duration and spanned 23 hr. The last half of orbit 15 in each visit was used for a gyro bias update and produced no usable science data. We took a direct image of the star with the F126N filter at the beginning of each orbit to determine the wavelength solution zero-point. The remainder of the orbit consisted of time-series spectroscopy with the G141 grism (1.1–1.7 µm) and the 256 × 256 pixel subarray. We used the SPARS10/NSAMP = 15 read-out mode, which has an exposure time of 103 s. To optimize the duty cycle of the observations, we used the spatial scan observing mode with a scan rate of 0.03 arcsec s<sup>-1</sup>, alternating between forward and backward scanning on the detector. The scan height was 25 pixels and the peak counts were $3.5 \times 10^4$ photoelectrons per pixel. We collected a total of 18 spatial scan exposures per orbit. The two eclipse observations from Program 13660 had a similar observing setup (described in detail in Cartier et al. 2017).

We reduced the data from both programs using a custom pipeline developed for past analyses of WFC3 data (for details, see Kreidberg et al. 2014a, 2014b, 2015). Briefly, we use the optimal extraction algorithm of Horne (1986) to separately extract each up-the-ramp sample (or “stripe”). The stripes are then summed to create the final spectrum. For each stripe, the extraction window is 24 pixels high and centered on the stripe midpoint. We estimate the background from the median of a region of the detector that is uncontaminated by the target spectrum (rows 5–50). The typical background counts are low (10–15 photoelectrons per pixel, roughly 0.03% of the peak counts from the target star). We note that the extracted spectrum includes flux from a nearby companion star, which is separated from WASP-103b by less than two pixels (0.12; Wöllert & Brandner 2015). We account for this contamination later in the analysis.

### 2.2. Spitzer

We also obtained Spitzer/IRAC observations with 3.6 and 4.5 µm photometric filters (referred to as Channel 1 and Channel 2, respectively). The observations had the following setup: each phase curve observation consisted of 30 hr of time series photometry, beginning three hours prior to one secondary eclipse and ending three hours after a second eclipse. We read out the full array and used 12 s exposures to maximize the duty cycle without saturating the detector. To minimize the intrapixel effect (variations in flux caused by imprecise pointing), we did not dither and also used PCRS peak-up<sup>17</sup> to improve the pointing accuracy. We began each observation with a 30 min position settling period, followed by three Astronomical Observation Requests (AORs) of equal duration. At the beginning of each AOR, the telescope was repointed to position the target in the “sweet spot” of the detector, where the response is fairly uniform over the pixel.

The data were reduced with the POET pipeline (Stevenson et al. 2012b; Cubillos et al. 2013). The pipeline starts by identifying and flagging bad pixels using a double-iteration four-sigma outlier rejection routine along the time axis. This is followed by performing 2D Gaussian centroiding on each frame, which has been shown to provide the most

<footnote n="17"><sup>17</sup> http://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/pcrsobs.shtml</footnote>
<table><caption>**Table 1**
Photometric Observations of WASP-103</caption><thead><tr><th>Observing Season</th><th>N<sub>obs</sub></th><th>Date Range (HJD-2,400,000)</th><th>Sigma (mag)</th><th>Seasonal Mean (mag)</th></tr></thead><tbody><tr><td>2014</td><td>59</td><td>56722-56972</td><td>0.0057</td><td>0.9546 ± 0.0007</td></tr><tr><td>2015</td><td>73</td><td>57028-57335</td><td>0.0062</td><td>0.9549 ± 0.0007</td></tr><tr><td>2016</td><td>26</td><td>57385-57451</td><td>0.0055</td><td>0.9485 ± 0.0011</td></tr></tbody></table>

precise centers for *Spitzer* data (Lust et al. 2014). The target remains centered near the sweet spot for the entire AOR in each observation, with a maximum drift of 0.1 pixels. Next, POET uses sub-pixel (5× interpolated) aperture photometry (Harrington et al. 2007) to subtract the background and sum the flux within a specified radius. Chosen from a grid of apertures between two and four pixels, we find that an aperture size of 2.75 pixels minimizes the residual noise in the light curve fits. For the background, we use an annulus with inner and outer radii of 7 and 15 pixels, respectively. The contaminating flux from the nearby star is within the same pixel as the target, so we included it in the photometry and corrected it in the light curve fits. A similarly strategy has been applied to successfully analyze dozens of *Spitzer* data sets (e.g., Stevenson et al. 2010, 2012a, 2012c, 2014a, 2014b, 2016, 2017; Campo et al. 2011; Nymeyer et al. 2011; Blecic et al. 2013, 2014; Cubillos et al. 2013; Diamond-Lowe et al. 2014).

## 2.3. Photometric Monitoring

To assess how stellar activity might impact the phase curve observations, we monitored WASP-103's photometric variability over 158 nights during 2014–2016 with the Tennessee State University Celestron 14 inch (C14) automated imaging telescope (AIT), located at Fairborn Observatory in southern Arizona (Henry 1999). The observations of WASP-103 were made in the Cousins R passband with an SBIG STL-1001E CCD camera. Each observation consisted of 4–10 consecutive exposures on WASP-103 along with several dozen comparison stars in the same field. The individual consecutive frames were co-added and reduced to differential magnitudes (i.e., WASP-103 minus the mean brightness of the six best comparison stars). The nightly observations were corrected for bias, flat-fielding, and differential atmospheric extinction. For each season, we determined extinction corrections with a linear least-squares fit to nightly differential magnitude as a function of airmass.

The photometric analyses are summarized for each observing season in Table 1. The standard deviations of a single observation with respect to the corresponding seasonal means are given in column 4; the mean of the three standard deviations is 0.0058 mag, suggesting there is little night-to-night variation in WASP-103. The three seasonal mean brightness values given in column 5 scatter about their grand mean with a standard deviation of 0.0036 mag, but we note that the most discrepant mean is from the third season, for which we have only partial coverage. Therefore, our results do not completely rule out low-level, year-to-year variability of <0.001 mag.

To maximize the possibility of detecting WASP-103's rotation, we normalized the photometry such that each observing season has the same mean, thereby removing any long-term variability in WASP-103 and/or the comparison stars (Figure 1, top panel). To estimate the stellar rotation

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Top: the normalized nightly Cousins R band photometric data set for WASP-103, acquired with the C14 automated imaging telescope at Fairborn Observatory. Vertical dashed lines denote separate observing seasons. Gaps are due to target visibility and the Arizona monsoon season (July–September). Middle: the frequency spectrum of the normalized data set suggests low-amplitude variability with a period of 6.814 days. Bottom: the normalized data set phased to the 6.814 day period, which we interpret as rotational modulation of a star spot or spots. A least-squares sine fit to the 6.814 day rotation period gives a peak-to-peak amplitude of just 0.005 mag.</figcaption></figure>

period, we performed a periodogram analysis of the normalized data set, based on least-squares fitting of sine curves. The resulting frequency spectrum and the phase curve computed with the best period are shown in the middle and lower panels of Figure 1, respectively. The best-fit period is 6.814 days, which agrees closely with the estimated stellar rotation period of 6.855 days (based on the projected stellar rotation velocity and stellar radius reported in Gillon et al. 2014). There are two nearby peaks in the periodogram (panel b of Figure 1) that are one-year aliases of each other, and we chose the peak that better matches the stellar rotation velocity. The peak-to-peak variability amplitude is 0.005 mag. Based on the formalism in Zellem et al. (2017), we calculate that this variability will bias the measured eclipse depth by ~10 parts per million (ppm) from epoch to epoch, which is well below the photon-limited precision of our measurements.

## 3. HST and *Spitzer* Light Curve Fits

We fit a two-component model to the light curves. One component models the astrophysical signal (the planet's thermal phase variation and transit), and the other component models the systematic noise introduced by time-dependent changes in instrument performance. For each light curve, we fit the physical and systematic components simultaneously, such
that the total observed flux as a function of time is given by $F(t) = F_{\text{physical}}(t) \times F_{\text{sys}}(t)$. For the HST data, where we observed two phase curves and two additional eclipses, we constrain the physical parameters to be the same for all visits, but allow some of the systematics parameters to vary (for more details, see Section 3.2.1). We fit the WFC3 band-integrated “white” light curve, as well as spectroscopic light curves created from 10 wavelength bins uniformly spaced at 0.05 µm intervals between 1.15 and 1.65 µm.

### 3.1. Astrophysical Signal

We assume the measured astrophysical signal $F_{\text{physical}}$ has the following form:

$$F_{\text{physical}}(\lambda, t) = T(\lambda, t) + c(\lambda, t) \times F_p/F_s(\lambda, t), \quad (1)$$

where $\lambda$ is wavelength, $T(\lambda, t)$ is the transit model (the fraction of the stellar disk that is visible at time $t$), $F_p/F_s(\lambda, t)$ is the disk-integrated planet-to-star flux, and $c$ is a correction factor for companion star dilution and the planet’s tidal distortion.

We calculated the transit model $T(t)$ with the `batman` package (Kreidberg 2015). Many of the physical parameters are tightly constrained by Southworth et al. (2015), so we fixed the orbital period, time of inferior conjunction, orbital inclination, and ratio of semimajor axis to stellar radius to the previously published values ($P = 0.925545613$ day, $t_0 = 2456836.2964455$ BJD$_{TDB}$, $i = 87.3^\circ$, and $a/R_s = 2.999$). As a test, we fit for these parameters with the `Spitzer` Channel 2 light curve, which has the best phase coverage and least systematic noise of the three data sets. We found that the transit parameters are consistent with the Southworth et al. (2015) results, so we proceeded with the remainder of the analysis holding those parameters fixed. The free parameters for the transit model were a wavelength-dependent transit depth $r_p(\lambda)$ and linear limb-darkening parameter $u(\lambda)$. More complex limb-darkening laws with additional free parameters were not merited according to the Bayesian Information Criterion (BIC). We initialized our MCMC chains on the least-squares best-fit parameters.

We modeled the planet-to-star flux $F_p/F_s$ in two different ways. First, we fit a sinusoid with a period equal to the planet’s orbital period. The free parameters were the sine curve amplitude and phase offset. For the second approach, we used the `spiderman` package (Louden & Kreidberg 2017) to model $F_p/F_s$. This package allows users to input a climate map (temperature or brightness as a function of latitude and longitude), and generate the corresponding flux ratio for an observation at time $t$. In our fit, we calculated the stellar flux with a NextGen model (Allard et al. 2012) interpolated to an effective temperature of 6110 K (Gillon et al. 2014), solar metallicity, and log g of 4.2 (in cgs units). For the planet flux, we tested three different maps: a two-temperature map, with a uniform dayside temperature $T_d$ and a uniform nightside temperature $T_n$; a spherical harmonics map of degree two (with four free parameters); and the physically motivated kinematic model from Zhang & Showman (2017), which has just three free parameters (the nightside temperature $T_n$, the change in temperature from day-to-night side $\Delta_T$, and the ratio of radiative to advective timescales $\xi$). In all cases, we assumed that the planet is tidally locked, such that each orbital revolution corresponds to one complete rotation on its spin axis.

We scaled the planet-to-star flux by a correction factor $c$ to account for dilution from the companion star and ellipsoidal

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Projected area of the planet as a function of orbital phase, normalized to unity at phase zero. The area variation was predicted analytically, using the model from Leconte et al. (2011a).</figcaption></figure>

variability due to the planet’s tidal distortion. The correction factor took the form:

$$c(\lambda, t) = [1 + \alpha(\lambda)]A(t), \quad (2)$$

where $\alpha(\lambda)$ is the additional fractional flux from the companion star and $A(t)$ is the sky-projected area of the planet. We estimated $\alpha(\lambda)$ based on the best fit spectral energy distribution from Cartier et al. (2017). The companion star contribution ranges from 10% at 1.1 µm to 16% at 4.5 µm. The uncertainty on the companion star flux contribution to the total system flux is less than 1%, which introduces negligible error in the estimated planet-to-star flux compared to the photon noise. We calculated $A(t)$ using the analytic formula from Equation B.9 of Leconte et al. (2011a), which computes the projected area of a triaxial ellipsoid. We estimated the ellipsoid properties using Table B.3 of Leconte et al. (2011b), assuming the planet radius is 1.5 $R_{\text{Jup}}$ and age is 5 Gyr. The predicted ellipsoidal variability is shown in Figure 2. At quadrature, the projected area is 8% larger than at phase zero (mid-transit). Using the analytic expression from Loeb & Gaudi (2003), we estimated the effect of Doppler beaming and found that it contributes less than 10 ppm to the measured flux.

### 3.2. Systematics

Both the HST and Spitzer phase curves have systematic noise caused by variations in the sensitivity of the instrument over time. For the HST/WFC3 data, the dominant systematic is an orbit-long exponential trend due to charge traps filling up over successive exposures (Long et al. 2015; Zhou et al. 2017). For Spitzer, the primary source of noise is the intrapixel sensitivity effect. The detector’s pixels do not have uniform sensitivity, so slight changes in telescope pointing cause the recorded flux to vary. In Figure 3, we show the raw light curves before systematic noise was removed. The systematics have comparable amplitude to the thermal phase variation signal, so they must be carefully corrected to recover the underlying planet-to-star flux.

#### 3.2.1. HST Systematics

We fit the WFC3 systematics using an analytic model of the form:

$$F_{\text{sys}}(t) = (c S(t) + v_1 t_v + v_2 t_v^2)(1 - \exp(-a t_{\text{orb}} - b)), \quad (3)$$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Raw light curves (points) for WASP-103b observed with HST/WFC3 (top four panels) and Spitzer/IRAC (bottom two panels). The black lines show the best fit models, which include the astrophysical signal and instrument systematics. The gray lines indicate the contribution from the instrument systematics alone (which would be observed for a source with constant brightness and no planet). For visual clarity, we corrected the HST data for the upstream-downstream effect and zoomed in on the phase variation, so the transits are not displayed in the panel.</figcaption></figure>

where $t_v$ is time elapsed since the first exposure in a visit and $t_{orb}$ is time since the first exposure in an orbit. Here, $S(t)$ is a scale factor equal to 1 for exposures with spatial scanning in the forward direction and $s$ for reverse scans, to account for the upstream-downstream effect (McCullough & MacKenty 2012). The orbit-long ramp parameters are consistent for all the visits, so we constrained $a, b$, and $s$ to have the same value for all visits in the final fit. The visit-long trends differ from visit to visit, so $c, v_1$, and $v_2$ were allowed to vary between visits. We fixed $v_2$ to zero for the two secondary eclipse observations from Program 13360 because the visit-long trend for shorter observations is fit well by a linear slope.

Some segments of the data exhibit stronger systematics than others, so we exclude these data in our final analysis. We drop the first orbit from every visit and the first exposure from every orbit (following common practice; see, e.g., Kreidberg et al. 2014b). We also discard exposures from the last half of orbit 15 from the phase curve observations, which were taken in staring mode to enable a gyro bias update. Because we observed two phase curves, we have complete orbital phase coverage of the planet despite discarding some data.

### 3.2.2. Spitzer Systematics

Warm Spitzer's primary systematic is intrapixel sensitivity variation, where the photometry depends on the precise location of the stellar center within its pixel. We fit this systematic using the Bilinearly Interpolated Subpixel Sensitivity (BLISS) mapping technique (Stevenson et al. 2012b). BLISS provides a flexible, non-analytic means to effectively weight the target flux by the spatial sensitivity variations within a pixel, while simultaneously fitting for other systematics and the physical parameters of the system. As demonstrated by Ingalls et al. (2016), the POET pipeline with BLISS mapping can accurately model simulated Spitzer light curves with known physical parameters and produce reliable results.

The BLISS sensitivity map is determined by bilinear interpolation over a grid of knots centered on the stellar flux. Each knot's sensitivity is calculated from the residuals to the light curve fit: the higher the flux values for data points near a given knot, the higher the detector sensitivity is at that position. To avoid overfitting, we chose the grid scale such that bilinear interpolation performed better than nearest-neighbor interpolation. For the 3.6 µm data, the grid scale was 0.008 pixel (0.0098 arcsec) in both $x$ and $y$. For 4.5 µm, the scale was 0.022 pixel (0.027 arcsec). In addition to the intrapixel sensitivity variation, we fit the data for a linear trend in time. We tested a quadratic trend, but did not find significant evidence for the additional model complexity based on the BIC.

### 3.3. Best Fits and Uncertainties

To determine the best fits, we performed a least-squares $\chi^2$ minimization for each wavelength and model. For a subset of these cases where we wish to calculate 68% confidence intervals, we also performed a Markov chain Monte Carlo (MCMC) analysis to estimate parameter uncertainties. These include the transit fits and the sine curve fit to the broadband phase curves. We used emcee (Foreman-Mackey et al. 2013) to fit the HST/WFC3 light curves and differential evolution Monte Carlo for the Spitzer fits (Braak 2006). We ran the MCMC until convergence according to the Gelman-Rubin statistic. We initialized the MCMC chain on the best-fit parameters and discarded the first 10% of the chain as burn-in. MCMC techniques only produce robust uncertainties when the noise is normally distributed and white, so to account for correlated noise in the 3.6 µm light curve (described in Section 3.4), we fit the wavelet model from Carter & Winn (2009) simultaneously with the other model parameters. We used a Haar wavelet and let the power spectral density of the red noise vary, following Diamond-Lowe et al. (2014). In our final fit, the noise power spectrum $1/f^\gamma$ had $\gamma = 1.1 \pm 0.1$ (implying an equal amount of white noise and correlated noise).

### 3.4. Goodness of Fit

We performed several tests of the quality of the light curve fits. First, we predicted the level of scatter in the light curves based on photon noise alone, then compared this value to the root-mean-square (rms) of the fit residuals. For the spherical harmonics fit to the phase variation, the Spitzer 4.5 µm light curve rms reaches the expected photon noise limit (637 versus 640 ppm). The 3.6 µm light curve has significantly larger rms (767 versus 470 ppm), due to time-correlated red noise (discussed below). The expected photon-limited rms for the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Root-mean-square variability in the light curves as a function of bin size (black lines) compared to the expected rms from photon noise (colored lines). The central wavelength of the light curve is indicated in the upper right corner of each panel. With the exception of the *Spitzer* 3.6 µm channel, the rms for the light curves bins down in agreement with predictions from the photon noise.</figcaption></figure>

WFC3 spectroscopic light curves range from 430–530 ppm, and the measured rms was typically within 5% of expectations for all spectroscopic channels. For the WFC3 band-integrated white light curve, the rms was slightly larger than predicted (172 versus 122 ppm). There are a number of possible origins for this discrepancy, including imperfect background subtraction, variation in the position of the spatial scan on the detector, and loss of flux outside the extraction window. In addition, the amplitude of the phase variation increases by 50% over the WFC3 wavelength range, which leads to a small increase in the noise in the white light curve. To make an order of magnitude estimate for the amplitude of this effect, we calculated the standard deviation of the secondary eclipse depths in all wavelength channels. It is ~100 ppm, which is comparable to the additional scatter we observed in the white light curve.

In addition to calculating the fit rms compared to the photon noise, we also tested for the presence of red noise based on whether the rms decreases as expected when the light curve in binned in time. If the noise is white (uncorrelated in time), the residuals are expected to decrease by a factor of $\sqrt{N}$, where N is the number of points in a bin. Figure 4 shows the binned residuals compared to expectations for white noise. The HST/WFC3 and Spitzer Channel 2 light curves agree well with expectations, whereas Spitzer Channel 1 shows higher noise than expected as bin size increases. This test confirms the presence of time-correlated noise in the Channel 1 light curve that can be seen by eye in the residuals in Figure 5. Both Spitzer channels use the same detector, but Channel 1 data are more susceptible to time-correlated noise because the point-

spread function is more undersampled at shorter wavelengths, making intrapixel sensitivity variations more pronounced.

## 4. Results

The fitted light curves are shown in Figure 5. This figure shows results from the spherical harmonics model for the thermal phase variation and has instrument systematics removed. Broadly speaking, the phase curves show large dayside planet-to-star flux values, ranging from 0.151% ± 0.015%, 0.446% ± 0.38%, and 0.569% ± 0.014% in the WFC3 white light curve, as well as the *Spitzer* 3.6 and 4.5 µm bandpasses, respectively. The planet flux changes significantly with orbital phase in all three of the data sets, suggesting a strong gradient from dayside to nightside temperature, and peak brightness occurs near phase 0.5. In this section, we quantitatively characterize the phase curve shape, split the data into phase resolved spectra, evaluate different temperature maps, compare with previous observations of the dayside thermal emission spectrum, and report the transmission spectrum.

### 4.1. Phase Curve Amplitudes and Hotspot Offsets

The shape of a phase curve can be summarized with two parameters: the amplitude of thermal phase variation (minimum to maximum brightness, divided by the secondary eclipse depth) and the location of peak brightness (typically called a “hotspot offset” and measured in degrees eastward of the substellar point). Table 2 lists the estimated amplitudes and hotspot offsets (median and $1\sigma$ credible interval) for the band-integrated WFC3 phase curve and both *Spitzer* channels. The estimates are from the sine curve model for the thermal phase variation. The advantage of using this model (even though it does not provide strictly the best fit) is that it directly fits the amplitude and offset as free parameters.

For all three phase curves, the hotspot offset is consistent with zero degrees, which could indicate a small ratio of radiative to advective timescales (the incident flux is reradiated to space faster than it is advected around to the nightside). Fast radiative timescales are predicted at high temperatures, and small hotspot offsets are also observed for other very hot Jupiters (Perez-Becker & Showman 2013; Komacek & Showman 2016; Komacek et al. 2017). The measured offsets are inconsistent with the trend reported in Zhang et al. (2018), which predicts the hotspot offset increases with planet temperature for irradiation temperatures greater than 3410 K. The Zhang et al. (2018) model predicts an eastward hotspot offset of 4.5 for WASP-103b, which is significantly larger than observed, hinting at diversity in the circulation patterns of the hottest planets. The phase curve amplitudes are large (near 0.8–0.9), as expected for an atmosphere with inefficient heat redistribution. In Section 6, we compare these results to expectations from 3D GCMs.

### 4.2. Phase-resolved Spectra

We used the best-fit phase curves (with systematics removed) to generate phase-resolved emission spectra. Because the spiderman thermal phase variation models fit the temperature of the planet rather than the eclipse depth directly, we estimated the dayside emission spectrum as follows. We used spiderman's *eclipse_depth* method to calculate the average planet-to-star flux for the best-fit model during secondary eclipse. To estimate uncertainties, we took the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** WASP-103b phase curve observations from HST/WFC3 (top) and Spitzer/IRAC (middle and bottom). For clarity, the data are phase-folded on the planet's orbital period and binned in 30 uniformly spaced bins between 0 and 1 (corresponding to 0.8 hr). The left column shows the phase curves with systematic noise removed (black points) compared to the best-fit spherical harmonics model (colored lines). The error bars denote $1\sigma$ uncertainties (in some cases, the errors are smaller than the data points). We include the transits in the fit, but they are not displayed in this figure. The right-hand column shows the binned residuals for the best-fit light curve. The gray error bars in the upper right of the left panels correspond to 500 ppm, to illustrate the changing y-axis scale.</figcaption></figure>

<table><caption>**Table 2**
Phase Curve Properties</caption><thead><tr><th>Bandpass</th><th>Source</th><th>Amplitude</th><th>Offset (Degrees)</th></tr></thead><tbody><tr><td rowspan="4">WFC3</td><td>data</td><td>0.91 ± 0.02</td><td>-0.3 ± 0.1</td></tr><tr><td>nominal GCM</td><td>0.89</td><td>15.32</td></tr><tr><td>[M/H] = 0.5 GCM</td><td>0.84</td><td>19.64</td></tr><tr><td>τ<sub>drag4</sub> GCM</td><td>0.97</td><td>2.34</td></tr><tr><td rowspan="4">Spitzer 3.6 µm</td><td>τ<sub>drag3</sub> GCM</td><td>0.99</td><td>0.18</td></tr><tr><td>data</td><td>0.86 ± 0.13</td><td>2.0 ± 0.7</td></tr><tr><td>nominal GCM</td><td>0.78</td><td>9.19</td></tr><tr><td>[M/H] = 0.5 GCM</td><td>0.72</td><td>12.79</td></tr><tr><td rowspan="4">Spitzer 4.5 µm</td><td>τ<sub>drag4</sub> GCM</td><td>0.86</td><td>0.90</td></tr><tr><td>τ<sub>drag3</sub> GCM</td><td>0.97</td><td>0.18</td></tr><tr><td>data</td><td>0.83 ± 0.05</td><td>1.0 ± 0.4</td></tr><tr><td>nominal GCM</td><td>0.79</td><td>8.11</td></tr><tr><td rowspan="4"></td><td>[M/H] = 0.5 GCM</td><td>0.73</td><td>11.35</td></tr><tr><td>τ<sub>drag4</sub> GCM</td><td>0.85</td><td>0.90</td></tr><tr><td>τ<sub>drag3</sub> GCM</td><td>0.93</td><td>0.18</td></tr></tbody></table>

standard deviation of the residuals of the in-eclipse data points, then added this value in quadrature to the standard deviation of the residuals of the out-of-eclipse data. This quadrature sum accounts for the uncertainty in the baseline flux. To account for red noise in the *Spitzer* 3.6 µm light curve, we use the approach of Pont et al. (2006) to determine the red noise contribution on the timescale of the eclipse. We add the estimated red noise in quadrature, which increases the uncertainty on planet-to-star flux by a factor of 2.5.

For the other orbital phases, we binned the light curve (with systematics removed) in eight intervals of about 0.1 in orbital

phase (2.2 hr), with endpoints at phases 0.06, 0.15, 0.25, 0.35, 0.44 and 0.56, 0.65, 0.75, 0.85, and 0.94. These endpoints were chosen to ensure that there is no contribution from in-transit or in-eclipse data. In each phase bin, we estimated the planet-to-star flux from the mean value of the data points in the bin. To estimate the uncertainty, we took the standard deviation of the points in the bin and added it in quadrature to the standard deviation of the data points during secondary eclipse (phase 0.46–0.54), to account for the uncertainty in baseline stellar flux. For the 3.6 µm data, we also add red noise on the timescale of a phase bin, following Pont et al. (2006). The phase-resolved emission spectra are shown in Figure 6 and listed in Table 3. We show the dayside spectrum in Figure 7.

To test that the phase-resolved planet-to-star flux values are robust to different approaches for fitting the phase curves, we compared the estimated planet-to-star flux for all four of the thermal phase variation models (sinusoid, kinematic, spherical harmonics, and two-temperature). Because the systematic noise is not strongly correlated with the astrophysical signal, the systematics-divided data are nearly identical for all the models. This point is illustrated in Figure 8 for the broadband WFC3 light curve. We found that the choice of model generally does not significantly change the estimated planet-to-star flux ratios. The estimates agree to better than one sigma for 90% of phase bins for the spherical harmonics, two-temperature, and physical models. For the WFC3 data, the sinusoid is higher than the other models by an average of $1.5\sigma$ for phases 0.5–1. This discrepancy may be due to the added flexibility in hotspot offset for the sinusoid model; other models do not allow for westward hotspot offsets.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Phase-resolved planet-to-star flux ratios (points) compared to the best-fit blackbody (blue line) and the GCM with $\tau_{\text{drag4}}$. Phases $\phi = 0.5$ and $0.0$ correspond to times of secondary eclipse and transit, when the substellar and antistellar points are respectively aimed at Earth. The blackbody model is fit to the HST/WFC3 data only. We also show a 100 ppm error bar in the lower right of each panel to emphasize the changing limits on the y-axes.</figcaption></figure>

### 4.3. Comparison with Previous Eclipse Observations

We compared our results to the dayside emission spectrum reported by Cartier et al. (2017), which is based on a Gaussian process analysis of two secondary eclipses from HST/WFC3. The shape of their spectrum is consistent with what we find, but their eclipse depths are 125 ppm smaller on average (a difference of about 10%). A likely explanation for this difference is that the Cartier et al. (2017) analysis does not include the planet's thermal phase variation, so that the Gaussian process models it as an instrument systematic. If the phase variation is absorbed into the systematic model, the measured eclipse depths would be biased low. By visual inspection of Figure 8, we estimate the amplitude of this effect is ~100 ppm, which is comparable to the offset between the two analyses. Our estimated uncertainties are a factor of four smaller than those reported in Cartier et al. (2017), which is consistent with photon-limited expectations (our data set includes two additional eclipses and an out-of-eclipse baseline that is longer by a factor of five, as well as 60% larger wavelength bins).

We also compared our dayside spectrum to the $z'$- and $K_S$-band secondary eclipse depths reported in Delrez et al. (2018). The $z'$ (0.9 µm) eclipse is 1.0σ lower than our best-fit blackbody spectrum (described in Section 5), and the $K_S$ (2.1 µm) measurement is higher than the model by 2.5σ. Because these results are consistent with (but less precise than) the WFC3 data, we do not include them in our analysis of the atmospheric composition, but we encourage additional measurements in the $K_S$ band in order to confidently determine whether an emission feature is present at those wavelengths.

### 4.4. Transmission Spectrum

Each phase curve observation includes a transit of WASP-103b. To measure the wavelength-dependent transit depths (the transmission spectrum), we select a subset of each phase curve. The subset includes the transit and additional baseline on either side, such that the total light curve has twice the duration of the transit. Over this short duration, there is negligible curvature in the light curve due to the planet's thermal phase variation. We fit the data with a transit model, which has free parameters for the planet-to-star radius and a linear limb-darkening parameter.

The orbital parameters (inclination, $a/R_s$, period, and time of central transit) were fixed at previously published values listed in Section 3. We fit for the instrument systematics using the same model as for the phase curve fits, except that we modeled the visit-long systematic as a linear trend in time (which is sufficient for the shorter duration). The advantage of fitting the transits separately from the full phase curves is that the resultant transit depths are not dependent on how the phase variation is modeled. The transit light curve fits are shown in Figure 9.

We show the measured transmission spectrum in Figure 10. The spectrum is biased by flux contamination from the planet's nightside; to correct for nightside emission, we estimate the average nightside planet-to-star flux ratio and subtract it from the transit depth. We calculate planet-to-star flux using a NextGen spectrum for the star (interpolated to $T_{\text{eff}} = 6110$ K), a blackbody for the planet, and a planet-to-star radius ratio of 0.1146. We assume a nightside temperature of 1700 K, which is near the median of the nightside temperatures estimated from the phase variation models (Table 4). We also show the uncorrected transit depths. The corrected and uncorrected transit depths are listed in Table 5. For the uncorrected data, there is an offset between the HST and Spitzer data of more than five atmospheric scale heights. The uncorrected spectrum is inconsistent with a flat line at 5.3σ confidence, whereas the nightside-corrected spectrum is consistent within 1σ. The corrected spectrum is also consistent with predictions from the $\tau_{\text{drag4}}$ GCM, which shows water features in the WFC3 bandpass. Future observations with higher precision could reveal these features.

## 5. Atmospheric Composition and Thermal Structure

We characterized the planet's atmospheric composition by fitting 1D models to the phase-resolved emission spectra. First, we modeled the planet flux as a simple blackbody to estimate the dayside brightness temperature and test for significant absorption or emission features. We then performed a more sophisticated grid-based retrieval to estimate the atmospheric metallicity, carbon-to-oxygen ratio, and thermal structure. We also evaluated the climate based on the best-fit spiderman temperature maps.
<table><caption>**Table 3**
Phase-resolved Emission Spectra</caption><thead><tr><th>λ</th><th>Dilution</th><th>φ = 0.1</th><th>φ = 0.2</th><th>φ = 0.3</th><th>φ = 0.4</th><th>φ = 0.5</th><th>φ = 0.6</th><th>φ = 0.7</th><th>φ = 0.8</th><th>φ = 0.9</th></tr></thead><tbody><tr><td>1.175</td><td>0.10</td><td>179 ± 79</td><td>411 ± 77</td><td>647 ± 80</td><td>1143 ± 65</td><td>1259 ± 47</td><td>1063 ± 64</td><td>710 ± 73</td><td>412 ± 78</td><td>177 ± 79</td></tr><tr><td>1.225</td><td>0.11</td><td>188 ± 76</td><td>398 ± 74</td><td>928 ± 77</td><td>1276 ± 62</td><td>1448 ± 46</td><td>1216 ± 62</td><td>888 ± 71</td><td>539 ± 75</td><td>280 ± 76</td></tr><tr><td>1.275</td><td>0.11</td><td>166 ± 76</td><td>379 ± 74</td><td>869 ± 77</td><td>1323 ± 62</td><td>1480 ± 46</td><td>1282 ± 62</td><td>814 ± 71</td><td>515 ± 75</td><td>247 ± 76</td></tr><tr><td>1.325</td><td>0.11</td><td>266 ± 75</td><td>432 ± 73</td><td>904 ± 76</td><td>1357 ± 62</td><td>1498 ± 45</td><td>1267 ± 61</td><td>925 ± 70</td><td>552 ± 74</td><td>333 ± 75</td></tr><tr><td>1.375</td><td>0.12</td><td>189 ± 81</td><td>514 ± 78</td><td>928 ± 82</td><td>1376 ± 66</td><td>1611 ± 48</td><td>1411 ± 65</td><td>954 ± 75</td><td>461 ± 79</td><td>292 ± 81</td></tr><tr><td>1.425</td><td>0.13</td><td>238 ± 79</td><td>532 ± 76</td><td>1198 ± 79</td><td>1431 ± 64</td><td>1718 ± 47</td><td>1511 ± 64</td><td>1063 ± 73</td><td>605 ± 77</td><td>338 ± 79</td></tr><tr><td>1.475</td><td>0.14</td><td>191 ± 81</td><td>527 ± 79</td><td>1068 ± 82</td><td>1460 ± 66</td><td>1667 ± 48</td><td>1392 ± 66</td><td>1090 ± 75</td><td>580 ± 80</td><td>268 ± 81</td></tr><tr><td>1.525</td><td>0.14</td><td>143 ± 84</td><td>478 ± 81</td><td>1048 ± 85</td><td>1429 ± 69</td><td>1623 ± 50</td><td>1367 ± 68</td><td>943 ± 77</td><td>607 ± 82</td><td>291 ± 84</td></tr><tr><td>1.575</td><td>0.15</td><td>367 ± 88</td><td>761 ± 85</td><td>1088 ± 89</td><td>1581 ± 72</td><td>1749 ± 52</td><td>1503 ± 71</td><td>1107 ± 81</td><td>754 ± 86</td><td>422 ± 88</td></tr><tr><td>1.625</td><td>0.16</td><td>359 ± 93</td><td>565 ± 90</td><td>1169 ± 94</td><td>1590 ± 76</td><td>1843 ± 56</td><td>1593 ± 75</td><td>1142 ± 86</td><td>542 ± 91</td><td>351 ± 93</td></tr><tr><td>3.6</td><td>0.17</td><td>982 ± 271</td><td>3474 ± 268</td><td>4309 ± 255</td><td>4060 ± 248</td><td>4458 ± 383</td><td>3524 ± 249</td><td>3725 ± 267</td><td>1865 ± 269</td><td>1116 ± 272</td></tr><tr><td>4.5</td><td>0.16</td><td>1560 ± 220</td><td>3347 ± 213</td><td>4150 ± 189</td><td>5240 ± 178</td><td>5686 ± 138</td><td>4995 ± 181</td><td>3677 ± 212</td><td>2403 ± 213</td><td>921 ± 219</td></tr></tbody></table>

Note. The planet-to-star flux in each phase bin $\phi$ is in units of ppm.

## 5.1. Blackbody Fits

A blackbody is the simplest model for the planet's thermal emission and provides a useful first evaluation of the atmospheric properties. To calculate the best-fit blackbody model, we assumed a planet-to-star radius ratio of 0.1146 and used a NextGen stellar spectrum interpolated to an effective temperature of 6110 K (Allard et al. 2012). We calculated the best fit with a least-squares fitting routine. To determine uncertainties on the planet brightness temperature, we performed an MCMC analysis with free parameters for the planet temperature and the stellar $T_{\text{eff}}$. We used a Gaussian prior on $T_{\text{eff}}$ of 6110 ± 160 K (Gillon et al. 2014).

The best-fit blackbodies are shown in Figures 6 and 7. The model fits the data fairly well: it is consistent with the data at the $3\sigma$ level for 70% of the orbital phases. However, the *Spitzer* data have higher brightness temperatures than the WFC3 data at phase 0.5, and lower brightness temperatures at phases 0.8–0.9. These differences suggest the presence of an emission feature on the dayside at *Spitzer* wavelengths, which transitions to an absorption feature on the nightside, perhaps indicating changes in thermal structure with longitude in the atmosphere.

We also fit independent blackbody models to the *HST*/WFC3 data and each *Spitzer* channel separately. The resulting brightness temperatures and $1\sigma$ uncertainties are listed in Table 6. The WFC3 data agree well with a blackbody model at all orbital phases except phase 0.5 (although still consistent at better than $1.5\sigma$). The dayside has higher signal-to-noise than the other orbital phases, thanks to the two secondary eclipse observations from Cartier et al. (2017). The more sophisticated grid-based retrieval (described in the next section) provides a better fit to the dayside.

We note that the uncertainties on the brightness temperatures in different bandpasses are correlated with each other because they include the uncertainty on the stellar temperature: i.e., if the stellar temperature increases, so do the brightness temperatures. To evaluate the significance of features in the emission spectra, we hold the stellar spectrum fixed in the retrieval analysis.

## 5.2. Grid-based Retrieval

To infer abundances from the dayside spectrum (phase 0.46–0.54), we use a self-consistent grid-based method (ScCHIMERA) similar to that employed in Arcangeli et al. (2018) and Mansfield et al. (2018). We generated a grid from

one-dimensional forward models of the atmosphere over a broad range of metallicities (M/H), carbon-to-oxygen ratios (C/O), and stellar irradiation ($f$). The $f$ parameter is a scaling factor for the stellar flux at the top of the atmosphere, where $f=1$ corresponds to full heat redistribution and $f=2$ is equivalent to only allowing the dayside to re-radiate.

At each point in the grid, we compute forward models to determine self-consistent, radiative-convective temperature-pressure (T-P) profiles. We determine the molecular abundances in each atmospheric layer, assuming thermochemical equilibrium (calculated with the NASA CEA routine; Gordon & McBride 1994). We include opacity from the major absorbers expected for a hot Jupiter atmosphere, including H$_2$O, CO, CO$_2$, TiO, VO, FeH, and H$_2$-H$_2$ CIA. Notably, in contrast to most prior atmospheric retrievals for the hottest planets, we also included opacity from H$^-$, which is an important absorber at temperatures above 2500 K (Arcangeli et al. 2018; Parmentier et al. 2018). Using these opacities and T-P profiles, we calculated thermal emission spectra over the full grid, using the CHIMERA retrieval suite (described in Line et al. 2013, 2014). We then explored the grid with an MCMC chain, using the emcee package (Foreman-Mackey et al. 2013), and interpolated in the grid to calculate the likelihood at each model step. The priors were uniform over the ranges $0.2 \le f \le 2.8$, $-1 \le \log M/H \le 2.5$, and $-1 \le \log C/O \le 0.95$. We fit this model to the dayside and nightside emission spectra (phases $\phi=0.5$ and $0.1$).

### 5.2.1. Dayside Spectrum

The main characteristics of the dayside emission spectrum are: (1) it is blackbody-like at WFC3 wavelengths, and (2) in the *Spitzer* bands, the planet-to-star flux is significantly higher than predicted for the best-fit blackbody, indicating an emission feature. The best-fit spectrum reproduces these data fairly well, with $\chi^2 = 1.77$ (for nine degrees of freedom). The largest contribution to the $\chi^2$ value is the 4.5 µm eclipse depth, which is larger than the best-fit model prediction by $2.9\sigma$. When the 4.5 µm point is removed, the fit has $\chi^2 = 1.17$ (eight degrees of freedom). The best-fit model has a moderately enhanced metallicity (23× solar), carbon-to-oxygen equal to 0.76, poor heat redistribution, and a thermal inversion (temperature increasing with altitude).

Figure 11 shows the opacity contributions of key absorbers for the best-fit model. In the optical (which we do not observe
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Dayside planet-to-star flux (points) compared to a blackbody model for the planet (dashed gray line) and the best-fit 1D model (blue line, with $1\sigma$ uncertainty shaded in orange) for the HST/WFC3 data (left) and Spitzer (right). The 1D model was fit to all the data simultaneously, but the blackbody was fit to the WFC3 data only. The best-fit 1D model binned at the resolution of the Spitzer data is indicated by blue squares (right).</figcaption></figure>

directly), there is strong absorption by TiO, VO, and FeH. In the near-infrared, H$_2$O, H$^-$, and hydrides/oxides all contribute to the opacity, leading to nearly constant opacity over the WFC3 wavelength range. In cooler atmospheres, water is the dominant absorber over this bandpass (e.g., Kreidberg et al. 2014a; Line et al. 2016), but in WASP-103b, H$_2$O is partially dissociated in the photosphere, leading to a drop in abundance by a factor of ~10 (see Figure 12). Water also has intrinsically weaker features at high temperature (e.g., Tinetti et al. 2012). On top of this, there is significant H$^-$ opacity from single H atoms bound with free electrons, which fills in the opacity at wavelengths shorter than 1.5 µm. Finally, the sharp vertical gradient in water abundance results in water becoming optically thick over a very narrow range in pressure, where temperature is nearly constant. Taken together, all these factors add up to produce a nearly featureless spectrum from 1.1–1.7 µm. Finally, in the infrared, the dominant absorber is CO, which produces the emission feature at Spitzer wavelengths.

Figure 12 shows a summary of the temperature–pressure profile and abundances for the best-fit model. The T–P profile is inverted, with temperature increasing from 2800 to 3500 K over the pressure range $10^{-2}$–$10^{-3}$ bar. The thermal inversion is probably driven by absorption of optical light by oxides and hydrides in the upper atmosphere and the absence of cooling by water molecules (which have dissociated). The observations are sensitive to pressures of ~0.01–0.001 bar, which spans the tropopause, where temperature begins to increase and the water abundance drops by more than an order of magnitude.

In Figure 13, we show the posterior distributions from the grid retrieval. We infer a range in metallicity of $23_{-13}^{+29}$ × solar, somewhat higher than expected based on Jupiter's metal enrichment ($3-5$ × solar; Wong et al. 2004) and the trend toward decreasing metallicity with increasing planet mass observed for the solar system and exoplanets (e.g., Kreidberg et al. 2014a). The metallicity is super-solar at $>3\sigma$ confidence. However, planet population synthesis models predict some scatter in atmospheric metallicity. Planets near WASP-103b's mass ($1.5 M_{Jup}$) are expected to have metallicities in the range of roughly $1-10 \times$ solar (Fortney et al. 2013; Mordasini et al. 2016). Our result for WASP-103b lies on the upper end of this range, and may be indicative of intrinsic scatter in the mass-atmospheric metallicity relation.

The retrieved C/O is consistent with solar, with a $1\sigma$ confidence interval of 0.54–0.85. We infer an upper limit on C/O of 0.9 at $3\sigma$ confidence, driven by the fact that the atmospheric chemistry is expected to change dramatically when C/O exceeds unity. For a carbon-rich composition, the equilibrium abundance of methane relative to CO increases by orders of magnitude compared to an oxygen-rich composition (e.g., Madhusudhan et al. 2011). Our Spitzer eclipse depths are sensitive to the relative abundance of these species, so we can confidently rule out a carbon-rich composition despite the lack of spectrally resolved features (assuming the atmosphere is in chemical equilibrium).

We infer a heat redistribution $f = 2.49_{-0.15}^{+0.14}$. An $f$ parameter of unity represents isotropic heat distribution, whereas $f = 2$ corresponds to dayside emission only. We estimate $f > 2$, indicating a thermal inversion and likely inefficient transport of heat to the nightside. The heat redistribution is strongly correlated with atmospheric metallicity because increasing metallicity shifts the T–P profile to lower pressures, resulting in hotter temperatures at a given pressure level (equivalent to less efficient heat redistribution).

Our analysis comes with several important caveats.

1. The best-fit model is not a perfect fit to the data (with $\chi^2 = 1.77$ for nine degrees of freedom), so the uncertainties produced by the MCMC may be underestimated.

2. The inferred C/O and metallicity are highly sensitive to the planet-to-star flux at Spitzer $4.5 \mu m$, which is the worst-fit data point. To fit this data point, the model favors super-solar metallicities and C/O, which drive up the CO abundance (the dominant absorber at $4.5 \mu m$).

3. The Spitzer $4.5 \mu m$ data is from broadband photometry, so the inferred CO feature is not spectrally resolved. It is possible that unknown absorbers or disequilibrium chemistry affect the $4.5 \mu m$ planet-to-star flux, but we cannot uniquely identify these features in our spectrum.

We therefore caution against over-interpreting these results until wider spectral coverage is available.

### 5.2.2. Nightside Spectrum

We also fit the nightside spectrum (phase 0.1) with the grid-based retrieval. The best-fit spectrum has a non-inverted
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left: fits to the broadband WFC3 phase curve compared to a GCM. The colored lines correspond to different temperature maps fit to the data, and the dashed gray line is from the $\tau_{drag4}$ GCM. We also show the measured planet-to-star flux for each map (points), which is model-dependent due to slight degeneracies with the instrument systematic model. Right: temperature maps from the best-fit models and the GCM at a pressure of 0.1 bar.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Best-fit transit light curves for WASP-103b (lines), compared to binned data (black points). From top-to-bottom, we show the band-integrated WFC3 light curve, *Spitzer* 3.6, and 4.5 µm.</figcaption></figure>

temperature pressure profile. At $1\sigma$ confidence, the metallicity is 15–240× solar and the C/O is unbounded over the full prior range. The atmospheric composition is consistent with results from the dayside spectrum.

This agreement is an encouraging sanity check; however, there are several model assumptions that may result in artificially tight constraints on the atmospheric properties on the nightside. One challenge in modeling the nightside spectrum is that the physical processes shaping the T–P profile are unknown. Our model assumes a scaled stellar irradiation at the top of the atmosphere, but in reality, the heat source is advection from the dayside. Another caveat is that the model is not self-consistent: the energy leaving the dayside is not constrained to equal the energy entering the nightside. Further work is needed to develop a fully self-consistent 2D retrieval method for phase curve observations.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** The transmission spectrum of WASP-103b, corrected for planet nightside emission at 1700 K (blue points) and uncorrected (red points). The dark gray line is the model transmission spectrum from the $\tau_{drag4}$ GCM, with squares indicating the model binned over the *Spitzer* bandpasses. Atmospheric scale height $H$ is shown on the right y-axis, where $H = 5.5 \times 10^6$ m (assuming a mean molecular weight of 2.3 atomic mass units, surface gravity $g = 15.9 \text{ m s}^{-2}$, and a temperature $T = 2410$ K).</figcaption></figure>

<table><caption>**Table 4**
Model Comparison</caption><thead><tr><th>Data</th><th>Model</th><th>T<sub>min</sub></th><th>T<sub>max</sub></th><th>T&#773;<sub>night</sub></th><th>T&#773;<sub>day</sub></th><th>&Delta;<sub>BIC</sub></th></tr></thead><tbody><tr><td rowspan="4">WFC3</td><td>Sph. Harmonics</td><td>1227</td><td>3237</td><td>1822</td><td>2636</td><td>0</td></tr><tr><td>Kinematic</td><td>1977</td><td>3953</td><td>1977</td><td>2769</td><td>14</td></tr><tr><td>Two Temp.</td><td>0</td><td>2879</td><td>0</td><td>2879</td><td>42</td></tr><tr><td>Sinusoid</td><td>...</td><td>...</td><td>...</td><td>...</td><td>17</td></tr><tr><td rowspan="4">Ch 1</td><td>Sph. Harmonics</td><td>1269</td><td>3391</td><td>1912</td><td>2741</td><td>0</td></tr><tr><td>Kinematic</td><td>1932</td><td>3630</td><td>1975</td><td>2614</td><td>34</td></tr><tr><td>Two Temp.</td><td>1418</td><td>2990</td><td>1418</td><td>2990</td><td>11</td></tr><tr><td>Sinusoid</td><td>...</td><td>...</td><td>...</td><td>...</td><td>25</td></tr><tr><td rowspan="4">Ch 2</td><td>Sph. Harmonics</td><td>888</td><td>3714</td><td>1729</td><td>2864</td><td>2</td></tr><tr><td>Kinematic</td><td>1614</td><td>3931</td><td>1621</td><td>2544</td><td>15</td></tr><tr><td>Two Temp.</td><td>1344</td><td>3241</td><td>1344</td><td>3241</td><td>0</td></tr><tr><td>Sinusoid</td><td>...</td><td>...</td><td>...</td><td>...</td><td>22</td></tr></tbody></table>

As a test, we also calculated the difference in brightness temperature between the *HST* and *Spitzer* 4.5 µm data (reported in Table 6) for both nightside phases ($\phi = 0.1$ and 0.9). At phase 0.1 and 0.9, the *Spitzer* temperature is lower by 2.7 and $5.6\sigma$, respectively. These values are a lower limit to the
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Abundance-weighted absorption cross-sections illustrating the important opacity sources at the photospheric pressure and temperature (5 mbar, 3036 K). The strong CO feature at 4.5 µm contributes to the high planet-to-star flux at that wavelength. Water, hydrides/oxides, and the H⁻ bound-free opacities all play a role in shaping the HST/WFC3 spectrum.</figcaption></figure>

<table><caption>**Table 5**
WASP-103b Transmission Spectrum</caption><thead><tr><th>Wavelength<br>(micron)</th><th>(R<sub>p</sub>/R<sub>s</sub>)<sup>2</sup> (%)<br>(T<sub>n</sub> = 0 K)</th><th>(R<sub>p</sub>/R<sub>s</sub>)<sup>2</sup> (%)<br>(T<sub>n</sub> = 1700 K)</th><th>Error<br>(%)</th></tr></thead><tbody><tr><td>1.175</td><td>1.3178</td><td>1.3115</td><td>0.0092</td></tr><tr><td>1.225</td><td>1.3144</td><td>1.3067</td><td>0.0087</td></tr><tr><td>1.275</td><td>1.3195</td><td>1.3103</td><td>0.0086</td></tr><tr><td>1.325</td><td>1.3454</td><td>1.3345</td><td>0.0087</td></tr><tr><td>1.375</td><td>1.3512</td><td>1.3385</td><td>0.0089</td></tr><tr><td>1.425</td><td>1.3588</td><td>1.3441</td><td>0.0091</td></tr><tr><td>1.475</td><td>1.3364</td><td>1.3197</td><td>0.0092</td></tr><tr><td>1.525</td><td>1.3522</td><td>1.3333</td><td>0.0093</td></tr><tr><td>1.575</td><td>1.3688</td><td>1.3476</td><td>0.0100</td></tr><tr><td>1.625</td><td>1.3549</td><td>1.3314</td><td>0.0104</td></tr><tr><td>3.6</td><td>1.4013</td><td>1.3238</td><td>0.0328</td></tr><tr><td>4.5</td><td>1.4329</td><td>1.3554</td><td>0.0264</td></tr></tbody></table>

**Notes.** Wavelength-dependent transit depths, corrected for companion dilution and nightside flux (assuming nightside temperatures of 0 and 1700 K for the second and third columns, respectively). The error corresponds to the 68% credible interval from an MCMC fit to the transit light curves.

significance, because the brightness temperatures noted in the table also include the uncertainty in the stellar $T_{eff}$ (which increases the uncertainty on the absolute planet temperature, but not the relative temperatures that are relevant for this calculation). The drop in brightness temperature is more significant at phase 0.9 than it is at 0.1, providing further evidence, in addition to the phase 0.1 retrieval, that the nightside temperature pressure profile is not inverted.

### 5.3. Climate

We fit three different models to characterize the planet's climate: a two-temperature map, the physically motivated kinematic model of Zhang & Showman (2017), and a spherical harmonics map. We also fit the thermal phase variation with a sinusoid, which can be inverted to map the climate (Cowan & Agol 2008; Cowan & Fujii 2017). All of the models provide

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Summary of the 1D self-consistent model atmosphere fits to the dayside emission spectrum. The temperature-pressure profile (top axis) is indicated by the $1\sigma$ spread of 500 randomly drawn T-P profiles from the posterior (light red) and a representative fit, with $f = 0.4$, [M/H] = 1.5, and C/O = 0.7 (dark red). The normalized thermal emission contribution functions for the *Spitzer* points are shown in solid red, the WFC3 in-water band (1.33–1.48 µm) in dark blue, and WFC3 out-of-water-band in light blue. The observations probe between ~0.01 and 0.001 bar, just above the tropopause region of the atmosphere where the temperature is increasing. The dashed curves are thermochemical equilibrium mixing ratios for important absorbers, computed along the representative fit's self-consistent T-P profile. Note the rapid dissociation of water above the ~10 mbar level where the inversion begins.</figcaption></figure>

<table><caption>**Table 6**
Phase-resolved Brightness Temperatures</caption><thead><tr><th>Orbital Phase</th><th>T<sub>b</sub><br>WFC3</th><th>T<sub>b</sub><br>Ch. 1</th><th>T<sub>b</sub><br>Ch. 2</th><th>χ<sup>2</sup><sub>ν</sub><br>(9 DoF)</th></tr></thead><tbody><tr><td>0.06–0.15</td><td>1883 ± 41</td><td>1523 ± 153</td><td>1589 ± 105</td><td>0.7</td></tr><tr><td>0.15–0.25</td><td>2208 ± 33</td><td>2612 ± 117</td><td>2299 ± 100</td><td>0.9</td></tr><tr><td>0.25–0.35</td><td>2587 ± 37</td><td>2926 ± 114</td><td>2592 ± 96</td><td>1.5</td></tr><tr><td>0.35–0.44</td><td>2831 ± 39</td><td>2834 ± 111</td><td>2976 ± 93</td><td>1.5</td></tr><tr><td>0.44–0.56</td><td>2933 ± 41</td><td>2995 ± 159</td><td>3154 ± 99</td><td>2.8</td></tr><tr><td>0.56–0.65</td><td>2811 ± 39</td><td>2631 ± 110</td><td>2891 ± 97</td><td>2.0</td></tr><tr><td>0.65–0.75</td><td>2572 ± 36</td><td>2708 ± 117</td><td>2421 ± 100</td><td>1.3</td></tr><tr><td>0.75–0.85</td><td>2263 ± 33</td><td>1952 ± 125</td><td>1939 ± 99</td><td>1.3</td></tr><tr><td>0.85–0.94</td><td>1987 ± 37</td><td>1594 ± 145</td><td>1288 ± 118</td><td>0.6</td></tr></tbody></table>

**Note.** $\chi^2$ values are for the fits to the WFC3 data only.

reasonable fits to the data, with $\chi^2$ near unity, but they yield significantly different temperature maps. Table 4 lists the best-fit minimum and maximum temperatures, as well as the mean day- and nightside temperatures. We also list the information criterion (BIC) values for the fits (a $\Delta$BIC value greater than 10 constitutes strong evidence against a given model; Kass & Raftery 1995).

The spherical harmonics map generally fits the data the best. It has a lower BIC value than all the other models for the broadband WFC3 and *Spitzer* 3.6 µm phase curves. For the 4.5 µm phase curve, the two-temperature map provides the best fit, but it only lowers the BIC value by two relative to the spherical harmonics map, which is not a statistically significant improvement (Kass & Raftery 1995). The spherical harmonics model yields a mean dayside temperature near 2700 K, whereas
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Posterior distributions for WASP-103b's atmospheric heat redistribution, metallicity, and C/O, from a grid-based fit to the dayside emission spectrum. The histograms on the diagonal show the marginalized distribution of each parameter, with dashed lines indicating the median and surrounding 68% credible interval. The blue lines correspond to solar metallicity (1) and C/O (0.54). The 2D histograms mark the $1\sigma$, $2\sigma$, and $3\sigma$ credible regions in dark, medium, and light blue, respectively.</figcaption></figure>

the nightside is closer to 1800 K, in good agreement with the blackbody fits to the phase-resolved spectra (see Section 5.1). The other models produce more extreme day-night temperature gradients. Between the models, there are substantial differences in the minimum and maximum temperatures (sometimes over 1000 K), whereas the day and nightside means are in better agreement (generally matching to within 250 K). This behavior is not surprising: a wide range of temperature gradients can yield similar average temperatures when integrated over the disk of the planet.

## 6. Comparison with GCMs

To explore the three-dimensional effects of atmospheric dynamics, we ran several GCMs to compare with the measured phase-resolved spectra. We simulated the atmospheric circulation and thermal structure of the planet using the combined SPARC/MITgcm model (Showman et al. 2009). The model solves the primitive equations in spherical geometry using the MITgcm (Adcroft et al. 2004) and the radiative transfer equations using a state-of-the-art one-dimensional radiative transfer model (Marley & McKay 1999). The code represents the opacities as correlated-k tables based on the line-by-line opacities described in Visscher et al. (2006) and Freedman et al. (2014). Our fiducial model assumes a solar composition with elemental abundances of Lodders & Fegley (2002) and the chemical equilibrium gas phase composition from Visscher et al. (2006). These calculations take into account the presence of $H^-$ opacities and the effect of molecular dissociation on the abundances. We used a timestep of 25s, ran the simulations for 300 days, and averaged all quantities over the last 100 days.

Overall, our setup is the same as described in Parmentier et al. (2016), but uses parameters specific to the WASP-103 system.

Several physical processes can reduce the ability of the atmosphere to transport heat efficiently through advection and change the overall circulation pattern. Among them, ohmic drag is thought to be an important phenomenon in the ionized environment of extremely hot Jupiters (Perna et al. 2010). We parameterize this effect as a Rayleigh drag with a drag constant $\tau_{drag}$ constant with pressure (Showman et al. 2013). Varying $\tau_{drag}$ from large values (i.e., weak drag) to small values (i.e., strong drag), the atmospheric circulation is expected to shift from a jet-dominated regime to a more axisymmetric circulation pattern going from the substellar to anti-stellar point. Moderate drag timescales are expected to change the circulation pattern—and thus significantly reduce the shift of the hottest point of the atmosphere—whereas short drag timescales are also expected to change the strength of the winds and thus the atmospheric day/night contrast (Komacek & Showman 2016; Komacek et al. 2017). Although Rayleigh drag is an incomplete representation of the complex magneto-hydrodynamic effects expected in these atmospheres (Batygin et al. 2013; Rogers & Komacek 2014; Rogers & Showman 2014; Rogers 2017), it nonetheless provides an estimate of the strength of the drag mechanism necessary to match the observations (Komacek et al. 2017; Parmentier & Crossfield 2017; Koll & Komacek 2018).

Our nominal GCM was a cloud-free, solar composition atmosphere with TiO/VO opacity and no added drag. Each GCM run is computationally intensive, so we ran a small number of additional models to see which parameters had the largest effect on the planet spectrum. We changed model parameters one at a time, considering cases with enhanced metallicity ([M/H] = 0.5), no TiO/VO, and added atmospheric drag with timescales $\tau_{drag} = 10^3$ and $10^4$ s, which we label $\tau_{drag3}$ and $\tau_{drag4}$, respectively. The GCM results are shown in Figure 14. To assess how well the GCM predictions reproduce the data, we calculated the amplitude and hot spot offset for all the models (listed in Table 2). The small observed hotspot offsets (–0.3–2.0 degrees) are best reproduced by the $\tau_{drag4}$ model, which has a smaller offset than the drag-free GCMs due to changes in wind pattern. In the drag models, the winds shift from a substellar to an antistellar flow rather than an equatorial jet. The $\tau_{drag4}$ model also provides a match to the observed phase curve amplitudes.

We also compared the TP profiles from the $\tau_{drag4}$ GCM to cloud condensation curves and the best-fit radiative-convective equilibrium models from the 1D retrieval (Figure 15). For the dayside photosphere, the TP profile slope and absolute temperature are in rough agreement between the 1D best fit and the GCM. At higher pressures, the GCM is systematically cooler, which is likely due to the effect of atmospheric circulation (at these pressures, the GCM mixes the temperature planet-wide). At lower pressures, the GCM is also cooler than the 1D fits, which may be due to metallicity differences between the models. The GCM has solar metallicity, whereas the best-fit 1D model has [M/H] ~ 1. Higher metallicity compositions have larger TiO/H$_2$O ratios, and because the pressure dependence of TiO dissociation is not as strong as for water dissociation, we expect stronger inversions for higher-metallicity atmospheres (Parmentier et al. 2018). On the nightside, we also find that the GCM is cooler than the 1D models. While in the 1D model, the day-to-night redistribution
<figure><img src="image_17.png" /><figcaption>**Figure 14.** GCM predictions (colored lines) compared to the best-fit spherical harmonics model for the WFC3 white light, Spitzer 3.6, and Spitzer 4.5 µm phase curves (black lines, top to bottom). The nominal model is solar composition and cloud-free, with TiO/VO opacity and no drag. The models are corrected for the predicted ellipsoidal variability of the planet.</figcaption></figure>

is fitted to the data, it is not a tunable parameter in the 3D GCM. There are several physical processes not included in the GCM that could contribute to a hotter nightside, including shocks, longitude-dependent drag, and latent heat released from H2 recombination (Bell & Cowan 2018). The best-fit nightside TP profile is hotter than the condensation curves through most of the photosphere, suggesting that the observable atmosphere is relatively free of clouds. This prediction could be tested with longer-wavelength phase curve observations.

We also compared the GCM output to temperature maps retrieved with *spiderman*. Figure 8 shows the 0.1 bar temperature map for the $\tau_{\text{drag4}}$ GCM compared to the best fit models. At this pressure, the GCM has minimum and maximum temperatures of 920 and 3360 K. The temperature gradient from the dayside to the terminator is intermediate between the kinematic and spherical harmonics models. The GCM predicts a cooler nightside than all models except the two-temperature model. We note that none of the models are perfect: there is degeneracy in the *spiderman* maps, with large differences in climate producing reasonably good fits to the phase curves (see Table 4), whereas the GCM is too cold on

<figure><img src="image_18.png" /><figcaption>**Figure 15.** TP profiles from the GCM and 1D retrieval, compared to condensation curves of potential cloud species. The thin blue/cyan lines correspond to a sample of randomly drawn nightside TP profiles from the GCM, and the thin red/orange lines correspond to dayside TP profiles. The darker colors (blue/red) indicate the extent of the contribution function (encompassing 80% of the emitted flux). The thick red and blue lines are the best-fit 1D TP profiles for the phase 0.5 (dayside average) and 0.1 (nightside average). The dashed lines are condensation curves for a range of possible cloud species. The dotted lines correspond to regions of constant H$_2$O abundance, with numbers indicating the log10 (H$_2$O volume mixing ratio).</figcaption></figure>

the nightside. Robustly constraining the climate will require more sophisticated GCMs and higher-precision phase curves/ eclipse mapping (e.g., de Wit et al. 2012).

The GCMs also provide insight into what molecules are present in which parts of the atmosphere. As discussed in Section 5, water dissociation and H$^-$ opacity are needed to explain the dayside emission spectrum. Figure 16 shows the photospheric abundances of H$_2$O and H$^-$ compared to the predicted temperature for the $\tau_{\text{drag4}}$ GCM. The water abundance drops by $\sim 10$ at the substellar point, and the H$^-$ opacity increases by $\sim 100$. By contrast, CO remains intact throughout the atmosphere. Our observations are not precise enough to detect water features on the nightside of the planet (see Section 7), but future high-precision data may be sensitive to these features and could help constrain the strength of horizontal transport in the atmosphere (Agúndez et al. 2014).

## 6.1. Constraints on the Planet's Magnetic Field

We show in this section that the small observed hotspot offset in the phase curves is best fit by a GCM that includes Rayleigh drag with a timescale $\tau_{\text{drag}} = 10^4$ s. This observation gives rise to the question: what magnetic field strength on the planet can produce drag with this timescale? Previous efforts to characterize exoplanet magnetic fields have mainly focused on magnetic interaction between the planet and its host star (e.g., Wright & Miller 2015 and references therein) and planetary radio emission (Grießmeier 2015). A complementary approach is to study the effect of the magnetic field on the planet's atmospheric dynamics.

Here, we make a simple order-of-magnitude estimate for the magnetic field required to produce a drag timescale of order $10^4$ s. We first computed the free electron abundance due to ionized metals for the $\tau_{\text{drag4}}$ GCM. Using the analytic expression from Perna et al. (2010) (Equation (12)), we
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Row (1): photospheric temperatures for the τdrag4 GCM for different viewing geometries. The orbital phase α = 0 corresponds to a secondary eclipse (when the substellar point faces Earth) and α = 180 corresponds to a transit (when the antistellar point faces Earth). Rows (2–4) show the abundances of CO, H₂O, and H⁻. Water dissociates at dayside temperatures, so the photospheric water abundance drops by ~10 and the H⁻ abundance increases by ~100. By contrast, the CO abundance is uniform throughout the photosphere.</figcaption></figure>

estimated a drag timescale at the substellar point of τdrag = 2 × 10⁴/B² s, where B is the magnetic field in Gauss. We assumed a temperature and pressure of 3359 K and 0.11 bar, and that the magnetic field is perpendicular to the flow. To reach a drag timescale of 10⁴ s, a magnetic field stronger than ~1 Gauss is required, comparable to Jupiter's magnetic field strength of 5–10 Gauss (Bagenal et al. 2004). To confirm this intriguing result, more detailed study is warranted, including a full magnetic hydrodynamic simulation of the atmospheric dynamics (e.g., Rogers 2017) that accounts for the possibility of complex magnetic field structure due to interactions between the magnetic fields of the planet and the star.

7. Comparison with Brown Dwarfs and Directly Imaged
Companions

WASP-103b is so highly irradiated that its photospheric
temperature (2000–3000 K) is comparable to that of low-mass
stars. However, the planet’s other properties (surface gravity,
rotation rate, irradiation) are different. To explore the effects of
varying these parameters, we selected spectra from WASP-
103b at three orbital phases: dayside (φ = 0.5), quadrature
(φ = 0.25), and nightside (φ = 0.1), and compared them to
brown dwarfs and young, directly imaged companions with
comparable brightness temperatures.

We also used three brown dwarfs/low-mass stars for
comparison. We chose these field sources: 2MASS
J13204427+0409045, (1320+0409), an optical L3; 2MASS
J04285096-2253227 (0428-2253), an optical L0.5; and
2MASS J00034227-2822410 (0003-2822), an optical M8
(see Filippazzo et al. 2015). We then used all currently
available photometric, astrometric, and spectroscopic data for
each source to evaluate fundamental parameters such as mass,
T<sub>eff</sub>, and log g (Filippazzo et al. 2015; Faherty et al. 2016) and
create flux-calibrated spectral energy distributions. For 1320
+0409, we used SDSS, WISE, and 2MASS photometry, along
with the optical spectrum from Reid et al. (2008), the near-

infrared spectrum from Bardalez Gagliuffi et al. (2014), and the
parallax reported in Faherty et al. (2012). For 0428–2253, we
used 2MASS, DENIS, and WISE photometry, along with the
optical spectrum from Kendall et al. (2003) the near-infrared
spectrum from Bardalez Gagliuffi et al. (2014) and the parallax
reported in Dieterich et al. (2014). For 0003–2822, we used
2MASS and WISE photometry, along with the optical
spectrum from Cruz et al. (2007), the near-infrared spectrum
from Cruz et al. (2018), and the parallax reported in Faherty
et al. (2010). Coincidentally, both 1320+0409 and 0003–2822
are widely separated (>2000 au) companions to a K2 and a G8
star, respectively. All data were gathered from the Brown
Dwarfs in New York City (BDNYC) database (Filippazzo
et al. 2015).<fnref n="18" /> At the assumed field ages of each source, 0003-
2822 would be above the nuclear burning boundary (star),
while 1320 + 0409 would be below (brown dwarf). It is likely
that 0428–2253 is a star, but it could be a brown dwarf at a
slightly younger field age.

The directly imaged spectra are for the sources CD-35 2722, USco 1610-1913B, and TWA 22A, and are taken from Wahhaj et al. (2011), Aller et al. (2013), Bonnefoy et al. (2014). They are young objects (aged 10-100 Myr), with lower surface gravities than brown dwarfs of comparable temperature. They also have gravitationally bound companions over a wide range of separations (67 ± 4, 840 ± 90, and 1.8 ± 0.1 au, respectively). The sources are calibrated in absolute flux using published H-band photometry (Bonnefoy et al. 2009; Wahhaj et al. 2011; Aller et al. 2013) and distances (Teixeira et al. 2009; Gaia Collaboration et al. 2016), a flux-calibrated spectrum from Vega (Hayes 1985; Mountain et al. 1985), and the corresponding filter passbands.

The system properties for all the objects are summarized in Table 7. Figure 17 shows the flux-calibrated spectra (assuming a distance of 10 pc). We compared the spectra over a

<footnote n="18">¹⁸ The BDNYC Database: http://database.bdnyc.org/.</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Flux-calibrated spectra for WASP-103b (left column) compared to observed spectra for brown dwarfs (middle) and directly imaged companions (right), assuming a distance of 10 pc. The WASP-103b spectra are from the nightside (phase 0.1, top row), quadrature (phase 0.25, middle), and the dayside (phase 0.5, bottom). Each row shows spectra from objects of comparable temperature. The dotted gray line corresponds to the best-fit blackbody. Effective temperatures are listed in the upper right corners. The wavelength bins used to calculate water feature amplitude $A_1$ are shaded in gray. Note that the panels do not have the same y-scale, because the objects have different radii, which leads to large variation in absolute flux.</figcaption></figure>

<table><caption>**Table 7**
Source Properties</caption><thead><tr><th></th><th>Object</th><th>T<sub>eff</sub> (K)</th><th>log g (cgs)</th><th>H<sub>2</sub>O A<sub>1</sub></th><th>H<sub>2</sub>O A<sub>2</sub></th></tr></thead><tbody><tr><td rowspan="3">Hot Jupiter</td><td>W103b night</td><td>1880 ± 40</td><td>3.2 ± 0.04</td><td>0.07 ± 1.8e-01</td><td>-0.00 ± 1.6e-01</td></tr><tr><td>W103b quadrature</td><td>2400 ± 40</td><td>3.2 ± 0.04</td><td>-0.14 ± 7.8e-02</td><td>-0.01 ± 6.7e-02</td></tr><tr><td>W103b dayside</td><td>2930 ± 40</td><td>3.2 ± 0.04</td><td>0.04 ± 1.4e-02</td><td>0.15 ± 1.2e-02</td></tr><tr><td rowspan="3">Brown Dwarf</td><td>2MASS J1320+0409</td><td>1880 ± 70</td><td>5.19 ± 0.16</td><td>0.21 ± 6.3e-04</td><td>-0.06 ± 5.0e-04</td></tr><tr><td>2MASS J0428-2253</td><td>2430 ± 80</td><td>5.22 ± 0.09</td><td>0.16 ± 1.2e-04</td><td>-0.03 ± 1.0e-04</td></tr><tr><td>2MASS J0003-2822</td><td>2890 ± 80</td><td>5.18 ± 0.04</td><td>0.26 ± 1.3e-04</td><td>0.10 ± 1.1e-04</td></tr><tr><td rowspan="3">Imaged Companion</td><td>CD-35 2722</td><td>1800 ± 100</td><td>4.5 ± 0.5</td><td>...</td><td>0.15 ± 1.0e-05</td></tr><tr><td>USco 1610-1913B</td><td>2400 ± 150</td><td>...</td><td>0.27 ± 2.0e-04</td><td>0.19 ± 2.2e-04</td></tr><tr><td>TWA 22A</td><td>3000 ± 100</td><td>4.5 ± 0.5</td><td>...</td><td>0.29 ± 2.0e-05</td></tr></tbody></table>

wavelength range of 1.1–1.7 µm. The most prominent spectral features over this range are expected to come from water, which has a forest of absorption lines near 1.4 µm. Spectra for the brown dwarfs and imaged companions have noticeable features in the water band, whereas WASP-103b does not. To quantitatively compare the water feature amplitude for different objects, we define an amplitude $A = (F_{1,2} - F_{3,4})/F_{1,2}$, where $F_{1,2}$ is the weighted mean flux in a wavelength bin $\lambda_1 - \lambda_2$. We calculated the water feature amplitude for two choices of wavelength bins. For the first, $A_1$, we considered data in and out of the water band, with $\lambda_{1,2,3,4} = 1.15, 1.3, 1.35$, and 1.5 µm. The ground-based direct imaging data do not span this entire wavelength range, so we define another amplitude, $A_2$, with $\lambda_{1,2,3,4} = 1.2, 1.35, 1.5$, and 1.65 µm. The estimated amplitudes and uncertainties are listed in Table 7. We note that a number of indices have been defined to characterize features in brown dwarf spectra, and these indices have revealed trends in the amplitude of a range of spectral features (water, sodium,

potassium, VO, and FeH) with surface gravity and effective temperature (Reid et al. 2001; Geballe et al. 2002; McLean et al. 2003). However, the WASP-103b data do not have high enough signal-to-noise or spectral resolution to meaningfully compare these indices. Instead we use a broader bandpass to address a simpler question: are those spectra with the same temperatures consistent with each other?

We find that the $A_1$ and $A_2$ values are significantly lower for WASP-103b at dayside and quadrature than for brown dwarfs and imaged companions of similar temperature. WASP-103b typically has $A_1$ and $A_2$ consistent with zero, indicating no water absorption (in agreement with the analysis in Section 4 that showed water is depleted in the photosphere). By contrast, the brown dwarfs and young companions have significant water features, with drops in flux of about 20% in the water band. This is not surprising: stars in the temperature range 2000–3000 K have well-known, prominent water features (Kirkpatrick et al. 1993). Based on the grid retrieval of
WASP-103b's atmospheric composition, there are several reasons WASP-103b would exhibit different behavior at the same temperature. WASP-103b is irradiated from above rather than below, changing the shape of the temperature-pressure profile. In addition, WASP-103b also has much lower surface gravity (logg = 3.2 versus 4–5 for stars), which pushes the photosphere to lower pressures, where water dissociates more readily (Arcangeli et al. 2018). These factors are not relevant on the nightside, and 3D models predict that WASP-103b has nightside water absorption features; however, the current data are not precise enough to distinguish between a blackbody spectrum versus water features like those seen in the other objects.

## 8. Summary

We observed thermal phase curves of the hot Jupiter WASP-103b measured with HST/WFC3 time series spectroscopy (1.15–1.65 µm) and Spitzer/IRAC broadband photometry (3.6 and 4.6 µm bands). Here, we summarize our conclusions about the atmosphere based on these measurements.

1. The dayside planet-to-star flux is 0.151% ± 0.015%, 0.446% ± 0.38%, and 0.569% ± 0.014% in the WFC3 bandpass, Spitzer 3.6, and Spitzer 4.5 µm, respectively. The best-fit blackbody to the WFC3 dayside spectrum has a brightness temperature of 2930 ± 40 K, making WASP-103b among the hottest exoplanets ever observed.

2. The phase curves have large amplitudes (0.8–0.9× the secondary eclipse depth), and small offsets in peak brightness from the substellar point (consistent with zero degrees at all wavelengths). These characteristics indicate inefficient redistribution of heat to the nightside, as seen in other very hot Jupiters (Komacek et al. 2017).

3. We fit the phase variation with the spiderman package (Louden & Kreidberg 2017) to evaluate different models of the planet's climate, including a two-temperature map, a physically motivated kinematic map, and spherical harmonics. The spherical harmonic temperature map generally provides the best fit to the data; however, all the maps produce reasonable fits ($\chi^2$ near unity), and there are large differences in temperature between them (up to 1000 K at a given latitude/longitude). Breaking the degeneracy between different climate maps will require higher-precision phase curves and/or secondary eclipse mapping (e.g., de Wit et al. 2012).

4. We calculated phase-resolved spectra in ten orbital phase bins. The HST/WFC3 spectra are consistent with black-body emission from the planet at all orbital phases. The best-fit brightness temperatures ranges from 1880 ± 40 K (phase $\phi$ = 0.1) to 2930 ± 40 K on the dayside. We attribute the absence of water features at WFC3 wavelengths to (1) H$_2$O dissociation on the dayside and (2) additional near-IR opacity from H$^-$, TiO/VO, and FeH.

5. The Spitzer data are not consistent with the best-fit blackbody to the WFC3 data: they have a higher brightness temperature at phases $\phi$ = 0.2–0.5, which transitions to a lower brightness temperature on the nightside ($\phi$ = −0.2–0.1). An atmospheric retrieval analysis suggests that these characteristics are likely due to CO features in the infrared, as well as to a temperature inversion on the dayside but not the nightside.

6. The measured transmission spectrum is featureless (after correcting for nightside emission from the planet). A 3D model predicts water features in transmission that could be detected with future high-precision observations.

7. We characterized the composition with a 1D grid-based retrieval that assumes thermochemical and radiative-convective equilibrium. The atmosphere is moderately metal-enriched ($23^{+29}_{-13}$× solar; and $>1$× solar at $3\sigma$ confidence). This value is somewhat higher than what is observed for other gas giants (e.g., Wong et al. 2004; Kreidberg et al. 2014a), but may be indicative of intrinsic scatter in the relationship between atmospheric metallicity and planet mass predicted by theoretical models (Fortney et al. 2013; Mordasini et al. 2016). However, the metallicity is strongly sensitive to the 4.5 µm Spitzer eclipse depth, and additional observations would be useful in confirming the metal enhancement. In addition to metallicity, we also infer an upper limit on the carbon-to-oxygen ratio of 0.9 ($3\sigma$ confidence). This estimate agrees with expectations from planet formation models, i.e., that pollution from water ice in planetesimals leads to C/O < 1 in gas giant atmospheres (Mordasini et al. 2016; Espinoza et al. 2017). The best-fit temperature pressure profile has a thermal inversion from ~$10^{-2}$–$10^{-3}$ bars due to TiO/VO absorption at high altitudes.

8. We ran several 3D GCMs to compare to the data, including a nominal model with a cloud-free solar composition, a metal-enriched model ([Fe/H] = 0.5), and two models with Lorentz drag. The GCM with a Lorentz drag timescale of $10^4$ s matches the data best. This model has an equator-to-pole wind pattern that reproduces the small observed hotspot offsets and large phase curve amplitudes. We made a simple order-of-magnitude estimate for the magnetic field strength required to produce this fast drag timescale, and found that it implies a magnetic field of ~1 Gauss.

9. We compared the spectra of WASP-103b at phases 0.5 (dayside), 0.25 (quadrature), and 0.1 (nightside) to brown dwarfs and directly imaged companions of similar temperature. We quantify the strength of the water feature and find that both brown dwarfs and imaged companions show evidence for water absorption at 1.4 µm, whereas the WASP-103b dayside and quadrature spectra do not. We attribute the difference to two factors: WASP-103b's irradiation environment, which changes the temperature pressure profile, and its low surface gravity, which pushes the photosphere to higher altitudes where water dissociates more easily. The WASP-103b nightside spectra have larger uncertainties and are consistent with the water feature amplitudes for other objects; higher-precision phase curves are needed to detect water on the nightside.

These results provide a first look at the global composition and thermal structure of WASP-103b. The planet is complex, showing changes in temperature profile with longitude, possible gradients in composition from dayside to nightside, and circulation patterns that may be influenced by the magnetic field. These findings highlight the 3D nature of exoplanets and illustrate the importance of phase curve observations to develop a comprehensive understanding of their atmospheric chemistry and physics.
Support for HST program GO-15050 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Support for Spitzer program 11099 was provided by NASA through an award issued by JPL/Caltech. The data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) and the NASA/IPAC Infrared Science Archive. The Infrared Science Archive is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work also made use of the Python packages SciPy and NumPy (Jones et al. 2001; Van Der Walt et al. 2011). The authors are grateful for helpful conversations with Caroline Morley, Thomas Beatty, Ming Zhao, Kimberly Cartier, Hannah Diamond-Lowe, and Nick Cowan. We also thank the organizers of the 2016 Santa Cruz Kavli Summer Program and the 2017 Ringberg Atmospheres of Disks and Planets meeting for facilitating productive discussion and collaboration. J.L.B. acknowledges support from the David and Lucile Packard Foundation. J.M.D. acknowledges support from the European Research Council (ERC) under the programme Exo-Atmos (grant agreement no. 679633). G.W.H. and M.H.W. also acknowledge long-term support from Tennessee State University and the State of Tennessee through its Centers of Excellence Program. M.R.L. acknowledges NASA XRP grant NNX17AB56G for partial support of the theoretical interpretation of the data, as well as the ASU Research Computing staff for support with the Saguaro and Agave compute clusters. Finally, we thank the anonymous referee for a thoughtful and detailed report that improved the quality of the paper.

## ORCID iDs

Laura Kreidberg https://orcid.org/0000-0003-0514-1147

Michael R. Line https://orcid.org/0000-0002-2338-476X

Vivien Parmentier https://orcid.org/0000-0001-9521-6258

Kevin B. Stevenson https://orcid.org/0000-0002-7352-7941

Tom Loudon https://orcid.org/0000-0002-2574-356X

Jacqueline K. Faherty https://orcid.org/0000-0001-6251-0573

Gregory W. Henry https://orcid.org/0000-0003-4155-8513

Keivan Stassun https://orcid.org/0000-0002-3481-9052

Jacob L. Bean https://orcid.org/0000-0003-4733-6532

Jonathan J. Fortney https://orcid.org/0000-0002-9843-4354

## References

Adcroft, A., Campin, J.-M., Hill, C., & Marshall, J. 2004, MWRv, 132, 2845
Agúndez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, A&A,
564, A73

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

Aller, K. M., Kraus, A. L., Liu, M. C., et al. 2013, ApJ, 773, 63

Arcangeli, J., Desert, J.-M., Line, M. R., et al. 2018, arXiv:1801.02489

Armstrong, D. J., de Mooij, E., Barstow, J., et al. 2016, NatAs, 1, 0004

Bagenal, F., Dowling, T. E., & McKinnon, W. B. 2004, Jupiter: The Planet,
Satellites and Magnetosphere (Cambridge: Cambridge Univ. Press)

Bardalez Gagliuffi, D. C., Burgasser, A. J., Gelino, C. R., et al. 2014, ApJ,
794, 143

Batygin, K., Stanley, S., & Stevenson, D. J. 2013, ApJ, 776, 53

Beatty, T. G., Madhusudhan, N., Pogge, R., et al. 2017a, AJ, 154, 242

Beatty, T. G., Madhusudhan, N., Tsiaras, A., et al. 2017b, AJ, 154, 158

Bell, T. J., & Cowan, N. B. 2018, arXiv:1802.07725

Blecic, J., Harrington, J., Madhusudhan, N., et al. 2013, ApJ, 779, 5

Blecic, J., Harrington, J., Madhusudhan, N., et al. 2014, ApJ, 781, 116

Bonnefoy, M., Chauvin, G., Dumas, C., et al. 2009, A&A, 506, 799

Bonnefoy, M., Chauvin, G., Lagrange, A.-M., et al. 2014, A&A, 562, A127

Borucki, W. J., Koch, D., Jenkins, J., et al. 2009, Sci, 325, 709

Braak, C. J. F. T. 2006, Statistics and Computing, 16, 239

Campo, C. J., Harrington, J., Hardy, R. A., et al. 2011, ApJ, 727, 125

Carter, J. A., & Winn, J. N. 2009, ApJ, 704, 51

Cartier, K. M. S., Beatty, T. G., Zhao, M., et al. 2017, AJ, 153, 34

Cowan, N. B., & Agol, E. 2008, ApJL, 678, L129

Cowan, N. B., & Fujii, Y. 2017, arXiv:1704.07832

Crossfield, I. J. M., Barman, T., Hansen, B. M. S., Tanaka, I., & Kodama, T. 2012, ApJ, 760, 140

Cruz, K. L., Núñez, A., Burgasser, A. J., et al. 2018, AJ, 155, 34

Cruz, K. L., Reid, I. N., Kirkpatrick, J. D., et al. 2007, AJ, 133, 439

Cubillos, P., Harrington, J., Madhusudhan, N., et al. 2013, ApJ, 768, 42

de Wit, J., Gillon, M., Demory, B.-O., & Seager, S. 2012, A&A, 548, A128

Delrez, L., Madhusudhan, N., Lendl, M., et al. 2018, MNRAS, 474, 2334

Demory, B.-O., de Wit, J., Lewis, N., et al. 2013, ApJL, 776, L25

Diamond-Lowe, H., Stevenson, K. B., Bean, J. L., Line, M. R., & Fortney, J. J. 2014, ApJ, 796, 66

Dieterich, S. B., Henry, T. J., Jao, W.-C., et al. 2014, AJ, 147, 94

Espinoza, N., Fortney, J. J., Miguel, Y., Thorngren, D., & Murray-Clay, R. 2017, ApJL, 838, L9

Evans, T. M., Sing, D. K., Wakeford, H. R., et al. 2016, ApJL, 822, L4

Faherty, J. K., Burgasser, A. J., Walter, F. M., et al. 2012, ApJ, 752, 56

Faherty, J. K., Burgasser, A. J., West, A. A., et al. 2010, AJ, 139, 176

Faherty, J. K., Riedel, A. R., Cruz, K. L., et al. 2016, ApJS, 225, 10

Filippazzo, J. C., Rice, E. L., Faherty, J., et al. 2015, ApJ, 810, 158

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419

Fortney, J. J., Mordasini, C., Nettelmann, N., et al. 2013, ApJ, 775, 80

Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Geballe, T. R., Knapp, G. R., Leggett, S. K., et al. 2002, ApJ, 564, 466

Gillon, M., Anderson, D. R., Collier-Cameron, A., et al. 2014, A&A, 562, L3

Gordon, S., & McBride, B. J. 1994, Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications, Vol. 1 (Washington, DC: NASA)

Grießmeier, J.-M. 2015, in Characterizing Stellar and Exoplanetary Environments, ed. H. Lammer & M. Khodachenko (Berlin: Springer), 213

Harrington, J., Hansen, B. M., Luszcz, S. H., et al. 2006, Sci, 314, 623

Harrington, J., Luszcz, S., Seager, S., Deming, D., & Richardson, L. J. 2007, Natur, 447, 691

Hayes, D. S. 1985, in IAU Symp. 111, Calibration of Fundamental Stellar Quantities, ed. D. S. Hayes, L. E. Pasinetti, & A. G. D. Philip (Dordrecht: Reidel), 225

Haynes, K., Mandell, A. M., Madhusudhan, N., Deming, D., & Knutson, H. 2015, ApJ, 806, 146

Henry, G. W. 1999, PASP, 111, 845

Horne, K. 1986, PASP, 98, 609

Hu, R., Demory, B.-O., Seager, S., Lewis, N., & Showman, A. P. 2015, ApJ,
802, 51

Hubeny, I., Burrows, A., & Sudarsky, D. 2003, ApJ, 594, 1011

Ingalls, J. G., Krick, J. E., Carey, S. J., et al. 2016, AJ, 152, 44

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Kataria, T., Sing, D. K., Lewis, N. K., et al. 2016, ApJ, 821, 9

Kendall, T. R., Mauron, N., Azzopardi, M., & Gigoyan, K. 2003, A&A,
403, 929

Kirkpatrick, J. D., Kelly, D. M., Rieke, G. H., et al. 1993, ApJ, 402, 643

Koll, D. D. B., & Komacek, T. D. 2018, ApJ, 853, 133

Komacek, T. D., & Showman, A. P. 2016, ApJ, 821, 16

Komacek, T. D., Showman, A. P., & Tan, X. 2017, ApJ, 835, 198

Kreidberg, L. 2015, PASP, 127, 1161

Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014a, ApJL, 793, L27

Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014b, Natur, 505, 69

Kreidberg, L., Line, M. R., Bean, J. L., et al. 2015, ApJ, 814, 66

Leconte, J., Lai, D., & Chabrier, G. 2011a, A&A, 536, C1

Leconte, J., Lai, D., & Chabrier, G. 2011b, A&A, 528, A41

Lendl, M., Cubillos, P. E., Hagelberg, J., et al. 2017, A&A, 606, A18

Line, M. R., Knutson, H., Wolf, A. S., & Yung, Y. L. 2014, ApJ, 783, 70

Line, M. R., Stevenson, K. B., Bean, J., et al. 2016, AJ, 152, 203

Line, M. R., Wolf, A. S., Zhang, X., et al. 2013, ApJ, 775, 137

Lodders, K., & Fegley, B. 2002, Icar, 155, 393
Loeb, A., & Gaudí, B. S. 2003, ApJL, 588, L17

Long, K. S., Baggett, S. M., & MacKenty, J. W. 2015, Persistence in the WFC3 IR Detector: an Improved Model Incorporating the Effects of Exposure Time, Tech. rep.

Lothringer, J. D., Barman, T., & Koskinen, T. 2018, arXiv:1805.00038

Louden, T., & Kreidberg, L. 2017, arXiv:1711.00494

Lust, N. B., Britt, D., Harrington, J., et al. 2014, PASP, 126, 1092

Madhusudhan, N., Harrington, J., Stevenson, K. B., et al. 2011, Natur, 469, 64

Mansfield, M., Bean, J. L., Line, M. R., et al. 2018, arXiv:1805.00424

Marley, M. S., & McKay, C. P. 1999, Icar, 138, 268

McCullough, P., & MacKenty, J. 2012, Considerations for Using Spatial Scans with WFC3, Tech. Rep. 2012-08 (Baltimore, MD: STScI)

McLean, I. S., McGovern, M. R., Burgasser, A. J., et al. 2003, ApJ, 596, 561

Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41

Mountain, C. M., Selby, M. J., Leggett, S. K., Blackwell, D. E., & Petford, A. D. 1985, A&A, 151, 399

Nymeyer, S., Harrington, J., Hardy, R. A., et al. 2011, ApJ, 742, 35

Parmentier, V., & Crossfield, I. 2017, arXiv:1711.07696

Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22

Parmentier, V., Line, M. R., Bean, J. L., et al. 2018, arXiv:1805.00096

Perez-Becker, D., & Showman, A. P. 2013, ApJ, 776, 134

Perna, R., Menou, K., & Rauscher, E. 2010, ApJ, 719, 1421

Pont, F., Zucker, S., & Queloz, D. 2006, MNRAS, 373, 231

Reid, I. N., Burgasser, A. J., Cruz, K. L., Kirkpatrick, J. D., & Gizis, J. E. 2001, AJ, 121, 1710

Reid, I. N., Cruz, K. L., Kirkpatrick, J. D., et al. 2008, AJ, 136, 1290

Rogers, T. M. 2017, NatAs, 1, 0131

Rogers, T. M., & Komacek, T. D. 2014, ApJ, 794, 132

Rogers, T. M., & Showman, A. P. 2014, ApJL, 782, L4

Sheppard, K. B., Mandell, A. M., Tamburo, P., et al. 2017, ApJL, 850, L32

Showman, A. P., Fortney, J. J., Lewis, N. K., & Shabram, M. 2013, ApJ, 762, 24

Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564

Showman, A. P., & Guillot, T. 2002, A&A, 385, 166

Southworth, J., Mancini, L., Ciceri, S., et al. 2015, MNRAS, 447, 711

Stevenson, K. B., Bean, J. L., Madhusudhan, N., & Harrington, J. 2014a, ApJ, 791, 36

Stevenson, K. B., Bean, J. L., Seifahrt, A., et al. 2014b, AJ, 147, 161

Stevenson, K. B., Bean, J. L., Seifahrt, A., et al. 2016, ApJ, 817, 141

Stevenson, K. B., Désert, J.-M., Line, M. R., et al. 2014c, Sci, 346, 838

Stevenson, K. B., Harrington, J., Fortney, J. J., et al. 2012a, ApJ, 754, 136

Stevenson, K. B., Harrington, J., Fortney, J. J., et al. 2012b, ApJ, 754, 136

Stevenson, K. B., Harrington, J., Lust, N. B., Lewis, N. K., et al. 2012c, ApJ, 755, 9

Stevenson, K. B., Harrington, J., Nymeyer, S., Madhusudhan, N., et al. 2010, Natur, 464, 1161

Stevenson, K. B., Line, M. R., Bean, J. L., et al. 2017, AJ, 153, 68

Teixeira, R., Ducourant, C., Chauvin, G., et al. 2009, A&A, 503, 281

Tinetti, G., Tennyson, J., Griffith, C. A., & Waldmann, I. 2012, RSPTA, 370, 2749

Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

Visscher, C., Lodders, K., & Fegley, B., Jr. 2006, ApJ, 648, 1181

Wahhaj, Z., Liu, M. C., Biller, B. A., et al. 2011, ApJ, 729, 139

Wöllert, M., & Brandner, W. 2015, A&A, 579, A129

Wong, M. H., Mahaffy, P. R., Atreya, S. K., Niemann, H. B., & Owen, T. C. 2004, Icar, 171, 153

Wright, J. T., & Miller, B. 2015, IAUGA, 22, 2258453

Zellem, R. T., Swain, M. R., Roudier, G., et al. 2017, ApJ, 844, 27

Zhang, M., Knutson, H. A., Kataria, T., et al. 2018, AJ, 155, 83

Zhang, X., & Showman, A. P. 2017, ApJ, 836, 73

Zhou, Y., Apai, D., Lew, B. W. P., & Schneider, G. 2017, AJ, 153, 243