An independent analysis of the Spitzer/IRAC phase curves of WASP43 b

TL;DR

This paper reanalyzes Spitzer phase curves of exoplanet WASP43 b using a novel blind signal separation technique, leading to updated planetary parameters and insights into atmospheric circulation.

Contribution

It introduces the wavelet pixel-Independent Component Analysis method for analyzing phase curves, improving parameter estimates and robustness over previous analyses.

Findings

Higher nightside temperatures detected

Smaller hotspot offsets observed

Results align better with atmospheric models

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 $\mu$m, and one visit at 4.5 $\mu$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 ($\sim$1 $\sigma$) between the two 3.6~$\mu$m visits, and (4) a greater similarity with the predictions of the atmospheric circulation models. Our parameter results are consistent within 1 $\sigma$ 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.

Tables

Table 1. Table 1: WASP43 system parameters

Stellar parameters
$T_{\text{eff}}$ (K)	4400$\pm$200
$\log g_{\ast }$ (cgs)	4.65${}_{}{}^{+0.06}{}_{-0.04}{}^{}$
$\left[\text{Fe/H}\right]$ (dex)	-0.05$\pm$0.17
$M_{\ast }$ ($M_{\odot }$)	0.58$\pm$0.05
$R_{\ast }$ ($R_{\odot }$)	0.60${}_{}{}^{+0.03}{}_{-0.04}{}^{}$
Planetary parameters
$M_p$ ($M_{\text{Jup}}$)	1.78$\pm$0.10
$R_p$ ($R_{\text{Jup}}$)	0.93${}_{}{}^{+0.07}{}_{-0.09}{}^{}$
$a$ (au)	0.0142$\pm$0.0004
Transit parameters
$p^2$	0.0255$\pm$0.0012
$b$	0.66${}_{}{}^{+0.04}{}_{-0.07}{}^{}$
$i$ (deg)	82.6${}_{}{}^{+1.3}{}_{-0.9}{}^{}$
$P$ (days)	0.813475$\pm$0.000001
$E.T.$ (HJD)	2455528.86774$\pm$0.00014

Table 2. Table 2: Spitzer /IRAC data sets analyzed for this study.

Obs.aaIRAC channel, visit number, and wavelength.	Prog. ID	AORsbbAstronomical Observation Requests.	UT Date	$\Delta$t (h)ccTotal duration of the AOR in hours.	ModeddReadout mode and frame time in seconds.	Pip. eePipeline version of the Basic Calibrated Data.
Ch1, visit 1	11001	52364544	2015 Mar 7	8.5	sub, 2.0	19.2.0
(3.6 $\mu$m)		52364800	2015 Mar 7	8.5	sub, 2.0	19.2.0
		52355312	2015 Mar 8	8.5	sub, 2.0	19.2.0
Ch1, visit 2	11001	57744384	2015 Sep 5	15.2	sub, 2.0	19.2.0
(3.6 $\mu$m)		57744640	2015 Sep 5	10.2	sub, 2.0	19.2.0
Ch2	10169	51777024	2014 Aug 27	8.5	sub, 2.0	19.2.0
(4.5 $\mu$m)		51777280	2014 Aug 27	8.5	sub, 2.0	19.2.0
		51777792	2014 Aug 28	8.5	sub, 2.0	19.2.0

Table 3. Table 3: Claret-4 limb-darkening coefficients for the WASP43 star in the 3.6 and 4.5  μ 𝜇 \mu m Spitzer /IRAC passbands.

Method	Channel	$a_1$	$a_2$	$a_3$	$a_4$
A17	Ch1, 3.6 $\mu$m	0.596193	-0.353618	0.234039	-0.070725
	Ch2, 4.5 $\mu$m	0.574190	-0.585735	0.550643	-0.199545
A100	Ch1, 3.6 $\mu$m	0.575555	-0.288784	0.154642	-0.037594
	Ch2, 4.5 $\mu$m	0.538245	-0.472568	0.411592	-0.141342
P100	Ch1, 3.6 $\mu$m	4.843472	-10.282954	10.828015	-4.187545
	Ch2, 4.5 $\mu$m	4.846365	-10.290074	10.835657	-4.190506
PQS	Ch1, 3.6 $\mu$m	0.763637	0.265362	-0.458262	0.090083
	Ch2, 4.5 $\mu$m	0.763610	0.265746	-0.458614	0.090143

Table 4. Table 4: Phase curve parameters of WASP43 b.

Parameter	3.6 $\mu$m	4.5 $\mu$m
$F_{\text{day}}^{\text{MAX}}$	(3.37$\pm$0.07)$\times$10-3	(3.90$\pm$0.12)$\times$10-3
$F_{\text{night}}^{\text{MIN}}$	(4.8$\pm$1.4)$\times$10-4	(3.0$\pm$1.5)$\times$10-4
$\Delta {\Phi }_{\text{day}}^{\text{MAX}}$	-0.028$\pm$0.013	-0.031$\pm$0.006
$\Delta {\Phi }_{\text{night}}^{\text{MIN}}$	-0.032$\pm$0.030	-0.009$\pm$0.013

Table 5. Table 5: Mean values of the transit parameters.

L-D	Mean	$b$	$T_0$ (s)	$p^2\times$10-2 (3.6 $\mu$m)	$p^2\times$10-2 (4.5 $\mu$m)
A17	a.	0.655$\pm$0.013	3480$\pm$14	2.501$\pm$0.019	2.504$\pm$0.019
	w.	0.657$\pm$0.007	3479$\pm$8	2.502$\pm$0.013	2.504$\pm$0.019
A100	a.	0.655$\pm$0.013	3479$\pm$14	2.501$\pm$0.019	2.504$\pm$0.019
	w.	0.657$\pm$0.007	3479$\pm$8	2.503$\pm$0.013	2.504$\pm$0.019
P100	a.	0.624$\pm$0.015	3543$\pm$14	2.456$\pm$0.019	2.434$\pm$0.022
	w.	0.629$\pm$0.008	3543$\pm$8	2.456$\pm$0.013	2.434$\pm$0.022
PQS	a.	0.630$\pm$0.014	3537$\pm$14	2.460$\pm$0.018	2.443$\pm$0.021
	w.	0.635$\pm$0.008	3538$\pm$8	2.460$\pm$0.013	2.443$\pm$0.021
all	a.	0.641$\pm$0.020	3510$\pm$32	2.480$\pm$0.026	2.471$\pm$0.033

Table 6. Table 6: Δ Δ \Delta BIC and Δ Δ \Delta AIC for the different long-term ramp models and types of observation.

		Ch2		Ch1, Visit 2		Ch1, Visit 1
Obs. type	Ramp model	$\Delta$BIC	$\Delta$AIC	$\Delta$BIC	$\Delta$AIC	$\Delta$BIC	$\Delta$AIC
full	Constant	0.0	0.0	+310.8	+324.0	0.0	+10.7
	Linear	+7.0	+0.4	+8.9	+15.6	+1.7	+5.7
	Quadratic	+13.2	0.0	0.0	0.0	+2.6	0.0
ecl1 + tr	Constant	0.0	0.0	+2.2	+8.3	0.0	+6.1
	Linear	+6.4	+0.1	0.0	0.0	+0.1	0.0
	Quadratic	+13.0	+0.5	+7.3	+1.2	+8.2	+1.8
tr + ecl2	Constant	0.0	0.0	0.0	4.1	0.0	+9.6
	Linear	+7.8	+1.6	+2.1	0.0	+5.8	+9.2
	Quadratic	+16.0	+3.5	+9.0	+0.7	+2.8	0.0

Equations

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

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

$$p^2=\frac{F_{\text{mbox}out}-F_{\text{mbox}in}}{F_{\text{mbox}out}}=1-\frac{F_{\text{mbox}in}}{F_{\text{mbox}out}},$$

$$\Delta F_{\text{mbox}in}=\frac{\sigma }{\sqrt{N_{\text{mbox}in}}},\Delta F_{\text{mbox}out}=\frac{\sigma }{\sqrt{N_{\text{mbox}out}}},$$

$$\Delta p^2=\sqrt{{\left(\frac{\partial p^2}{\partial F_{\text{mbox}in}}\Delta F_{\text{mbox}in}\right)}^2+{\left(\frac{\partial p^2}{\partial F_{\text{mbox}out}}\Delta F_{\text{mbox}out}\right)}^2}=$$

$$=\sqrt{{\left(-\frac{\Delta F_{\text{mbox}in}}{F_{\text{mbox}out}}\right)}^2+{\left(\frac{F_{\text{mbox}in}}{F_{\text{mbox}out}^2}\Delta F_{\text{mbox}out}\right)}^2}.$$

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

$$\Delta p^2\approx \frac{\sigma }{F}\sqrt{\frac{1}{N_{\text{mbox}in}}+\frac{1}{N_{\text{mbox}out}}}=\frac{\sigma }{F}\sqrt{\frac{N_{\text{mbox}tot}}{N_{\text{mbox}in}{N}_{\text{mbox}out}}}.$$

$${\sigma }_{par}={\sigma }_{par,0}\sqrt{\frac{{\sigma }_0^2+{\sigma }_{ICA}^2}{{\sigma }_0^2}},$$

$${\sigma }_{ICA}^2=\sum _j{o}_j^{2}IS{R}_j,$$

Full text

planets and satellites: individual (WASP43 b) - planets and satellites: atmospheres - planets and satellites: fundamental parameters - stars: atmospheres - techniques: photometric - techniques: spectroscopic

1 INTRODUCTION

WASP43 b is a hot Jupiter orbiting around a K7 V star in $\sim$19.5 hours (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, hereinafter 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 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: $\varepsilon=$0.002${}_{-0.002}^{+0.01}$, where $\varepsilon$ 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 $\mu$m data set, which presented discrepant results, and larger correlated noise in the light curve residuals. S17 also discarded a $\sim$2 hr interval from the second 3.6 $\mu$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 by 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, $\varepsilon\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. (2018), hereinafter M18, reanalyzed the three Spitzer/IRAC phase curves without discarding any data. M18 found a better agreement between the two 3.6 $\mu$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, 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

[TABLE]

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-coefficient limb-darkening (Claret, 2000), hereinafter claret-4, for the WASP43 star in the 3.6 and 4.5 $\mu$m Spitzer/IRAC passbands, using the code provided by Espinoza & Jordán (2015) at GitHub111http://www.github.com/nespinoza/limb-darkening/. The code adopts two grids of stellar-atmosphere intensity models, i.e., ATLAS9222http://kurucz.harvard.edu/grids.html (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, 2013, 2013b; Espinoza & Jordán, 2015). We tested the following fitting options:

• •

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

• •

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

• •

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

• •

PQS, i.e., least-squares fit to the PHOENIX model intensities with $\mu\geq$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_{\mbox{\footnotesize 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\geq$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$\times$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 $\sim$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 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 $\mu$m and first 3.6 $\mu$m visits. For the second 3.6 $\mu$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 $\mu$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 $\mu$m visit, 1.52$\times$10*-3* for the second 3.6 $\mu$m visit, and 1.87$\times$10*-3* for the 4.5 $\mu$m visit. We estimate them to be $\sim$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 $\mu$m residuals show no significant deviations from the theoretical behavior of white noise, different from the 3.6 $\mu$m residuals. The amount of residual correlated noise in the second 3.6 $\mu$m visit is notably smaller than in the first visit.

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_{\mbox{\footnotesize day}}^{\mbox{\footnotesize MAX}}$ and $F_{\mbox{\footnotesize night}}^{\footnotesize MIN}$), and their offsets relative to the mid-eclipse and mid-transit times ($\Delta\Phi_{\mbox{\footnotesize day}}^{\mbox{\footnotesize MAX}}$ and $\Delta\Phi_{\mbox{\footnotesize night}}^{\footnotesize MIN}$) respectively.

The results obtained for the second 3.6 $\mu$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 $\mu$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 $\mu$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.3${}^{\circ}\pm$2.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.2∘ east and 1.5∘ west of the anti-stellar point.

The 3.6 $\mu$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 $\mu$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 $\mu$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 $\mu$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.4${}^{\circ}\pm$2.3∘ and 5.6${}^{\circ}\pm$2.7∘ east of the substellar point. The phase curve minimum occurs at +6$\pm$14 minutes (quadratic ramp), or 0${}_{-14}^{+16}$ minutes (linear ramp), relative to the mid-transit time. These offsets correspond to shifts of 2${}^{\circ}\pm$4∘ and 0${}_{-4}^{+5}$∘ west of the anti-stellar point.

The first 3.6 $\mu$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$\pm$4∘ and 32${}^{\circ}\pm$6∘ 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 $\mu$m observations may be larger than the nominal error bars. For example, our estimate of the dayside peak offset for the first 3.6 $\mu$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 $\mu$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 $\mu$m. The results at 3.6 $\mu$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 (unphyisical) results obtained for the second 3.6 $\mu$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 parameters obtained by using the ATLAS and PHOENIX sets of coefficients. In particular, the PHOENIX models lead to smaller transit depths by $\sim$400–700 parts per million (ppm), smaller impact parameters by $\sim$0.02–0.04, and longer transit durations by 55–70 s, at 3.6 and 4.5 $\mu$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_{0}$, across the three visits for each limb-darkening model, and the mean transit depths at 3.6 and 4.5 $\mu$m. Table 5 also reports the global mean values over all the different limb-darkening models. While the absolute transit depths are model-dependent, the difference between the 3.6 and 4.5 $\mu$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 $\mu$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 by 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 $\mu$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 $\sim$0.6). The correlation with the linear ramp coefficient is much smaller (PCC$\sim$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 $\mu$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)

[TABLE]

where $(R_{p}/R_{*})^{2}$ is the sky-projected planet-to-star area ratio, $T_{p}$ is the phase-dependent exoplanet brightness temperature, and $T_{*}$ is the star brightness temperature. We computed the brightness temperatures for a star with $T_{\mbox{\footnotesize 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_{\mbox{\footnotesize day}}^{\mbox{\footnotesize MAX}}$, to be 1660$\pm$21 K for the first 3.6 $\mu$m visit, 1643$\pm$19 K (quadratic ramp) or 1639$\pm$19 K (linear ramp) for the second 3.6 $\mu$m visit, and 1502$\pm$18 K at 4.5 $\mu$m. The corresponding minimum nightside temperatures, $T_{\mbox{\footnotesize night}}^{\mbox{\footnotesize MIN}}$, are 1016${}_{-\mbox{63}}^{+\mbox{56}}$ K, $<$568 K, 837${}_{-\mbox{103}}^{+\mbox{79}}$ K, and 700${}_{-\mbox{93}}^{+\mbox{68}}$ K. These temperature estimates can be visualized in Figure 8.

The choice of a linear or a quadratic ramp parameterization for the second 3.6 $\mu$m data set changes dramatically the inferred astrophysical scenario. The former leads to consistent measurements between the two 3.6 $\mu$m visits within 1.5 $\sigma$, and lower brightness temperatures at 4.5 $\mu$m. The lower brightness temperatures suggest higher absorption/scatter at 4.5 $\mu$m within the WASP43 b atmosphere, assuming 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 $\mu$m visit, suggesting some variability with 2.5 $\sigma$ 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 $\mu$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 (Venot et al., in preparation). The magnitude of the zonal wind is imposed at the substellar point at 4 km s*-1* 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/\alpha$; the wind is therefore purely longitudinal and meridional if $\alpha\rightarrow\infty$ or purely longitudinal and vertical for $\alpha\rightarrow 0$. As in Tremblin et al. (2017), a relatively low value of $\alpha$ 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 $\alpha$ will produce a cold deep interior as in the standard 1D models. In this study, we used a simulation with $\alpha$=104 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 m2 kg*-1* with a fixed bottom pressure of 0.1 bar. We explored different metallicities (1$\times$, 3$\times$, and 10$\times$ solar) and different top pressure levels for the cloud deck (0.1 bar, i.e., no clouds, 0.02, 0.01, 10*-3*, 10*-4* and 10*-7* 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 $\mu$m. The best matches are the models with 3$\times$ or 1$\times$ solar metallicity and cloud top pressure of 10*-3* bar. In both cases, the discrepancies between the fitted and the model profiles are smaller than 200 ppm with rms amplitudes of $\sim$100 ppm. The corresponding model phase curve parameters are all consistent with the measured values within 1 $\sigma$. In general, all of the models with cloud top pressure lower than 10*-2* bar are in good agreement with the 4.5 $\mu$m observation, but the models with 10$\times$ solar metallicity tend to predict a smaller dayside peak offset. The models with no clouds can be ruled out at the 4-8 $\sigma$ level in the nightside flux, and they also tend to predict significantly larger peak offsets, depending on the metallicity.

The results at 3.6 $\mu$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$\times$ solar metallicity and cloud top pressure of 2$\times$10*-2* bar. In this case, the discrepancies between the fitted and the model profiles are within $\sim$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 $\sigma$. The models with 3$\times$ and 1$\times$ solar metallicity predict smaller than measured dayside fluxes by $\sim$150 ppm (2 $\sigma$) and 300 ppm (4 $\sigma$), respectively.

The observations at 3.6 and 4.5 $\mu$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 $\sim$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 $\mu$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 $\mu$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 $\mu$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 $\mu$m visit, our best-fit amplitudes are consistent with those reported by M18 within 1 $\sigma$, though our central values are $\sim$300 ppm larger. S17 obtained larger peak-to-peak amplitudes for the second 3.6 $\mu$m and 4.5 $\mu$m visits. When taking the quadratic ramp model for the second 3.6 $\mu$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 $\sigma$ 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 $\sigma$. Our estimates of the nightside temperatures are higher than the 2 $\sigma$ upper limits reported by S17. Figure 6 shows that we obtained a significantly smaller dayside peak offset compared to S17 at 4.5 $\mu$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 top pressure of 10*-2* bar and enhanced carbon dioxide (CO2) 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 $\sim$400–500 ppm, i.e., equal or larger than those obtained in this study (see Section 5.3). Furthermore, the non-equilibrium CO2 was introduced ad hoc by M18 to reproduce the low nightside flux observed at 4.5 $\mu$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 by using ATLAS limb-darkening coefficients are marginally consistent within 1 $\sigma$ 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 by using ATLAS and PHOENIX limb-darkening coefficients, respectively. These discrepancies are above 3 $\sigma$. 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 $\mu$m transit depths with those obtained by S17. Unsurprisingly, our estimates using ATLAS (claret-4) limb-darkening coefficients best match the results obtained by S17, which also adopted ATLAS (quadratic) 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 HD209458 b have similar irradiation temperatures of $\sim$2000–2100 K. Schwartz et al. (2017) reported $\varepsilon=$0.49${}_{-0.14}^{0.15}$ for HD209458 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 $\varepsilon\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 $\mu$m dayside temperatures (see Section 5.2). We propose a more direct comparison between the 4.5 $\mu$m phase curves of WASP43 b and HD209458 b. Zellem et al. (2014) reported $T_{\mbox{\footnotesize day}}=$1499$\pm$15 K and $T_{\mbox{\footnotesize night}}=$972$\pm$44 K for HD209458 b at 4.5 $\mu$m. We obtained the same dayside temperature (within 0.2 $\sigma$) and $\sim$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 HD209458 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 pros and cons of the different observation types for exoplanet characterization. We tested two phase intervals, $|\Phi|\leq 0.1$ and $|\Phi|\leq 0.07$, which correspond to $\sim$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.

Figures 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_{\mbox{\footnotesize tot}}/N_{\mbox{\footnotesize out}}}$, where $N_{\mbox{\footnotesize tot}}$ is the total number of data points, and $N_{\mbox{\footnotesize out}}$ is the number of out-of-transit points. The mathematical derivation of this result is reported in Appendix A (Equation A5). More specifically, the error bars in transit depth obtained for the first 3.6 and 4.5 $\mu$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 A5 (8–16$\%$ and 13–29$\%$, see Appendix A). For the second 3.6 $\mu$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 $\mu$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 than those obtained with the longer transit-only configuration for all of the observations, in good agreement with the range predicted by Equation A5 (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 differences 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 $\mu$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 onboard 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 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|\leq 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 $\mu$m (two observations) and 4.5 $\mu$m using a blind signal-source separation method, i.e., the wavelet pixel-ICA. We assessed the robustness of the results by analysing 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 $\sim$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 ($\sim$1 $\sigma$) between the two 3.6 $\mu$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. (2018), 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 Union’s 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

[TABLE]

where $F_{\mbox{\footnotesize in}}$ and $F_{\mbox{\footnotesize 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:

[TABLE]

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

We calculate the error bar in transit depth, $p^{2}$, by using the “law of propagation of error” (Taylor, 1996):

[TABLE]

By injecting Equation A2 into Equation A, we obtain

[TABLE]

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

[TABLE]

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

The formula in Equation A5 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 $\sim$20–50$\%$ larger than those estimated using Equation A5 with $N_{\mbox{\footnotesize 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 A5 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 vs 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

[TABLE]

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

[TABLE]

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 $\mu$m observation. Instead, the models obtained for the 3.6 $\mu$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 $\mu$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 wavelet pixel-ICA within less than 0.5 $\sigma$ (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 $\mu$m eclipses (least correlated noise), the two methods lead to consistent results within 1 $\sigma$.

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 $\mu$m visit, we split the first AOR in 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 $\Delta$BIC and $\Delta$AIC for the different ramp models.

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 $\mu$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 $\mu$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|\leq$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 $\mu$m.

Figure 26 compares the two transit models obtained with the P100 and A100 limb-darkening coefficients at 4.5 $\mu$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 eclipses333Even if, the stellar limb-darkening does not affect the eclipse shape, the eclipse duration is affected, as it is equal to the transit duration. and the maximum peaks are $\sim$600 ppm. The rms amplitude of the residuals is 1870 ppm, i.e., more than three times larger than the maximum difference.