# An Independent Analysis of the Spitzer/IRAC Phase Curves of WASP43 b

Giuseppe Morello, Camilla Daniëlski, D. Dickens, Pascal Tremblin, Pierre-Olivier Lagage

► To cite this version:

Giuseppe Morello, Camilla Daniëlski, D. Dickens, Pascal Tremblin, Pierre-Olivier Lagage. An Independent Analysis of the Spitzer/IRAC Phase Curves of WASP43 b. The Astronomical Journal, 2019, 157 (5), pp.205. 10.3847/1538-3881/ab14e2 . hal-04488791

HAL Id: hal-04488791

https://hal.science/hal-04488791v1

Submitted on 5 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

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

CrossMark

# An Independent Analysis of the *Spitzer*/IRAC Phase Curves of WASP43 b

G. Morello¹, C. Daniëlski¹, D. Dickens¹, P. Tremblin², and P.-O. Lagage¹

¹AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France; giuseppe.morello@cea.fr

²Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

Received 2018 October 17; revised 2019 March 4; accepted 2019 March 28; published 2019 May 1

## Abstract

We present here a reanalysis of the *Spitzer Space Telescope* phase curves of the hot Jupiter WASP43 b, using the wavelet pixel-independent component analysis, a blind signal-source separation method. The data analyzed were recorded with the Infrared Array Camera and consisted of two visits at 3.6 µm, and one visit at 4.5 µm, each visit covering one transit and two eclipse events. To test the robustness of our technique we repeated the analysis on smaller portions of the phase curves, and by employing different instrument ramp models. Our reanalysis presents significant updates of the planetary parameters compared to those reported in the original phase curve study of WASP43 b. In particular, we found (1) higher nightside temperatures, (2) smaller hotspot offsets, (3) a greater consistency (~1σ) between the two 3.6 µm visits, and (4) a greater similarity with the predictions of the atmospheric circulation models. Our parameter results are consistent within 1σ with those reported by a recent reanalysis of the same data sets. For each visit we studied the variation of the retrieved transit parameters as a function of various sets of stellar limb-darkening coefficients, finding significant degeneracy between the limb-darkening models and the analysis output. Furthermore, we performed the analysis of the single transit and eclipse events, and we examined the differences between these results with the ones obtained with the whole phase curve. Finally we provide a formula useful to optimize the trade-off between precision and duration of observations of transiting exoplanets.

**Key words:** planets and satellites: atmospheres – planets and satellites: fundamental parameters – planets and satellites: individual (WASP43 b) – stars: atmospheres – techniques: photometric – techniques: spectroscopic

## 1. Introduction

WASP43 b is a hot Jupiter orbiting around a K7 V star in ~19.5 hr (Hellier et al. 2011). Table 1 reports the stellar, planetary, and transit parameters taken from the WASP43 b discovery paper (Hellier et al. 2011). The ultra-short orbital period of WASP43 b has inspired multiple observational programs of its full phase curve using the *Hubble Space Telescope* (HST; Stevenson et al. 2014), and the *Spitzer Space Telescope* (Stevenson et al. 2017, hereafter S17). Another full phase curve observation of WASP43 b is planned as part of the Transiting Exoplanet Community Early Release Science Program of the *James Webb Space Telescope* (JWST) using the Mid-infrared Instrument (Bean et al. 2018).

Exoplanet phase curves are measurements of the flux coming from a star+exoplanet system as a function of the orbital phase. If the exoplanet is transiting, its phase curve includes (usually) both transit and eclipse events. The flux modulations observed in the mid-infrared are attributed to the thermal emission from the exoplanet with varying phase angle (Cooper & Showman 2005; Fortney et al. 2006; Cowan et al. 2007). Exoplanets with short orbital periods are expected to be tidally locked to their host star (Showman & Guillot 2002), therefore exhibiting a hotter dayside and a cooler nightside. The day-night temperature contrast depends on the heat recirculation efficiency of the exoplanetary atmosphere. Numerical simulations also predict a (model-dependent) hotspot offset from the substellar point (Showman & Guillot 2002; Cooper & Showman 2005; Kataria et al. 2015; Schwartz et al. 2017; Zhang & Showman 2017).

Stevenson et al. (2014) and S17 claimed extremely low circulation efficiency for the atmosphere of WASP43 b: $ε = 0.002^{+0.01}_{-0.002}$, where $ε$ is the night-day bolometric flux ratio. They also detected a wavelength-dependent eastward

hotspot offset, i.e., their phase curve models peak prior to secondary eclipses. However, S17 discarded the first 3.6 µm data set, which presented discrepant results, and larger correlated noise in the light curve residuals. S17 also discarded a ~2 hr interval from the second 3.6 µm data set, corresponding to an unexpected flux decrement in their detrendend light curve, that the authors attributed to unmodeled instrumental or astrophysical red noise.

The low nightside fluxes and the large eastward offsets measured by S17 in the *Spitzer*/Infrared Array Camera (IRAC) passbands could not be reproduced using the SPARC/MITgcm code of Kataria et al. (2015). The SPARC/MITgcm is a 3D global circulation model coupled to a non-gray radiative transfer code. Keating & Cowan (2017) pointed out that the atmosphere of WASP43 b should have a much higher circulation efficiency, $ε \sim 0.5$, based on the inverse correlation between the day-night temperature contrast and stellar irradiation (Cowan & Agol 2011; Perez-Becker & Showman 2013; Schwartz & Cowan 2015; Komacek & Showman 2016).

Mendonça et al. (2018a hereafter M18) reanalyzed the three *Spitzer*/IRAC phase curves without discarding any data. M18 found a better agreement between the two 3.6 µm observations, and higher nightside fluxes than those reported by S17.

In this paper we present an independent reanalysis of the three *Spitzer*/IRAC phase curves of WASP43 b using the wavelet pixel-independent component analysis (ICA) pipeline (Morello et al. 2016). We repeated the analysis by adopting different stellar limb-darkening models, which affect the retrieved transit parameters. We compare our results with those reported by S17 and M18, and with theoretical expectations. In addition to the full phase curve analyses, we explore the ability to constrain the different parameters from shorter observations, nominally half phase curves, transit-only,
<table><caption>**Table 1**
WASP43 System Parameters</caption><thead><tr><th></th><th>Stellar Parameters</th></tr></thead><tbody><tr><td>T<sub>eff</sub> (K)</td><td>4400 ± 200</td></tr><tr><td>log g<sub>*</sub> (cgs)</td><td>4.65<sup>+0.06</sup><sub>-0.04</sub></td></tr><tr><td>[Fe/H] (dex)</td><td>-0.05 ± 0.17</td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>0.58 ± 0.05</td></tr><tr><td>R<sub>*</sub> (R<sub>☉</sub>)</td><td>0.60<sup>+0.03</sup><sub>-0.04</sub></td></tr></tbody><thead><tr><th></th><th>Planetary Parameters</th></tr></thead><tbody><tr><td>M<sub>p</sub> (M<sub>Jup</sub>)</td><td>1.78 ± 0.10</td></tr><tr><td>R<sub>p</sub> (R<sub>Jup</sub>)</td><td>0.93<sup>+0.07</sup><sub>-0.09</sub></td></tr><tr><td>a (au)</td><td>0.0142 ± 0.0004</td></tr></tbody><thead><tr><th></th><th>Transit Parameters</th></tr></thead><tbody><tr><td>p<sup>2</sup></td><td>0.0255 ± 0.0012</td></tr><tr><td>b</td><td>0.66<sup>+0.04</sup><sub>-0.07</sub></td></tr><tr><td>i (deg)</td><td>82.6<sup>+1.3</sup><sub>-0.9</sub></td></tr><tr><td>P (days)</td><td>0.813475 ± 0.000001</td></tr><tr><td>E.T. (HJD)</td><td>2455528.86774 ± 0.00014</td></tr></tbody></table>

Note. From Hellier et al. (2011).

and eclipse-only. This kind of study will be useful for planning future JWST proposals, and optimizing the time schedule of the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) mission, in order to maximize their scientific return.

## 2. Observations

We reanalyzed three Spitzer/IRAC observations of the phase curve of WASP43 b. Each visit consists of two to three consecutive astronomical observation requests (AORs) over a 25.4 hr interval, including one transit and two eclipse events. Observational and detector information for the individual data sets is given in Table 2.

## 3. Data Analysis

### 3.1. The Phase Curve Model

In our model the stellar flux is constant in time, and normalized to 1. The exoplanetary flux is given by

$$c_0 + c_1 \cos[2\pi(\Phi - \Delta\Phi - c_2)] + c_3 \cos[4\pi(\Phi - \Delta\Phi - c_4)], \quad (1)$$

where $\Phi$ is the so-called orbital phase, i.e., the time from the reference epoch of transit (E.T.) in units of the orbital period (P), $\Delta\Phi$ is the mid-transit phase offset, and $c_0-c_4$ are free parameters used to model the phase curve modulations. Equation (1) is equivalent to the formula adopted by S17 and M18. We used the formalism of Mandel & Agol (2002) for modeling the exoplanetary transit and eclipses.

### 3.2. Stellar Limb-darkening Coefficients

We calculated multiple sets of four limb-darkening coefficients (Claret 2000) hereafter claret-4, for the WASP43 star in the 3.6 and 4.5 µm Spitzer/IRAC passbands, using the code provided by Espinoza & Jordán (2015) at GitHub.<fnref n="3" /> The code

adopts two grids of stellar-atmosphere intensity models, i.e., ATLAS9⁴ (Kurucz 1979) and PHOENIX (Husser et al. 2013). The intensities in the models are given as a function of $\mu = \cos \theta$, where $\theta$ is the angle between the surface normal and the line of sight. The ATLAS models adopt a plane-parallel approximation for the stellar atmosphere, while the PHOENIX models use spherical geometry. As a consequence, the PHOENIX models show a characteristic steep drop-off in intensity at small, but finite $\mu$ values, which is not well approximated by any of the standard parametric laws (Claret et al. 2012, 2013; Morello et al. 2017). The limb-darkening coefficients also depend on the sampling of the intensities (Howarth 2011; Neilson & Lester 2013a, 2013b; Espinoza & Jordán 2015). We tested the following fitting options:

1. A17, i.e., least-squares fit to the ATLAS model intensities calculated at 17 angles;

2. A100, i.e., least-squares fit to the ATLAS intensities interpolated at 100 angles, uniformly sampled in $\mu$, with a cubic spline;

3. P100, i.e., least-squares fit to the PHOENIX intensities interpolated at 100 angles, uniformly sampled in $\mu$, with a cubic spline;

4. PQS, i.e., least-squares fit to the PHOENIX model intensities with $\mu \ge 0.1$ (quasi-spherical models, as defined by Claret et al. 2012).

We discarded the least-squares fit to all the PHOENIX model intensities, because it led to anomalous (non-monotonic) limb-darkening profiles. The most likely cause of the anomalous results was that the PHOENIX model intensities are more finely sampled near the steep drop-off, which is then overweighted in the fit. We interpolated the limb-darkening coefficients in $T_{eff}$ and log g to the WASP43 parameter values reported in Table 1. Table 3 reports the four sets of claret-4 limb-darkening coefficients obtained with the different fitting options. Figure 1 shows the corresponding intensity profiles. We note that the ATLAS limb-darkening profiles, A17 and A100, overlap in the plot. The PHOENIX profiles, P100 and PQS, indicate stronger limb-darkening than the ATLAS profiles. The P100 profiles reach zero intensity at the stellar limb, while the PQS profiles remain significantly above zero. Note that the PQS profiles are not accurate at the stellar limb, given that their behavior is extrapolated from the model intensities with $\mu \ge 0.1$.

### 3.3. Detrending Spitzer/IRAC Data

For our analysis we used the basic calibrated data (BCD) provided by the Spitzer Heritage Archive (Wu et al. 2010). BCD are flat-fielded and flux-calibrated frames (Fazio et al. 2004; IRAC Instrument & Instrument Support Teams 2015). We extracted the individual pixel time series from a 5 × 5 array having the stellar centroid at its center, and computed the sum-of-pixel time series, here referred to as raw light curves. We binned the time series by a factor of 8, i.e., temporal bin size of 16 s, in order to reduce the computational time for the data analysis. The chosen bin size is smaller than the time scales of interest, e.g., it is ~1/63 of the transit ingress duration.

Then, we applied the wavelet pixel-ICA technique (Morello et al. 2016) to simultaneously fit the phase curve model and the instrumental effects. We also tested the time pixel-ICA

<footnote n="3">³ http://www.github.com/nespinoza/limb-darkening/</footnote>

<footnote n="4">⁴ http://kurucz.harvard.edu/grids.html</footnote>
<table><caption>**Table 2**
Spitzer/IRAC Data Sets Analyzed for This Study</caption><thead><tr><th>Obs.<sup>a</sup></th><th>Prog. ID</th><th>AORS<sup>b</sup></th><th>UT Date</th><th>Δt (hr)<sup>c</sup></th><th>Mode<sup>d</sup></th><th>Pip.<sup>e</sup></th></tr></thead><tbody><tr><td rowspan="3">Ch1, visit 1<br>(3.6 μm)</td><td rowspan="3">11001</td><td>52364544</td><td>2015 Mar 7</td><td>8.5</td><td>sub, 2.0</td><td>19.2.0</td></tr><tr><td>52364800</td><td>2015 Mar 7</td><td>8.5</td><td>sub, 2.0</td><td>19.2.0</td></tr><tr><td>52355312</td><td>2015 Mar 8</td><td>8.5</td><td>sub, 2.0</td><td>19.2.0</td></tr><tr><td rowspan="2">Ch1, visit 2<br>(3.6 μm)</td><td rowspan="2">11001</td><td>57744384</td><td>2015 Sep 5</td><td>15.2</td><td>sub, 2.0</td><td>19.2.0</td></tr><tr><td>57744640</td><td>2015 Sep 5</td><td>10.2</td><td>sub, 2.0</td><td>19.2.0</td></tr><tr><td rowspan="3">Ch2<br>(4.5 μm)</td><td rowspan="3">10169</td><td>51777024</td><td>2014 Aug 27</td><td>8.5</td><td>sub, 2.0</td><td>19.2.0</td></tr><tr><td>51777280</td><td>2014 Aug 27</td><td>8.5</td><td>sub, 2.0</td><td>19.2.0</td></tr><tr><td>51777792</td><td>2014 Aug 28</td><td>8.5</td><td>sub, 2.0</td><td>19.2.0</td></tr></tbody></table>

**Notes.**

a IRAC channel, visit number, and wavelength.

b Astronomical observation requests.

c Total duration of the AOR in hours.

d Readout mode and frame time in seconds.

e Pipeline version of the basic calibrated data.

<table><caption>**Table 3**
Claret-4 Limb-darkening Coefficients for the WASP43 Star in the 3.6 and 4.5 μm Spitzer/IRAC Passbands</caption><thead><tr><th>Method</th><th>Channel</th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>a<sub>4</sub></th></tr></thead><tbody><tr><td rowspan="2">A17</td><td>Ch1, 3.6 μm</td><td>0.596193</td><td>-0.353618</td><td>0.234039</td><td>-0.070725</td></tr><tr><td>Ch2, 4.5 μm</td><td>0.574190</td><td>-0.585735</td><td>0.550643</td><td>-0.199545</td></tr><tr><td rowspan="2">A100</td><td>Ch1, 3.6 μm</td><td>0.575555</td><td>-0.288784</td><td>0.154642</td><td>-0.037594</td></tr><tr><td>Ch2, 4.5 μm</td><td>0.538245</td><td>-0.472568</td><td>0.411592</td><td>-0.141342</td></tr><tr><td rowspan="2">P100</td><td>Ch1, 3.6 μm</td><td>4.843472</td><td>-10.282954</td><td>10.828015</td><td>-4.187545</td></tr><tr><td>Ch2, 4.5 μm</td><td>4.846365</td><td>-10.290074</td><td>10.835657</td><td>-4.190506</td></tr><tr><td rowspan="2">PQS</td><td>Ch1, 3.6 μm</td><td>0.763637</td><td>0.265362</td><td>-0.458262</td><td>0.090083</td></tr><tr><td>Ch2, 4.5 μm</td><td>0.763610</td><td>0.265746</td><td>-0.458614</td><td>0.090143</td></tr></tbody></table>

**Note.** Using ATLAS9 and PHOENIX stellar-atmosphere models, calculated with the code by Espinoza & Jordán (2015), available at http://www.github.com/nespinoza/limb-darkening/.

technique (Morello 2015), obtaining similar or less robust results which we report in Appendix B, together with a detailed comparison of both techniques. Both algorithms rely on ICA, i.e., a blind source separation technique, to extract the instrumental components from the light curves. Such blind approaches have proven to perform as well as or better than other state-of-the-art pipelines to detrend *Spitzer*/IRAC observations of exoplanetary transits and eclipses (Morello et al. 2015; Ingalls et al. 2016).

In this work, the pixel-ICA pipelines were applied to full phase curve observations, which may be affected by detector systematics with longer time scales compared to the transit-only and eclipse-only observations. We checked for residual long-trend systematics by adding a linear or quadratic function of time in the light curve fits, and by comparing the differences in the Bayesian information criterion (BIC; Schwarz 1978) obtained with these various ramp models (constant, linear, or quadratic), as suggested by S17. Then, following the Occam's Razor principle, we confirmed the solution obtained with the constant ramp, if it had the lowest BIC. In an alternative case, the model selection was based on a number of considerations that will be explained in the following sections.

## 4. Results

The BIC favored the pure wavelet ICA + phase curve (constant ramp) models for the 4.5 μm and first 3.6 μm visits. For the second 3.6 μm visit the lowest BIC was obtained with the quadratic ramp model, while the BIC obtained with the constant ramp model was the highest (see Table 6). We observed that the best-fit astrophysical parameters do not significantly depend on the choice of the ramp parameterization, except for the phase curve parameters of the second 3.6 μm visit (see Section 4.1 and Appendix C).

Figures 2 and 3 show the raw light curves, the relevant best-fit models (with the minimum BIC), and the residuals. The rms amplitudes of the normalized residuals are $1.56 \times 10^{-3}$ for the first 3.6 μm visit, $1.52 \times 10^{-3}$ for the second 3.6 μm visit, and $1.87 \times 10^{-3}$ for the 4.5 μm visit. We estimate them to be ~24%, 22%, and 4% above the photon noise limit. Figure 4 shows how the rms amplitudes of the fitting residuals scale as a function of the bin size. The 4.5 μm residuals show no significant deviations from the theoretical behavior of white noise, different from the 3.6 μm residuals. The amount of residual correlated noise in the second 3.6 μm visit is notably smaller than in the first visit.
<figure><img src="image_3.png" /><figcaption>Stellar limb-darkening profiles</figcaption></figure>

**Figure 1.** Left panel: stellar limb-darkening profiles of the WASP43 star in the 3.6 µm Spitzer/IRAC channel, computed using the code provided by Espinoza & Jordán (2015) at http://www.github.com/nespinoza/limb-darkening/, with different settings: A17 (orange), A100 (red), P100 (blue), and PQS (cyan). Right panel: analogous plot for the 4.5 µm Spitzer/IRAC channel.

## 4.1. Phase Curve Models and Parameters

Figure 5 shows the best-fit phase curve models. Figure 6 reports the corresponding estimates of the planet dayside maximum and nightside minimum flux normalized to the stellar flux ($F_{\text{day}}^{\text{MAX}}$ and $F_{\text{night}}^{\text{MIN}}$), and their offsets relative to the mid-eclipse and mid-transit times ($\Delta\Phi_{\text{day}}^{\text{MAX}}$ and $\Delta\Phi_{\text{night}}^{\text{MIN}}$), respectively.

The results obtained for the second 3.6 µm data set with the quadratic ramp parameterization appear to be unphysical, yielding negative nightside fluxes, but still consistent with zero at the $1\sigma$ level. The results obtained with the linear ramp parameterization are more plausible, because they can be explained by a simpler physical model (see Section 5.2). We present here the two sets of results for the second 3.6 µm data set, together with the selected results for the other data sets. A more detailed discussion about the model selection criteria is reported in Section 5.1.

The (normalized) planet dayside flux at 4.5 µm is $(3.90 \pm 0.12) \times 10^{-3}$. The nightside flux is $(3.0 \pm 1.5) \times 10^{-4}$. The maximum dayside flux occurs $37 \pm 7$ minutes prior to the mid-eclipse time, which corresponds to a shift of $11^\circ 3 \pm 2^\circ 1$ east of the substellar point. The minimum nightside flux occurs within the interval $-10 \pm 15$ minutes relative to the mid-transit time, i.e., between $7^\circ 2$ east and $1^\circ 5$ west of the anti-stellar point.

The 3.6 µm phase curve models have remarkably different amplitudes and shapes, but with similar dayside fluxes: $(3.43 \pm 0.11) \times 10^{-3}$ for the first visit, and $(3.34 \pm 0.10) \times 10^{-3}$ (quadratic ramp) or $(3.32 \pm 0.10) \times 10^{-3}$ (linear ramp) for the second visit. The three estimates are mutually consistent within $0.5\sigma$ and smaller than the 4.5 µm maximum with $2-2.5\sigma$ significance level. The planet flux minima are $(6.9 \pm 1.6) \times 10^{-4}$ for the first 3.6 µm visit, and $(-1.6 \pm 1.9) \times 10^{-4}$ (quadratic ramp) or $(3.0 \pm 1.5) \times 10^{-4}$ (linear ramp) for the second visit.

For the second 3.6 µm visit, the phase curve maximum occurs $14 \pm 7$ minutes (quadratic ramp), or $18 \pm 9$ minutes

(linear ramp), earlier than the mid-eclipse time. These offsets correspond to hotspot shifts of $4^\circ 4 \pm 2^\circ 3$ and $5^\circ 6 \pm 2^\circ 7$ east of the substellar point. The phase curve minimum occurs at $+6 \pm 14$ minutes (quadratic ramp), or $0^{+16}_{-14}$ minutes (linear ramp), relative to the mid-transit time. These offsets correspond to shifts of $2^\circ \pm 4^\circ$ and $0^{+5}_{-4}^\circ$ west of the anti-stellar point.

The first 3.6 µm phase curve model is strongly asymmetric, with peaks occurring $64 \pm 13$ (maximum) and $103 \pm 18$ (minimum) minutes earlier than the mid-eclipse and mid-transit time, or, equivalently, $20^\circ \pm 4^\circ$ and $32^\circ \pm 6^\circ$ east of the substellar and anti-stellar points.

However, the tests reported in the Appendix C suggest that the true uncertainties in the peak offsets estimated for the 3.6 µm observations may be larger than the nominal error bars. For example, our estimate of the dayside peak offset for the first 3.6 µm visit becomes $18 \pm 9$ minutes before mid-eclipse when considering only the first two out of three AORs, which is identical to the estimate from the second visit (linear ramp). The corresponding dayside and nightside fluxes are consistent with those obtained from the full data set analysis within $1\sigma$. The same tests confirm the robustness of the parameter estimates for the 4.5 µm observation within their nominal error bars.

Table 4 reports our final measurements of the day and nightside fluxes and peak offsets at 3.6 and 4.5 µm. The results at 3.6 µm are the weighted averages between those obtained for the two visits, with inflated error bars for those parameters which were not consistent within $1\sigma$. We discarded the (unphysical) results obtained for the second 3.6 µm visit with a quadratic ramp, for reasons that will be further elaborated in Sections 5.1 and 5.2.

## 4.2. Transit Parameters

Figure 7 reports the best-fit transit depth ($p^2$), impact parameter ($b$), and transit duration ($T_0$) obtained with the four sets of limb-darkening coefficients reported in Table 3 (see Section 3.2). There appear to be systematic offsets between the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top panels: raw light curves (blue dots) obtained for the *Spitzer*/IRAC observations at 3.6 µm and relevant best-fit models (red line). Bottom panels: residuals from the above light curves and models (blue points), and standard deviations (black lines).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top panel: raw light curve (blue dots) obtained for the *Spitzer*/IRAC observations at 4.5 µm, and relevant best-fit model (red line). Bottom panel: residuals from the above light curve and model (blue points), and standard deviations (black lines).</figcaption></figure>

parameters obtained using the ATLAS and PHOENIX sets of coefficients. In particular, the PHOENIX models lead to smaller transit depths by ~400–700 parts per million (ppm),

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Normalized rms of residuals as function of bin size for the first 3.6 µm visit (green), second 3.6 µm visit using a quadratic (blue) or linear (dodger blue) ramp model, and 4.5 µm visit (red). The black dashed line shows the theoretical behavior for Gaussian residuals.</figcaption></figure>

smaller impact parameters by ~0.02–0.04, and longer transit durations by 55–70 s, at 3.6 and 4.5 µm, respectively. These differences correspond to two to five times the respective parameter error bars.

We did not find strong evidence in favor of one specific model (see Appendix D). Table 5 reports the arithmetic and weighted mean values of the geometric parameters, *b* and *T*₀, across the three visits for each limb-darkening model, and the mean transit depths at 3.6 and 4.5 µm. Table 5 also reports the global mean values over all the different limb-darkening models. While the absolute transit depths are model-dependent,
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left panel: best-fit phase curve models for the first 3.6 µm visit (green), second 3.6 µm visit using a quadratic (blue) or linear (dodger blue) ramp model, and 4.5 µm visit (red). The black horizontal line indicates the stellar flux level. Right panel: zoom-in of the left panel.</figcaption></figure>

the difference between the 3.6 and 4.5 µm transit depths is always consistent with zero within $1\sigma$.

# 5. Discussion

## 5.1. Reliability of the Model Selection Criteria

The minimum BIC solution for the second 3.6 µm visit (quadratic ramp) includes negative nightside flux values, which are unphysical. Even if the minimum nightside flux is consistent with being positive within $1\sigma$, the low upper limit poses a challenge for the modeling of exoplanetary atmospheres (e.g., Kataria et al. 2015). The solution obtained using a linear ramp parameterization, instead of quadratic, appears to be less problematic, as it is discussed in Section 5.2.

We tested model selection tools other than the BIC, which all agreed on the choice of the quadratic ramp model, although with different strengths of evidence. In particular, the $\Delta$BIC = 8.9 between the linear and quadratic parameterizations (see Table 6) denotes a strong, but not conclusive, preference for the latter according to Raftery (1995). The Akaike information criterion (AIC; Akaike 1974) and the consistent Akaike information criterion (CAIC; Bozdogan 1987) favor the quadratic ramp model more/less strongly than BIC ($\Delta$AIC = 15.6, $\Delta$CAIC = 7.9), because of a smaller/larger penalty for the number of model parameters. The deviance information criterion (Spiegelhalter et al. 2002) gives results similar to the AIC, while the Bayesian evidences calculated with MultiNest (Buchner et al. 2014) are consistent with the BIC estimates. Therefore, any information criterion weighted average of the alternative models, e.g., the marginalization method proposed by De Wit et al. (2016), would be driven by the quadratic ramp. A more sophisticated approach consists of marginalization over the hyperparameters of a Gaussian process (GP; Gibson 2014; Evans et al. 2015). In this paper, we did not pursue the GP method, because of its high computational cost and unclear performances in previous analyses of the *Spitzer*/IRAC data (Ingalls et al. 2016).

All of the tests discussed above rely on the relative amplitudes of the residuals. However, it is not guaranteed that smaller residuals correspond to more reliable parameter estimates. The potential errors in the instrumental systematics model may be compensated in part by biasing the retrieved astrophysical parameters, especially if the two sets of parameters (instrumental systematics and astrophysical) are correlated. We observed that, in the second 3.6 µm data set, the minimum nightside flux is strongly correlated with the two quadratic ramp coefficients, as measured by the absolute value of the Pearson correlation coefficients (PCCs ~0.6). The correlation with the linear ramp coefficient is much smaller (PCC ~ 0.1).

In conclusion, simple statistical criteria based on the amplitude of the best-fit residuals can provide useful guidelines to model selection, but they should not be considered alone. Physical plausibility may pose important constraints to the model selection, especially when the competing models have similar scores (Ingalls et al. 2016). We highly recommend to perform some self-consistency tests on the data, e.g., checking that the best-fit parameters do not vary dramatically if analyzing smaller portions rather than the whole data set (see Appendix C).

## 5.2. Exoplanet Disk-integrated Brightness Temperatures

For the *Spitzer*/IRAC 3.6 and 4.5 µm channels, the exoplanet flux contribution is predominantly thermal emission. By neglecting the nonthermal contributions, the observed planet-to-star flux ratio in a given passband is (Charbonneau et al. 2005)

$$ \frac{F_{p,\lambda}}{F_{*,\lambda}} = \left( \frac{R_p}{R_*} \right)^2 \frac{B_{\lambda}(T_p)}{B_{\lambda}(T_*)}, \qquad (2) $$

where $(R_p/R_*)^2$ is the sky-projected planet-to-star area ratio, $T_p$ is the phase-dependent exoplanet brightness temperature, and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left, top panel: maximum exoplanetary flux relative to the stellar flux, for the first 3.6 µm visit (green square), second 3.6 µm visit using a quadratic (blue) or linear (dodger blue) ramp model, weighted average between the first and second visit with a linear ramp (olive), and 4.5 µm visit (red square). Left, bottom panel: orbital phase of the maximum exoplanetary flux relative to mid-eclipse, including the values reported by S17 (black triangles). Right panels: analogous plots for the minimum exoplanetary to the stellar flux.</figcaption></figure>

<table><caption>**Table 4**
Phase Curve Parameters of WASP43 b</caption><thead><tr><th>Parameter</th><th>3.6 µm</th><th>4.5 µm</th></tr></thead><tbody><tr><td>F<sup>MAX</sup><sub>day</sub></td><td>(3.37 ± 0.07) × 10<sup>-3</sup></td><td>(3.90 ± 0.12) × 10<sup>-3</sup></td></tr><tr><td>F<sup>MIN</sup><sub>night</sub></td><td>(4.8 ± 1.4) × 10<sup>-4</sup></td><td>(3.0 ± 1.5) × 10<sup>-4</sup></td></tr><tr><td>ΔΦ<sup>MAX</sup><sub>day</sub></td><td>-0.028 ± 0.013</td><td>-0.031 ± 0.006</td></tr><tr><td>ΔΦ<sup>MIN</sup><sub>night</sub></td><td>-0.032 ± 0.030</td><td>-0.009 ± 0.013</td></tr></tbody></table>

**Note.** The 3.6 µm parameters are the weighted mean values over the two visits (using the linear ramp for the second visit); the error bars on $F_{\text{night}}^{\text{MIN}}$, $\Delta\Phi_{\text{day}}^{\text{MAX}}$, and $\Delta\Phi_{\text{night}}^{\text{MIN}}$ are inflated by the difference between the individual estimates in units of $\sigma$ (factors of 1.26, 2.15, and 3.05, respectively).

$T_*$ is the star brightness temperature. We computed the brightness temperatures for a star with $T_{\text{eff}} = 4400$ K, $\log g = 4.65$ (see Table 1), in the two IRAC passbands, by interpolating on a grid of PHOENIX stellar-atmosphere models (Husser et al. 2013). By inverting Equation (2), we calculated the exoplanet maximum dayside temperature, $T_{\text{day}}^{\text{MAX}}$, to be 1660 ± 21 K for the first 3.6 µm visit, 1643 ± 19 K (quadratic ramp) or 1639 ± 19 K (linear ramp) for the second 3.6 µm visit, and 1502 ± 18 K at 4.5 µm. The corresponding minimum nightside temperatures, $T_{\text{night}}^{\text{MIN}}$, are $1016_{-63}^{+56}$ K, $<568$ K, $837_{-103}^{+79}$ K, and $700_{-93}^{+68}$ K. These temperature estimates can be visualized in Figure 8.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top panel: transit depth estimates obtained with different sets of limb-darkening coefficients: A100 (red squares), A17 (orange circles), P100 (blue, upward triangles), and PQS (cyan, downward triangles). Middle and bottom panels: analogous plots for the impact parameter and for the transit duration.</figcaption></figure>

The choice of a linear or a quadratic ramp parameterization for the second 3.6 µm data set changes dramatically the inferred astrophysical scenario. The former leads to consistent measurements between the two 3.6 µm visits within 1.5σ, and lower brightness temperatures at 4.5 µm. The lower brightness temperatures suggest higher absorption/scatter at 4.5 µm within the WASP43 b atmosphere, assuming a non-inverted thermal profile (Blecic et al. 2014).

The quadratic ramp model implies lower nightside flux, and brightness temperature (upper limit), for the second 3.6 µm visit, suggesting some variability with the 2.5σ significance level. The wavelength-temperature trend is inverted between the exoplanet dayside and nightside, which indicates different atmospheric opacities between the two sides. These results might be explained with the appearance of high-altitude clouds in the exoplanet atmosphere during the second 3.6 µm visit. In conclusion, we cannot rule out this solution as physically impossible, but it is most likely biased by the strong parameter correlations, as mentioned in Section 5.1. This idea is reinforced by the simpler physical interpretation (and smaller parameter correlations) associated with the alternative solution using a linear ramp model.

### 5.3. Atmospheric Circulation Models

We used the 2D-ATMO code (Tremblin et al. 2017) to compute a grid of phase curve models for WASP43 b. The 2D-ATMO is an extension of 1D-ATMO (Tremblin et al. 2015) that takes into account the circulation induced by the irradiation from the host star at the equator of the planet. The atmospheric model for WASP43 b has been computed as part of a model comparison performed for future JWST observations (O. Venot et al. 2019, in preparation). The magnitude of the zonal wind is imposed at the substellar point at 4 km s⁻¹ and is computed accordingly to the momentum conservation law in the rest of the equatorial plane. The vertical mass flux is assumed to be proportional to the meridional mass flux with a proportionality constant 1/α; the wind is therefore purely longitudinal and meridional if α → ∞ or purely longitudinal and vertical for α → 0. As in Tremblin et al. (2017), a relatively low value of α drives the vertical advection of entropy/potential temperature in the deep atmosphere that can produce a hot interior, which can explain the inflated radii of hot Jupiters. A high value of α will produce a cold deep interior as in the standard 1D models. In this study, we used a simulation with α = 10⁴ that should be representative of WASP43 b since the planet is not highly inflated. In order to reproduce the low fluxes on the nightside of the planet, we added a simple cloud model consisting of a gray absorbing cloud deck with an absorption of 2.5 m² kg⁻¹ with a fixed bottom pressure of 0.1 bar. We explored different metallicities (1×, 3×, and 10× solar) and different top pressure levels for the cloud deck (0.1 bar, i.e., no clouds, 0.02, 0.01, 10⁻³, 10⁻⁴, and 10⁻⁷ bar).

Figure 9 compares the measured day and nightside fluxes and peak offsets (from Table 4) with those predicted by the atmospheric models. Figure 10 compares the whole phase curve profiles obtained from the data with the best matching profiles from the atmospheric models.

A number of atmospheric phase curve models are in excellent agreement with the observed profile at 4.5 µm. The best matches are the models with 3× or 1× solar metallicity and cloud top pressure of 10⁻³ bar. In both cases, the discrepancies between the fitted and the model profiles are smaller than 200 ppm with rms amplitudes of ~100 ppm. The corresponding model phase curve parameters are all consistent with the measured values within 1σ. In general, all of the models with cloud top pressure lower than 10⁻² bar are in good agreement with the 4.5 µm observation, but the models with 10× solar metallicity tend to predict a smaller dayside peak offset. The models with no clouds can be ruled out at the 4–8σ level in the nightside flux, and they also tend to predict significantly larger peak offsets, depending on the metallicity.

The results at 3.6 µm are more problematic, as the measured dayside flux is higher than predicted by the atmospheric models. The best match to the observed profile (average of the two observations) is the model with 10× solar metallicity and cloud top pressure of 2 × 10⁻² bar. In this case, the discrepancies between the fitted and the model profiles are within ~400 ppm with rms amplitudes below 300 ppm, which are larger than the error bars of the best-fit profile at certain orbital phases. The corresponding model phase curve parameters are consistent with the measured values within 1.5σ. The models with 3× and 1× solar metallicity predict smaller than measured dayside fluxes by ~150 ppm (2σ) and 300 ppm (4σ), respectively.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left panel: maximum dayside temperatures obtained in this work for the first 3.6 µm visit (green), second 3.6 µm visit using a quadratic (blue) or linear (dodger blue) ramp model, weighted average between the first and second visit with a linear ramp (olive), 4.5 µm visit (red), and temperatures reported by S17 (black triangles). Right panel: analogous plot for the minimum nightside temperatures. Note that S17 only reported 2σ upper limits for the nightside temperatures.</figcaption></figure>

<table><caption>**Table 5**
Mean Values of the Transit Parameters</caption><thead><tr><th>L-D</th><th>Mean</th><th>b</th><th>T<sub>0</sub> (s)</th><th>p<sup>2</sup> × 10<sup>-2</sup> (3.6 µm)</th><th>p<sup>2</sup> × 10<sup>-2</sup> (4.5 µm)</th></tr></thead><tbody><tr><td rowspan="2">A17</td><td>a.</td><td>0.655 ± 0.013</td><td>3480 ± 14</td><td>2.501 ± 0.019</td><td>2.504 ± 0.019</td></tr><tr><td>w.</td><td>0.657 ± 0.007</td><td>3479 ± 8</td><td>2.502 ± 0.013</td><td>2.504 ± 0.019</td></tr><tr><td rowspan="2">A100</td><td>a.</td><td>0.655 ± 0.013</td><td>3479 ± 14</td><td>2.501 ± 0.019</td><td>2.504 ± 0.019</td></tr><tr><td>w.</td><td>0.657 ± 0.007</td><td>3479 ± 8</td><td>2.503 ± 0.013</td><td>2.504 ± 0.019</td></tr><tr><td rowspan="2">P100</td><td>a.</td><td>0.624 ± 0.015</td><td>3543 ± 14</td><td>2.456 ± 0.019</td><td>2.434 ± 0.022</td></tr><tr><td>w.</td><td>0.629 ± 0.008</td><td>3543 ± 8</td><td>2.456 ± 0.013</td><td>2.434 ± 0.022</td></tr><tr><td rowspan="2">PQS</td><td>a.</td><td>0.630 ± 0.014</td><td>3537 ± 14</td><td>2.460 ± 0.018</td><td>2.443 ± 0.021</td></tr><tr><td>w.</td><td>0.635 ± 0.008</td><td>3538 ± 8</td><td>2.460 ± 0.013</td><td>2.443 ± 0.021</td></tr><tr><td>all</td><td>a.</td><td>0.641 ± 0.020</td><td>3510 ± 32</td><td>2.480 ± 0.026</td><td>2.471 ± 0.033</td></tr></tbody></table>

**Note.** The uncertainties in the overall arithmetic means are the standard deviations of the individual parameter values.

<table><caption>**Table 6**
ΔBIC and ΔAIC for the Different Long-term Ramp Models and Types of Observation</caption><thead><tr><th rowspan="2">Obs. Type</th><th rowspan="2">Ramp Model</th><th colspan="2">Ch2</th><th colspan="2">Ch1, Visit 2</th><th colspan="2">Ch1, Visit 1</th></tr><tr><th>ΔBIC</th><th>ΔAIC</th><th>ΔBIC</th><th>ΔAIC</th><th>ΔBIC</th><th>ΔAIC</th></tr></thead><tbody><tr><td rowspan="3">full</td><td>Constant</td><td><strong>0.0</strong></td><td>0.0</td><td>+310.8</td><td>+324.0</td><td><strong>0.0</strong></td><td>+10.7</td></tr><tr><td>Linear</td><td>+7.0</td><td>+0.4</td><td>+8.9</td><td>+15.6</td><td>+1.7</td><td>+5.7</td></tr><tr><td>Quadratic</td><td>+13.2</td><td>0.0</td><td><strong>0.0</strong></td><td>0.0</td><td>+2.6</td><td>0.0</td></tr><tr><td rowspan="3">ecl1 + tr</td><td>Constant</td><td><strong>0.0</strong></td><td>0.0</td><td>+2.2</td><td>+8.3</td><td><strong>0.0</strong></td><td>+6.1</td></tr><tr><td>Linear</td><td>+6.4</td><td>+0.1</td><td><strong>0.0</strong></td><td>0.0</td><td>+0.1</td><td>0.0</td></tr><tr><td>Quadratic</td><td>+13.0</td><td>+0.5</td><td>+7.3</td><td>+1.2</td><td>+8.2</td><td>+1.8</td></tr><tr><td rowspan="3">tr + ecl2</td><td>Constant</td><td><strong>0.0</strong></td><td>0.0</td><td><strong>0.0</strong></td><td>4.1</td><td><strong>0.0</strong></td><td>+9.6</td></tr><tr><td>Linear</td><td>+7.8</td><td>+1.6</td><td>+2.1</td><td>0.0</td><td>+5.8</td><td>+9.2</td></tr><tr><td>Quadratic</td><td>+16.0</td><td>+3.5</td><td>+9.0</td><td>+0.7</td><td>+2.8</td><td>0.0</td></tr></tbody></table>

**Note.** “full” = full phase curve; “ecl1 + tr” = half phase curve including the eclipse prior transit; “tr + ecl2” = half phase curve including the eclipse after transit. The boldface values correspond to the minimum BIC configurations.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Measured phase curve parameters as reported in Table 4, where the best parameter values are represented by continuous horizontal lines and the limits of the $1\sigma$ interval are represented by dashed horizontal lines. Phase curve parameters obtained from the atmospheric models with $1\times$ (black squares), $3\times$ (dark gray left-pointing triangles), and $10\times$ (light gray downward-pointing triangles) solar metallicity as a function of the cloud top pressure ($10^{-1}$ bar corresponds to the cloud-free models). The highlighted points correspond to the best matching phase curve models represented in Figure 10.</figcaption></figure>

The observations at 3.6 and 4.5 µm are best described by atmospheric models with different metallicity and cloud top pressure, although a range of models with metallicity higher than solar and cloud top pressure of ~$10^{-2}$ bar can reproduce all of the measured phase curve parameters within less than $2\sigma$. It is likely that a chemical composition different than scaled solar abundances could provide a better match to the data, without the need to find a compromise at the edges of the acceptable parameter ranges for the observations at the two wavelengths. In this paper, we refrain from speculating about the possible nature of the non-standard chemistry in the atmosphere of WASP43 b, which cannot be probed with the current data.

As the observations were not taken simultaneously, we cannot exclude some variability of the nightside clouds over the different visits. In Section 4.1 we noted that the parameter results derived from the individual 3.6 µm visits were not fully consistent at the $1\sigma$ level, but the apparent discrepancies might be caused by correlated noise in the fitting residuals. In general,

the observations with 4.5 µm channel are much less affected by correlated noise (Krick et al. 2016), which is also confirmed by the analyses in this paper (see Section 4 and Appendix C). Therefore, new observations at 4.5 µm would help to assess the level of variability in the atmosphere of WASP43 b.

### 5.4. Comparison with Previous Analyses of the Same Observations

The rms amplitudes of the light-curve fitting residuals obtained with our wavelet pixel-ICA are within 1% of those reported by M18, using an extension of the BLISS mapping algorithm (Stevenson et al. 2012).

Figure 11 reports the phase curve peak-to-peak amplitudes that we computed to compare with those reported in the previous literature. When taking the linear ramp model for the second 3.6 µm visit, our best-fit amplitudes are consistent with those reported by M18 within $1\sigma$, though our central values are
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left, top panel: observed 3.6 µm phase curve profile (continuous line and points with error bars), i.e., average of the best-fit profiles for the two visits (corresponding to the green and dodger blue curves in Figure 5), the best match from our grid of atmospheric phase curve models (dashed line), and the best match with solar metallicity (dotted line). The larger error bars for $\Phi < 0.5$ take into account the discrepancy between the best-fit profiles for the two visits. Left, bottom panel: residuals between the observed and two model phase curves. Right panels: analogous plots for the 4.5 µm phase curve.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Peak-to-peak phase curve amplitudes obtained in this work (squares), and reported by S17 (circles) and M18 (triangles). Same choice of colors as in Figures 6 and 8.</figcaption></figure>

~300 ppm larger. S17 obtained larger peak-to-peak amplitudes for the second 3.6 and 4.5 µm visits. When taking the quadratic ramp model for the second 3.6 µm visit, we also obtain a larger peak-to-peak amplitude, in agreement with S17. Overall, the different estimates for each visit are consistent at the 2σ level.

We observed that, in the analyses discussed here, a smaller peak-to-peak amplitude corresponds to a higher nightside flux (at a given wavelength). Figure 8 compares the brightness temperatures obtained in this work with those reported by S17. The dayside temperatures are consistent within 0.5σ. Our estimates of the nightside temperatures are higher than the 2σ upper limits reported by S17. Figure 6 shows that we obtained a significantly smaller dayside peak offset compared to S17 at 4.5 µm.

S17 could not find adequate approximations to the observed phase curve profiles with the cloud-free atmospheric models of Kataria et al. (2015). M18 computed new global circulation models with THOR (Mendonça et al. 2016); their best match to the Spitzer/IRAC data is a model with a nightside cloud deck with a top pressure of $10^{-2}$ bar and enhanced carbon dioxide (CO$_2$) with a longitudinal gradient. A visual inspection of the Figure 6 in M18 reveals that the maximum discrepancies between the fitted and the best model phase curves are ~400–500 ppm, i.e., equal or larger than those obtained in this study (see Section 5.3). Furthermore, the non-equilibrium CO$_2$ was introduced ad hoc by M18 to reproduce the low nightside flux observed at 4.5 µm, with a lower limit for the cloud top pressure of $10^{-2}$ bar (from Kreidberg et al. 2014). Mendonça et al. (2018b) could not physically explain that chemical disequilibrium.

Figure 12 reports our final estimates of the geometric parameters averaged over the three observations, and the analogous parameters derived from those reported by Stevenson et al. (2014) using HST/WFC3 data. Our estimates of the impact parameter using ATLAS limb-darkening coefficients are marginally consistent within 1σ with the value reported by Stevenson et al. (2014). We found longer transit duration than Stevenson et al. (2014) by 42 s and 101–106 s using ATLAS and PHOENIX limb-darkening coefficients, respectively. These discrepancies are above 3σ. Note that the small uncertainties in the parameters obtained from Stevenson et al. (2014) do not account for the degeneracy with the stellar limb-darkening.

Figure 13 compares the 3.6 and 4.5 µm transit depths with those obtained by S17. Unsurprisingly, our estimates using ATLAS (claret-4) limb-darkening coefficients best match the
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Left panel: impact parameter estimates obtained in this work (Table 5) and in the previous literature. This work: weighted mean values over the three observations using A100 (red square), A17 (orange circle), P100 (blue, upward triangle), and PQS (cyan, downward triangle) limb-darkening coefficients, and overall arithmetic mean (green diamond). The S14 result (black, leftward triangle) has been calculated, in this work, from the other parameters reported by Stevenson et al. (2014). Right panel: analogous plot for the transit duration.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** 3.6 and 4.5 µm transit depths estimates obtained in this work (Table 5) and by S17. Same choice of colors and symbols as in Figure 12.</figcaption></figure>

results obtained by S17, which also adopted ATLAS (quad-ratic) limb-darkening coefficients.

## 5.5. Comparison with Other Observations and with Other Exoplanets

Some authors suggested the existence of a simple relationship between the irradiation temperature and circulation efficiency of the exoplanetary atmospheres (Cowan & Agol 2011; Perez-Becker & Showman 2013; Schwartz & Cowan 2015; Komacek & Showman 2016; Keating & Cowan 2017). The S17 claim of zero circulation efficiency for the WASP43 b atmosphere injected an apparent outlier to the expected trend. WASP43 b and HD

209458 b have similar irradiation temperatures of ~2000–2100 K. Schwartz et al. (2017) reported $\epsilon = 0.49_{-0.14}^{+0.15}$ for HD 209458 b, based on visible-to-infrared observations, often limited to the secondary eclipses.

In this work, we obtained significantly higher nightside temperatures than the previous estimates by S17 for WASP43 b in the Spitzer/IRAC passbands. Assuming a blackbody-like emission, the circulation efficiency goes up to $\epsilon \sim 0.1–0.3$. We emphasize that this is just a broad estimate of the circulation efficiency in the WASP43 b atmosphere. In fact, the blackbody assumption is not valid, as revealed by the $>4.5\sigma$ difference between the 3.6 and 4.5 µm dayside temperatures (see Section 5.2). We propose a more direct comparison between the 4.5 µm phase
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Top panel: transit depth estimates obtained from the full phase curve, and from two transit-only analyses with different baselines (see Section 6) using the wavelet pixel-ICA (darker colors) and time pixel-ICA techniques (lighter colors). The horizontal dashed lines correspond to the full phase curve results. Middle and bottom panels: analogous plots for the impact parameter and transit duration.</figcaption></figure>

curves of WASP43 b and HD 209458 b. Zellem et al. (2014) reported $T_{day} = 1499 \pm 15$ K and $T_{night} = 972 \pm 44$ K for HD 209458 b at 4.5 µm. We obtained the same dayside temperature (within $0.2\sigma$) and ~200–300 K lower nightside temperature for WASP43 b at the same wavelength. These comparisons suggest the WASP43 b atmosphere may have a lower circulation efficiency than HD 209458 b, but the difference appears to be significantly smaller than from the original estimates reported in the literature. Furthermore, there are some hints of variability in the nightside cloud deck of WASP43 b (see Section 5.3), that would affect the temperature measurements. A new set of observations is desirable to test this hypothesis.

## 6. Transit-only Analyses

We analyzed smaller portions of the data sets as transit-only observations, in order to evaluate the pros and cons of the different observation types for exoplanet characterization. We tested two phase intervals, $|\Phi| \le 0.1$ and $|\Phi| \le 0.07$, which correspond to ~3.3 and 2.3 times the full transit duration. The fitted transit model does not include any phase curve

modulation or exoplanet nightside pollution, as is common practice in transit-only observations.

Figure 14 compares the transit parameters obtained from the phase curve and transit-only analyses. We found that the parameter error bars scale approximately as $\sqrt{N_{tot}/N_{out}}$, where $N_{tot}$ is the total number of data points, and $N_{out}$ is the number of out-of-transit points. The mathematical derivation of this result is reported in Appendix A (Equation (7)). More specifically, the error bars in transit depth obtained for the first 3.6 and 4.5 µm visits with the transit-only analyses are 9% and 8% (longer configuration), and 17% and 22% (shorter configuration) larger than those obtained from the respective full phase curves, which are within the ranges predicted by Equation (7) (8%–16% and 13%–29%, see Appendix A). For the second 3.6 µm observation, the transit-only error bars are larger than expected, i.e., 46% and 53% larger than the respective full phase curve error bars. It is worth noting that the second 3.6 µm observation is the only one that required a non-constant ramp model in the phase curve analysis. The error bars (in transit depth) obtained with the shorter transit-only configuration are 6%–12% larger
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Eclipse depth estimates from the eclipse-only analyses using the wavelet pixel-ICA (darker colors) and time pixel-ICA techniques (lighter colors). The horizontal dashed lines act as upper limits, i.e., the flux maxima obtained from the full phase curve analyses. The coordinates 1 and 2 refer to the first and second eclipse of each visit.</figcaption></figure>

than those obtained with the longer transit-only configuration for all of the observations, in good agreement with the range predicted by Equation (7) (4%–10%, see Appendix A). The other transit parameters have similar differences between the error bars obtained with the various configurations.

Interestingly, the transit depth estimated from the (wavelet) transit-only analyses, and especially from those with the longer phase interval, are slightly larger than those obtained from the full phase curve analyses in all cases. Although the three transit depth estimates for the same observation are mutually consistent within $1\sigma$, we found that the observed systematic behavior can be caused by the flat baseline approximation. In fact, the phase-dependent exoplanetary flux is higher before and after transit than during the transit, therefore increasing the apparent transit depth if this effect is not taken into account. Based on the best-fit phase curve model, the differences between the time-averaged exoplanetary flux out-of-transit and in-transit are in the range of 40–120 ppm. The difference in transit depths are of the same order of magnitude, but not identical because of the correlations with the other free parameters in the fit, which are also slightly biased. Consequently, the largest differences in transit depth are obtained for the 4.5 µm visit, i.e., +271 and +185 ppm for the longer and shorter transit-only analyses, respectively. The parameter offsets decrease significantly if the best-fit phase curve parameters are fixed in the transit-only analyses. While these potential bias are not statistically significant with the current error bars, they might become significant with the smaller error bars that are expected to be achieved with the next-generation instruments, such as those on board JWST and ARIEL. The potential bias are expected to be smaller for exoplanets with larger orbits, both because of the smaller day-night temperature contrast (phase curve amplitude) and the longer orbital period relative to the transit duration.

## 7. Eclipse-only Analyses

We also analyzed a portion of the data sets as eclipse-only observations. We fixed the orbital parameters to the estimates

obtained from the corresponding transit-only observation, then fitted for the eclipse depth and timing. There are no analogous eclipse depth estimates for the full phase curve analyses, as the out-of-eclipse flux is not constant, but the dayside maxima should represent upper limits for the eclipse depths. Figure 15 reports the eclipse depth estimates for the phase interval $|\Phi - 0.5| \le 0.1$, and the dayside maxima from the full phase curve analyses. The eclipse depths obtained with wavelet pixel-ICA at the same wavelength are mutually consistent within $1.5\sigma$, and they are all below the respective phase curve upper limits.

## 8. Conclusions

We analyzed three Spitzer/IRAC phase curves of the exoplanet WASP43 b at 3.6 µm (two observations) and 4.5 µm using a blind signal-source separation method, i.e., the wavelet pixel-ICA. We assessed the robustness of the results by analyzing both the full and the half phase curves with different instrument ramp models. We revealed a significant degeneracy between stellar limb-darkening and transit parameters, and a potential bias obtained analyzing only the transit portion of the phase curve. This bias is of the order of ~100 ppm in transit depth for WASP43 b in the mid-infrared, and it is expected to be highly significant for the upcoming JWST and ARIEL observations. We found higher nightside temperatures, smaller hotspot offsets, and greater consistency (~1σ) between the two 3.6 µm visits than those reported by Stevenson et al. (2017). Our results point toward a greater circulation efficiency of the WASP43 b atmosphere, in agreement with an empirical trend between irradiation temperature and circulation efficiency. Additionally, we compared the observed phase curves with a grid of atmospheric models, enabling quantitative estimates of the cloud top pressure. Our phase curve parameter results are consistent within $1\sigma$ with those reported in a recent reanalysis by Mendonça et al. (2018a), but we provide an alternative interpretation with a lower cloud top pressure instead of invoking a strong
disequilibrium chemistry. Furthermore, we proposed a simple formula for estimating how the error bars scale with the duration of the observations. Such formula can be used for optimizing the trade-off between parameters precision and duration of the observations.

This work was supported by the LabEx P2IO and the French ANR contract 05-BLAN-NT09-573739. The research leading to these results has received funding from the European Unions Horizon 2020 Research and Innovation Programme, under grant Agreement n 776403. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

## Appendix A
Scaling Relation for the Error Bars in Transit Depth

We derive here a simple analytical formula to estimate the error bar relative to the transit depth as a function of the time spent observing the out-of-transit. We consider a simplified case, with flat out-of-transit and in-transit, no stellar limb-darkening, and neglecting the transit ingress and egress. In this case, the transit depth is equal to

$$p^2 = \frac{F_{\text{out}} - F_{\text{in}}}{F_{\text{out}}} = 1 - \frac{F_{\text{in}}}{F_{\text{out}}}, \quad (3)$$

where $F_{\text{in}}$ and $F_{\text{out}}$ are the constant flux values in-transit and out-of-transit, respectively. If the data are only affected by Gaussian noise, the flux values can be estimated with the following error bars:

$$\Delta F_{\text{in}} = \frac{\sigma}{\sqrt{N_{\text{in}}}}, \quad \Delta F_{\text{out}} = \frac{\sigma}{\sqrt{N_{\text{out}}}}, \quad (4)$$

where $\sigma$ is the standard deviation of the Gaussian noise, and $N_{\text{in}}$ and $N_{\text{out}}$ are the numbers of in-transit and out-of-transit data points, respectively.

We calculate the error bar in transit depth, $p^2$, using the law of propagation of error (Taylor 1996),

$$\begin{aligned} \Delta p^2 &= \sqrt{\left(\frac{\partial p^2}{\partial F_{\text{in}}} \Delta F_{\text{in}}\right)^2 + \left(\frac{\partial p^2}{\partial F_{\text{out}}} \Delta F_{\text{out}}\right)^2} = \\ &= \sqrt{\left(-\frac{\Delta F_{\text{in}}}{F_{\text{out}}}\right)^2 + \left(\frac{F_{\text{in}}}{F_{\text{out}}^2} \Delta F_{\text{out}}\right)^2}. \end{aligned} \quad (5)$$

By injecting Equation (4) into Equation (5), we obtain

$$\Delta p^2 = \frac{\sigma}{F_{\text{out}}} \sqrt{\frac{1}{N_{\text{in}}} + \frac{F_{\text{in}}^2}{F_{\text{out}}^2 N_{\text{out}}}}. \quad (6)$$

Now, we make the approximation $F_{\text{in}} \approx F_{\text{out}} = F$, obtaining

$$\Delta p^2 \approx \frac{\sigma}{F} \sqrt{\frac{1}{N_{\text{in}}} + \frac{1}{N_{\text{out}}}} = \frac{\sigma}{F} \sqrt{\frac{N_{\text{tot}}}{N_{\text{in}} N_{\text{out}}}}. \quad (7)$$

We estimate that, for typical transit depth values up to ~3%, the impact of this approximation is less than 0.1% in $\Delta p^2$.

The formula in Equation (7) should provide a lower limit for the error bars. In a more realistic case, the error bars will be larger because of the non-flatness introduced by the stellar

limb-darkening and the phase curve modulations, and, in general, because of correlations between the free parameters in the fit. In this work, we found that the error bars in transit depth are ~20%–50% larger than those estimated using Equation (7) with $N_{\text{in}}$ ranging from the number of data points between the second to third contact points and the number of data points between the first to fourth contact points (Seager & Mallén-Ornelas 2003).

Equation (7) provides a useful tool to predict how the error bars can scale with the longer observations, then to optimize the trade-off between observing time and precision with the future missions.

## Appendix B
Time versus Wavelet Pixel-ICA

The core of the pixel-ICA method is the ICA transform of a set of pixel time series into maximally independent components, i.e., a linear transformation that minimizes the mutual information (Hyvärinen & Oja 2000). In the wavelet pixel-ICA algorithm the pixel time series undergo discrete wavelet transform (DWT) before the ICA separation, and the independent components are transformed back into the time domain. More specifically, we adopt one-level DWT with mother wavelet Daubechies-4 (Daubechies 1992).

One of the independent components has the morphology of the astrophysical signal (transit, eclipse, or phase curve), the other components represent the instrumental systematics. We model fit the sum-of-pixel time series, the so-called raw light curve, as a linear combination of a parametric model of the astrophysical signal (instead of the relevant independent component) and the other independent components.

The MCMC error bars are then rescaled as

$$\sigma_{\text{par}} = \sigma_{\text{par},0} \sqrt{\frac{\sigma_0^2 + \sigma_{\text{ICA}}^2}{\sigma_0^2}}, \quad (8)$$

where $\sigma_0^2$ is the likelihood variance, approximately equal to the variance of the residuals, and $\sigma_{\text{ICA}}^2$ is a term accounting for the uncertainty in the ICA components. The latter term is calculated as (Morello et al. 2016)

$$\sigma_{\text{ICA}}^2 = \sum_j o_j^2 \text{ISR}_j, \quad (9)$$

where ISR is the so-called interference-to-signal-ratio matrix (Tichavský et al. 2008), and $o_j$ are the best-fit coefficients of the linear combination.

Figure 16 compares the best-fit phase curve model obtained using pixel-ICA in the time or wavelet domains. The two approaches led to similar phase curve models for the 4.5 µm observation. Instead, the models obtained for the 3.6 µm visits using time pixel-ICA are less reliable, as they assume strong negative emission from the exoplanet nightside. Figure 17 shows that the worse phase curve models are associated with higher levels of correlated noise in the fitting residuals, though, in some cases, the rms amplitudes are smaller. This study suggests that the alternative pixel-ICA algorithms are equivalent below a certain level of correlated noise (e.g., at 4.5 µm), otherwise the wavelet-based approach outperforms the analysis in the time domain.

For the transit-only analyses, the parameters obtained with the time pixel-ICA are consistent with those obtained with the
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Best-fit phase curve models obtained using the wavelet pixel-ICA (darker colors) and time pixel-ICA techniques (lighter colors). The right panels are a zoom-in of the left panels.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Left panels: normalized rms of residuals as a function of bin size obtained using the wavelet pixel-ICA (darker colors) and time pixel-ICA techniques (lighter colors). The black dashed lines show the theoretical behavior for Gaussian residuals. Right panels: ratio between the rms of residuals obtained using the wavelet pixel-ICA and time pixel-ICA techniques as a function of bin size. The black dashed lines denotes the separation (ratio = 1); the points below the lines correspond to the case of smaller residuals obtained using the wavelet pixel-ICA (and vice versa).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Top, left panel: maximum exoplanetary flux, relative to the stellar flux, for the first 3.6 µm visit from the full and half phase curve analyses using the different ramp models (see Appendix C). The letters “A” and “B” indicate the minimum AIC and BIC solutions among the different ramp models. Top, right panel: analogous plot for the minimum exoplanetary flux. Bottom, left panel: analogous plot for the offset of the phase curve maximum relative to mid-eclipse. Bottom, right panel: analogous plot for the offset of the phase curve minimum relative to mid-transit.</figcaption></figure>

wavelet pixel-ICA within less than 0.5σ (see Figure 14). It is reasonable to expect that the impact of low-frequency noise is smaller over the transit timescale, therefore explaining the apparent equivalence of the two methods.

For the eclipse-only analyses, the eclipse depths obtained with the time pixel-ICA are systematically larger, and with larger error bars than those obtained with the wavelet pixel-ICA. In some cases, the eclipse depth estimates obtained with the time pixel-ICA are above their phase curve upper limits. Therefore, the wavelet pixel-ICA outperforms the time pixel-ICA in eclipse-only analyses. This fact was already observed in the previous literature (Morello et al. 2016), and it is attributed to the smaller signal-to-noise ratio of the eclipse signal. For the 4.5 µm eclipses (least correlated noise), the two methods lead to consistent results within 1σ.

## Appendix C: Half Phase Curves

We discuss here the analysis of the so-called half phase curves, i.e., continuous observations including one transit and either the eclipse immediately before or after. In this study, we obtain two half phase curves from each visit by considering two out of three consecutive AORs (where applicable). For the second 3.6 µm visit, we split the first AOR into two parts in order to get three pieces as in the other visits. Note that the half phase curves within the same visit share the same transit event.

Figures 18–20 report the phase curve parameters for the full and half phase curve analyses with a constant, linear, and quadratic ramp model. Figures 21–23 report the corresponding transit parameters. Table 6 reports the relative ΔBIC and ΔAIC for the different ramp models.
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Top, left panel: maximum exoplanetary flux, relative to the stellar flux, for the second 3.6 µm visit from the full and half phase curve analyses using the different ramp models (see Appendix C). The letters “A” and “B” indicate the minimum AIC and BIC solutions among the different ramp models. Top, right panel: analogous plot for the minimum exoplanetary flux. Bottom, left panel: analogous plot for the offset of the phase curve maximum relative to mid-eclipse. Bottom, right panel: analogous plot for the offset of the phase curve minimum relative to mid-transit.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Top, left panel: maximum exoplanetary flux, relative to the stellar flux, for the 4.5 µm visit from the full and half phase curve analyses using the different ramp models (see Appendix C). The letters “A” and “B” indicate the minimum AIC and BIC solutions among the different ramp models. Top, right panel: analogous plot for the minimum exoplanetary flux. Bottom, left panel: analogous plot for the offset of the phase curve maximum relative to mid-eclipse. Bottom, right panel: analogous plot for the offset of the phase curve minimum relative to mid-transit.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 21.** Top panel: transit depth estimates for the first 3.6 µm visit from the full and half phase curve analyses using the different ramp models (see Appendix C). The letters “A” and “B” indicate the minimum AIC and BIC solutions among the different ramp models. Middle and bottom panels: analogous plots for the impact parameter and for the transit duration.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Top panel: transit depth estimates for the second 3.6 µm visit from the full and half phase curve analyses using the different ramp models (see Appendix C). The letters “A” and “B” indicate the minimum AIC and BIC solutions among the different ramp models. Middle and bottom panels: analogous plots for the impact parameter and for the transit duration.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 23.** Top panel: transit depth estimates for the 4.5 µm visit from the full and half phase curve analyses using the different ramp models (see Appendix C). The letters “A” and “B” indicate the minimum AIC and BIC solutions among the different ramp models. Middle and bottom panels: analogous plots for the impact parameter and for the transit duration.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 24.** Normalized rms of residuals as a function of the bin size for the first 3.6 µm visit. The full phase curve analysis is represented as the dark green line and the half phase curve analysis including the eclipse prior transit as the light green line. The black dashed line shows the theoretical behavior for Gaussian residuals.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 25.** Relative chi-square obtained using different sets of limb-darkening coefficients: A100 (red squares), A17 (orange circles), P100 (blue, upward triangles), and PQS (cyan, downward triangles). The full markers refer to the full phase curve residuals. The empty markers refer to the subinterval of the same residuals centered on the transit. The horizontal lines delimit the significance levels according to Raftery (1995): $\Delta\chi^2 \le 2$ is not significant, $2 < \Delta\chi^2 \le 6$ denotes positive evidence against the model with higher $\chi^2$, and $6 < \Delta\chi^2 \le 10$ denotes strong evidence against the model with higher $\chi^2$.</figcaption></figure>

The phase curve parameters are more degenerate with the ramp parameters in the half phase curve models, as suggested by the much larger and asymmetric error bars. The transit parameters are independent on the choice of ramp model, half or full phase curve, i.e., their dispersion are much smaller than their error bars. Only for the first 3.6 µm visit, the analysis of the half phase curve, including the eclipse before the transit, outperforms the full phase curve analysis. Figure 24 shows that the correlated noise in the half phase curve residuals is significantly smaller than in the full phase curve residuals. Also, the dayside shift resulting from the half phase curve is in better agreement with the results from the other observations.

## Appendix D
Limb-darkening Coefficients

Figure 25 shows the $\chi^2$ differences between the light curve fits with the different sets of limb-darkening coefficients. In all cases, the $\chi^2$ differences are smaller than 2, except the 2.7 difference between the P100 and A100 models of the 4.5 µm light curve. Such differences are not significant or barely significant according to Raftery (1995). Given that the limb-darkening coefficients only affect the points during the transit, we recalculated the $\chi^2$ differences over the phase interval $|\Phi| \le 0.1$. The new differences only indicate a strong evidence ($\Delta\chi^2 = 7.1$) against the P100 model, and a positive evidence ($\Delta\chi^2 = 4.5$) against the PQS model at 4.5 µm.

Figure 26 compares the two transit models obtained with the P100 and A100 limb-darkening coefficients at 4.5 µm, that led to the largest $\Delta\chi^2$ in the light curve residuals. Figure 27 shows the difference between the corresponding light curve residuals. The difference is non-zero only during the transit and the two eclipses<fnref n="5" /> and the maximum peaks are ~600 ppm. The rms amplitude of the residuals is 1870 ppm, i.e., more than three times larger than the maximum difference.

<footnote n="5">⁵ Even if the stellar limb-darkening does not affect the eclipse shape, the eclipse duration is affected, as it is equal to the transit duration.</footnote>
<figure><img src="image_28.png" /><figcaption>**Figure 26.** Left panel: best-fit transit models for the 4.5 µm visit obtained using A100 (red) and P100 (blue) limb-darkening coefficients. Right panel: difference between the alternative transit models.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 27.** Top panel: light curve residuals of the 4.5 µm visit obtained using A100 (red) and P100 (blue) limb-darkening coefficients. Bottom panel: difference between the residual time series above. Note that the only differences occur during the transit and eclipses. The difference between the residuals obtained with different limb-darkening coefficients is smaller for the other observations.</figcaption></figure>

## ORCID iDs

G. Morello https://orcid.org/0000-0002-4262-5661

C. Daniëlski https://orcid.org/0000-0002-3729-2663

P. Tremblin https://orcid.org/0000-0001-6172-3403

## References

* Akaike, H. 1974, ITAC, 19, 716

* Bean, J. L., Stevenson, K. B., Batalha, N. M., et al. 2018, PASP, 130, 114402

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

* Bozdogan, H. 1987, Psychometrika, 52, 345

* Buchner, J., Georgakakis, A., Nadra, K., et al. 2014, A&A, 564, A125

* Charbonneau, D., Allen, L. E., Thomas Megeath, S., et al. 2005, ApJ, 626, 523

* Claret, A. 2000, A&A, 363, 1081

* Claret, A., Hauschildt, P. H., & Witte, S. 2012, A&A, 546, A14

* Claret, A., Hauschildt, P. H., & Witte, S. 2013, A&A, 552, A16

* Cooper, C. S., & Showman, A. P. 2005, ApJL, 629, L45

* Cowan, N. B., & Agol, E. 2011, ApJ, 729, 54

* Cowan, N. B., Agol, E., & Charbonneau, D. 2007, MNRAS, 379, 641

* Daubechies, I. 1992, Ten Lectures on Wavelets (Philadelphia, PA: SIAM)

* De Wit, J., Wakeford, H. R., Gillon, M., et al. 2016, Natur, 537, 69

* Espinoza, N., & Jordán, A. 2015, MNRAS, 450, 1879

* Evans, T. M., Aigrain, S., Gibson, N., et al. 2015, MNRAS, 451, 680

* Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

* Fortney, J. J., Cooper, C. S., Showman, A. P., Marley, M. S., & Freedman, R. S. 2006, ApJ, 652, 746

* Gibson, N. P. 2014, MNRAS, 445, 3401

* Hellier, C., Anderson, D. R., Collier Cameron, A., et al. 2011, A&A, 535, L7

* Howarth, I. D. 2011, MNRAS, 418, 1165

* Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

* Hyvärinen, A., & Oja, E. 2000, NN, 13, 411

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

* IRAC Instrument & Instrument Support Teams 2015, IRAC Instrument Handbook, v2.1.2, http://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/iracinstrumenthandbook/

* Kataria, T., Showman, A. P., Fortney, J. J., et al. 2015, ApJ, 801, 86

* Keating, D., & Cowan, N. B. 2017, ApJL, 849, L5

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

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

* Krick, J. E., Ingalls, J., Carey, S., et al. 2016, ApJ, 824, 27

* Kurucz, R. L. 1979, ApJS, 40, 1

* Mandel, K., & Agol, E. 2002, ApJL, 580, L171

* Mendonça, J. M., Grimm, S. L., Grosheintz, L., & Heng, K. 2016, ApJ, 829, 115

* Mendonça, J. M., Malik, M., Demory, B.-O., & Heng, K. 2018a, AJ, 155, 150

* Mendonça, J. M., Tsai, S.-M., Malik, M., Grimm, S. L., & Heng, K. 2018b, ApJ, 869, 107

* Morello, G. 2015, ApJ, 808, 56

* Morello, G., Tsiaras, A., Howarth, I. D., & Homeier, D. 2017, AJ, 154, 111

* Morello, G., Waldmann, I. P., Tinetti, G., et al. 2015, ApJ, 802, 117

* Morello, G., Waldmann, I. P., & Tinetti, G. 2016, ApJ, 820, 86

* Neilson, H. R., & Lester, J. B. 2013a, A&A, 554, A98

* Neilson, H. R., & Lester, J. B. 2013b, A&A, 556, A86

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

* Raftery, A. E. 1995, Sociological Methodology, 25, 111

* Schwartz, J. C., & Cowan, N. B. 2015, MNRAS, 449, 4192

* Schwartz, J. C., Kashner, Z., Jovmir, D., & Cowan, N. B. 2017, ApJ, 850, 154

* Schwarz, G. 1978, AnSta, 6, 461

* Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038

* Showman, A. P., & Guillot, T. 2002, A&A, 385, 166
Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, A. 2002, *Journal of the Royal Statistical Society Series B*, **64**, 583

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

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

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

Taylor, J. R. 1996, An Introduction to Error Analysis (2nd ed.; Mill Valley, CA: Univ. Science Books)

Tichavský, P., Koldovský, Z., Yeredor, A., Gómez-Herrero, G., & Doron, E. 2008, *ITNN*, **19**, 421

Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, *ApJL*, **804**, L17

Tremblin, P., Chabrier, G., Mayne, N. J., et al. 2017, *ApJ*, **841**, 30

Wu, X., Roby, T., & Ly, L. 2010, *Proc. SPIE*, **7737**, 773716

Zellem, R. T., Lewis, N. K., Knutson, H. A., et al. 2014, *ApJ*, **790**, 53

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