Repeatability of Spitzer/IRAC exoplanetary eclipses with Independent Component Analysis
Abstract
The research of effective and reliable detrending methods for Spitzer data is of paramount importance for the characterization of exoplanetary atmospheres. To date, the totality of exoplanetary observations in the mid and farinfrared, at wavelengths 3 m, have been taken with Spitzer. In some cases, in the past years, repeated observations and multiple reanalyses of the same datasets led to discrepant results, raising questions about the accuracy and reproducibility of such measurements. Morello et al. (2014, 2015) proposed a blindsource separation method based on the Independent Component Analysis of pixel time series (pixelICA) to analyze IRAC data, obtaining coherent results when applied to repeated transit observations previously debated in the literature. Here we introduce a variant to pixelICA through the use of wavelet transform, wavelet pixelICA, which extends its applicability to lowS/N cases. We describe the method and discuss the results obtained over twelve eclipses of the exoplanet XO3b observed during the “Warm Spitzer” era in the 4.5 m band. The final results are reported, in part, also in (Ingalls et al., 2016), together with results obtained with other detrending methods, and over ten synthetic eclipses that were analyzed for the “IRAC Data Challenge 2015”. Our results are consistent within 1 with the ones reported in Wong et al. (2014) and with most of the results reported in Ingalls et al. (2016), which appeared on the arXiv while this paper was under review. Based on many statistical tests discussed in Ingalls et al. (2016), the wavelet pixelICA method performs as well as or better than other stateofart methods recently developed by other teams to analyze Spitzer/IRAC data, and, in particular, it appears to be the most repeatable and the most reliable, while reaching the photon noise limit, at least for the particular dataset analyzed. Another strength of the ICA approach is its highest objectivity, as it does not use prior information about the instrument systematics, making it a promising method to analyze data from other observatories. The selfconsistency of individual measurements of eclipse depth and phase curve slope over a span of more than three years proves the stability of Warm Spitzer/IRAC photometry within the error bars, at the level of 1 part in 10 in stellar flux.
1 Introduction
Observations of transits and eclipses of exoplanets have been used, in the last decade, to characterize the nature of more than 100 of these alien worlds with unprecedented detail. Molecular, atomic and ionic signatures have been detected in the atmospheres of exoplanets through transmission spectroscopy, i.e. multiwavelength measurements of the transit depth, showing differential absorption/scatter of the stellar light through the planetary limb (e.g. Charbonneau et al., 2002; Tinetti et al., 2010; Deming et al., 2013). For the brightest targets, emission spectra have been measured during (planetary) eclipses to constrain the atmospheric chemistry, pressuretemperature profile, albedo, and global circulation patterns (e.g. Charbonneau et al., 2005; Swain et al., 2009). Many datasets were obtained using the InfraRed Array Camera (IRAC, Fazio et al., 2004) onboard Spitzer Space Telescope. Since depletion of the telescope’s cryogen in 2009, IRAC continues to operate in the 3.6 and 4.5 m bands, as the “Spitzer Warm Mission”.
A precision of 0.01 is required to study the atmospheres of giant exoplanets through transmission and/or emission spectroscopy (Brown, 2001; Tessenyi et al., 2013; Waldmann et al., 2015, 2015b). Detrending instrumental systematics in raw data is necessary to achieve the target precision. It is not always obvious how to decorrelate the data using auxiliary information of the instrument and, in some cases, different methods lead to significantly different results (see e.g. Waldmann, 2012; Morello et al., 2015). The majority of exoplanets’ transit and eclipse multiband photometric data adopted in the literature are obtained by combining data at different wavelengths from separate epochs years apart. Repeated observations are necessary to estimate the overall level of variability, due to astrophysical variations and possible uncorrected instrument effects. If consistent, compared to single observations, repeated observations can provide more accurate parameter values, leading to higher signaltonoise atmospheric spectra for exoplanets, and potentially allowing the characterization of smaller exoplanets around fainter stars.
The main instrumental effect affecting Spitzer/IRAC data at 3.6 and 4.5 m is due to intrapixel gain variations and spacecraftinduced motion, socalled pixelphase effect. The measured flux from the star is correlated with its position on the detector, hence the idea of correcting the data with a polynomial function of the stellar centroid as proposed in the literature (Charbonneau et al., 2005; MoralesCaldéron et al., 2006; Stevenson et al., 2010; Beaulieu et al., 2011). Multiple reanalyses of the same datasets with the polynomial method show that, in some cases, results can be sensitive to some specific options/variants, such as the degree of the polynomial adopted, partial data rejection, including temporal or other decorrelations (e.g. Beaulieu et al., 2008; Désert et al., 2009; Stevenson et al., 2010; Beaulieu et al., 2011; Knutson et al., 2011). Recently, several alternative methods have been proposed to decorrelate Spitzer/IRAC data: gain mapping (Ballard et al., 2010; Cowan et al., 2012; Knutson et al., 2012; Lewis et al., 2013; Zellem et al., 2014), bilinearly interpolated subpixel sensitivity mapping (BLISS, Stevenson et al., 2012), Independent Component Analysis using pixel time series (pixelICA, Morello et al., 2014, 2015), pixellevel decorrelation (PLD, Deming et al., 2015), and Gaussian Process models (Evans et al., 2015). A comparison between pixelICA and PLD methods is reported in Morello (2015). The discussion about the detrending methods, their performances, reliability and potential biases for Spitzer/IRAC data is a hot topic. In the “IRAC Data Challenge 2015” different methods have been tested over synthetic data created by the IRAC team, which contain ten simulated eclipse observations of the exoplanet XO3b (Ingalls et al., 2016). Researchers were also encouraged to reanalyze a similar set of real observations obtained in the 4.5 m band.
In this paper, we describe an evolution of the pixelICA method proposed in Morello et al. (2014, 2015); Morello (2015), and present the results obtained by applying said method to the analysis of twelve eclipses of the exoplanet XO3b taken with Warm Spitzer/IRAC in the 4.5 m band. PixelICA method differs from the other detrending methods proposed in the literature, as it is an unsupervised machine learning algorithm. The lack of any prior assumptions about the instrument systematics and astrophysical signals ensures a high degree of objectivity, and indicates that the same method could be applied to detrend data taken with different instruments. PixelICA method gave coherent results when applied over multiple transit observations of the exoplanets HD189733b and GJ436b, for which the previous literature reported discrepant results (Morello et al., 2014, 2015), and over simulated observations with a variety of instrumental systematics (Morello, 2015). Similar techniques have been used in the literature to detrend Spitzer/IRS and Hubble/NICMOS data, the main difference being in the choice of the input time series (Waldmann, 2012, 2014; Waldmann et al., 2013). The ability of ICA to decorrelate nongaussian signals is inherently limited to a low gaussian white noise amplitude relative to the nongaussian signals. In this paper, we propose a wavelet pixelICA algorithm, which outperforms the traditional pixelICA algorithm in lowS/N observations, extending applicability to planetary eclipses taken during the “Warm Spitzer” era.
XO3b is a hot Jupiter (11.7 0.5 , JohnsKrull et al., 2008; Hirano et al., 2011) in an eccentric orbit (0.283 0.003, Knutson et al., 2014) with a period of 3.19 days and orbital semimajor axis of AU (Winn et al., 2008). The host star is F5V with 6760 80 K, and 4.24 0.03 (Winn et al., 2008; Torres et al., 2012). A previous analysis of the 12 eclipses reported an average eclipse depth of 1.58 10 relative to the outofeclipse flux (star+planet), and a phase curve slope of 6.0 10 days (Wong et al., 2014). Here we compare our results with the ones reported in Wong et al. (2014).
2 Data Analysis
2.1 Observations
We analyze twelve eclipse observations of XO3b taken with Spitzer/IRAC in the 4.5 m band. Ten individual eclipses were observed over 6 months (Nov 11, 2012  May 24, 2013), including two sets of three consecutive eclipses, another eclipse is contained within a fullorbit observation on May 5, 2013 (PID: 90032). Each individual observation consists of 14,912 frames over 8.4 hr using IRAC’s subarray readout mode with 2.0 s integration time. In subarray mode 64 frames are taken consecutively, the reset time is 1 s. We extracted 14,912 frames from the fullorbit observation to analyze the lightcurve of the eclipse over a time interval similar to other observations. The last eclipse was extracted out of a 66 hr observation on April 8, 2010 (PID: 60058). Table 1 reports the dates in which the eclipses were observed.
Obs. Number  UT Date  Orbit Number 

1  2010 Apr 8  0 
2  2012 Nov 11  297 
3  2012 Nov 17  299 
4  2012 Nov 20  300 
5  2012 Nov 23  301 
6  2012 Dec 2  304 
7  2012 Dec 9  306 
8  2013 Apr 22  348 
9  2013 May 5  352 
10  2013 May 18  356 
11  2013 May 21  357 
12  2013 May 24  358 
2.2 The eclipse model
In our model the stellar flux is constant in time, and normalized to 1. We adopt the formalism of Mandel & Agol (2002) for the eclipse model, accounting for the planet being occulted by the star. We approximate the planet’s phase curve in the region of the eclipse as a linear function of the time, the slope is called “phase constant”, adopting the same terminology of Wong et al. (2014). The slope is only due to planet’s flux variations. While the planet is completely occulted by the star, the flux is constantly 1. The eclipse depth is defined as the (unseen) planet’s flux at the center of eclipse, in units of stellar flux (see Figure 1). When fitting the eclipse models, the orbital parameters are fixed to the values reported in Table 2, taken from Wong et al. (2014), while the center of eclipse, eclipse depth and phase constant are free parameters to determine. We also considered models with zero phase curve’s slope (see Appendix A.5).
orbital period, (days)  3.19153285 

scaled semimajor axis,  7.052 
inclination, (deg)  84.11 
eccentrity,  0.2833 
argument of periastron, (deg)  346.8 
2.3 Wavelet ICA
Continuous Wavelet Transform (CWT)
The wavelet transform (WT) decomposes a given signal, , into its frequency components. This is done by convolving the time signal with a basis of highly localized impulses or “wavelets”. To fix the ideas, we assume that is a time series, although this is not necessary, as the DWT can be applied to different kinds of signals. The individual wavelet functions are derived from a single “mother wavelet”, , through translation and dilation of the mother wavelet. The mathematical definition of the CWT is
(1) 
(2) 
where is the mother wavelet for a given scaling and translation , and is the wavelet coefficient with respect to and .
If the wavelet basis is orthogonal, the inverse wavelet transform can be used to reconstruct the original time series:
(3) 
The mother wavelet can be chosen among a variety of wavelet families with different properties. For more details we refer to the relevant literature, e.g. Daubechies (1992); Percival & Walden (2000).
Discrete Wavelet Transform (DWT)
Astronomical data are usually in the form of discrete time series. For the DWT the mother wavelet is denoted by , and the scaling function, also called “father wavelet”, is denoted by . The mother and father wavelets act as highpass and lowpass frequency filters, respectively. They are related by
(4) 
where is the filter length and corresponds to the number of points in the time series .
The onelevel DWT is defined by
(5) 
(6) 
The time series approximates the underlying lowfrequency trend of (average coefficients), while the time series represents a higher frequency component (detail coefficients). They are downsampled by a factor of 2 (“” in Equations 5 and 6) with respect to the original time series, because of the Nyquist theorem.
It is possible to apply the and filters to the time series, then obtaining new sets of coefficients, and , and iterate the process. The nlevel DWT includes the series of average coefficients, downsampled by a factor of 2, and n series  of detail coefficients, representing bands of higher frequencies.
The original data can be reconstructed by reversing the process:
(7) 
Wavelet ICA
In this section we describe the wavelet ICA algorithm, which is used in a variety of contexts, such as medical sciences (e.g. La Foresta et al., 2006; Inuso et al., 2007; Mammone et al., 2012), engineering (e.g. Lin & Zhang, 2005), acoustic (e.g. Moussaoui et al., 2006; Zhao et al., 2006), image denoising (e.g. Karande & Talbar, 2014), and astrophysics (e.g. Waldmann, 2014).
Be the column vector of observed signals. In this paper, are individual pixel time series, socalled pixellightcurves, which are mixtures of different source signals, astrophysical or instrumental in nature, and gaussian noise. The formalisms adopted in this subsection is valid in a more general context, where the can be any kind of mixed signals. ICA is a linear transformation of the observed (mixed) signals which minimizes the mutual information to decorrelate the independent components:
(8) 
where is the column vector of the original source signals, A is the matrix of mixing coefficients, and . We refer the reader to Waldmann (2012); Morello et al. (2014, 2015) for additional details.
The ability of ICA to decorrelate nongaussian signals is inherently limited to a low gaussian noise amplitude relative to the nongaussian signals. Wavelet ICA algorithms are designed to be less sensitive to white noise compared to the simple ICA separation, described above. In wavelet ICA algorithms, the DWT is applied to the observed signals:
(9) 
(10) 
The ICA separation is performed with the series of coefficients:
(11) 
The independent components series of coefficients are:
(12) 
They can be converted into the time domain through inverse DWT (Equation 7).
The DWT preliminarly separates the highfrequency components from the lowfrequency trend, enhancing the ability of ICA to disentangle the lowfrequency independent components. This step is particularly important in cases where the gaussian noise is dominant. Additional processing options/variants have been proposed in the literature to further improve the ICA performances in specific contexts, e.g. coefficients’ thresholding (Stein, 1981; Donoho, 1995), suppression of some frequency ranges (Lin & Zhang, 2005), taking individual levels as input to ICA (Inuso et al., 2007; Mammone et al., 2012). In this paper, we aim to provide the most objective analysis of the datasets, with minimal prior assumptions, hence those variants are not considered. The impact of those variants and other operations to the data will be carefully investigated in future studies.
2.4 Detrending method, lightcurve fitting and error bars
In this Section we list the main steps of the wavelet pixelICA method, followed by a more accurate description and comments:

Selecting an array of pixels. The raw lightcurve is the sum of the individual pixel time series within the selected array.

Removing outliers.

Subtracting the background from the raw lightcurve.

Computing the wavelet transforms of the time series from the pixels within the selected array, hereafter called pixellightcurves.

Performing ICA decomposition of the wavelettransformed pixellightcurves.

Computing the inverse wavelet transforms of the independent components.

Simultaneous fitting of the components (except the eclipse one) and astrophysical model on the raw lightcurve.

Estimating parameter error bars.
Selecting the pixelarray
We use squared arrays of pixels as photometric apertures; in this paper, we tested 55 and 77 arrays with the stellar centroid at their centers. By default, all pixellightcurves within the selected array are also used to decorrelate the instrument systematics through ICA. Our previous analyses of transit observations indicated the 55 and 77 arrays to be optimal choices, and, in general, results were very little affected by the choice of different arrays (Morello et al., 2014, 2015).
Outlier rejection
We flag and correct outliers in the flux time series. First, we calculate the standard deviations of any set of five consecutive points and take the median value as the representative standard deviation. We define the expected value in one point as the median of the four closest points, i.e. two before and two after. Points differing from their expected values by more than five times the standard deviation are flagged as outliers. They are then replaced with the mean value of the points immediately before and after, or, in case of two consecutive outliers, with a linear interpolation between the closest points which are not outliers. We checked that the outliers removed after this process are coincident with outliers that would have been spotted by eye. They are less than 0.35 the number of points in each observation.
Background subtraction
The background is estimated by taking the mean flux over four arrays of pixels with the same size of the selected array (55 or 77) near the corners of the subarray area. In Appendix A.3 we discuss how this preliminary step improves the results. Here we anticipate that the impact on individual measurements of the eclipse depth is at the level of 10, well below the error bars. However, this difference might become significant when averaging results from multiple observations, and the error bars are reduced.
Wavelet transforms
The main novelty of the algorithm proposed in this paper is that the pixellightcurves are wavelettransformed before performing the ICA separation. More specifically, we adopt onelevel DWT with mother wavelet Daubechies4 (Daubechies, 1992). We found that different choices of the mother wavelet, among different families and numbers, lead to exactly the same results, and higherlevel DWTs are not useful. We also investigated the effect of binning the time series (see Appendix A.1).
ICA decomposition
It is performed on the wavelettransformed pixellightcurves (see Equation 11).
Inverse wavelet transforms
The independent components are transformed in the time domain through inverse DWTs (see Equation 7).
MCMC fitting
The raw lightcurves are linear combinations of the independent components. One of the components is the eclipse signal (with some residual systematics), other components may be instrumental systematics and/or other astrophysical signals. Instead of fitting an eclipse model to the eclipse component, more robust estimates of the eclipse parameters are obtained by fitting a linear combination of the eclipse model and the noneclipse components, hereafter full models, to the raw lightcurves. The free parameters of the eclipse model are fitted together with the scaling coefficients of the independent components. First estimates of the parameters and scaling coefficients are obtained through a NelderMead optimization algorithm (Lagarias et al., 1998); they are then used as optimal starting points for an Adaptive Metropolis algorithm with delayed rejection (Haario et al., 2006), generating chains of 300,000 values. The output chains are samples of the posterior (gaussian) distributions. We adopt the mean values of the chains as final best estimates of the parameters, and the standard deviations as zeroorder error bars, .
Final error bars
The final parameter error bars are:
(13) 
is the sampled likelihood variance, approximately equal to the variance of the residuals for the best transit model; is a term accounting for the potential bias of the components obtained with ICA.
(14) 
where ISR is the socalled InterferencetoSignalRatio matrix, and are the coefficients of the noneclipse components. The term relative to the goodness of the fit of the components does not appear in Equation 14, as it is automatically included in when the components’ coefficients and astrophysical parameters are fitted simultaneously.
2.5 Results
Figure 2 shows the raw lightcurves for the twelve observations, the correspondent detrended eclipses and best models, obtained using the 55 array and binning over 32 frames, i.e. 64 s. The individual measurements of eclipse depth and phase constant are reported in Figure 3 and Table 3.
Obs. Number  Depth (10)  1 error  Phase constant (10 days)  1 error 
1  1.66  0.11  6  7 
2  1.72  0.11  6  8 
3  1.54  0.10  8  5 
4  1.56  0.10  9  7 
5  1.52  0.10  9  7 
6  1.56  0.13  2  10 
7  1.64  0.11  5  10 
8  1.57  0.12  0  11 
9  1.54  0.11  0  5 
10  1.52  0.10  11  8 
11  1.50  0.12  16  7 
12  1.48  0.12  8  6 
The results from all epochs are consistent within the error bars, suggesting the lack of any detectable astrophysical variability for this system, and residual instrument variability. By taking the weighted means of the individual measurements, we obtain global best estimates of (1.57 0.03)10 for the eclipse depth, and (4.4 2.0)10 days for the phase constant.
3 Discussion
3.1 Reduced chisquared tests
The underlying assumption for the weighted mean to be a valid parameter estimate is that individual measurements of that parameter are normally distributed around the same mean value with variances , and there are no systematic errors. The reduced chisquared can be used to test, in part, this hypothesis:
(15) 
where are the individual measurements, is the weighted mean value, 12 is the number of measurements, and 1 is the number of calculated parameters. Ideally, if the assumption is valid, we should expect 1. Conventionally, the hypotesis is rejected if , where is the critical value corresponding to a probability of less than 5 for the hypothesis to be valid. We found 0.42 for the eclipse depth, and 1.0 for the phase constant, confirming the nondetection of any interepoch variability. 0.42 may suggest that the error bars for the eclipse depth are overestimated, but this is not totally surprising, given that we actively increased them to account for potential uncorrected systematics and biases in the detrending method. Note that the reduced chisquared tests whether the actual dispersion in the measurements is consistent within their error bars, but it is not sensitive to a uniform bias for all measurements, e.g. a constant shift. Hence, it is not sufficient alone to justify the use of the weighted mean as global estimate of a parameter. Additional tests, reported in the Appendices, show that the weighted mean result is very stable for the eclipse depth. The phase constant appears to be more dependent on certain detrending options, in particular background subtraction. In this case, the adopted weighted mean error bar of 2.010 days is a lower limit, valid under some caveats. In the worstcase scenario, the maximum error bar, calculated without scaling when combining multiple measurements, is 710 days.
3.2 Comparison with a previous analysis of the same observations
Our results are consistent within 1 with the ones from a previous analysis reported in Wong et al. (2014) (see Figures 3 and 4). Our error bars are generally larger by a factor 0.81.5 for the eclipse depth (smaller in 1 case) and 1.02.0 for the phase constant compared to the ones in the literature. The factors for the weighted mean eclipse depth and phase constant are 0.9 and 1.4, respectively. Slightly larger error bars are a worthwhile tradeoff for much higher objectivity, which derives from the lack of assumptions about the origin of instrument systematics and their functional forms in our detrending method. We also note that, despite the larger nominal error bars, the dispersions in our best parameter estimates are slightly smaller than the ones calculated from the results reported in Wong et al. (2014) (see Table 4).
The reduced chisquared values inferred from their individual parameter estimates are 0.86 for the eclipse depth, and 4.3 for the phase constant. While the first value is consistent with the hypothesis of a constant transit depth within the quoted error bars, the second value indicates that the analogous hypothesis for the phase constant can be rejected (less than 0.1 probability of being true). This may suggest either that they were able to detect some astrophysical variability of the phase curve’s slope, or that their individual error bars are underestimated by a factor 2. Given that the astrophysical slope is degenerate with other instrumental trends, such as longterm position drift of the telescope and possible thermal heating, it is possible that their individual error bars do not fully accounts for these degeneracies, as the authors themselves state. If this is the case, we observe that their final error bar on the phase constant, derived by a joint fit of all eclipses, could be equally underestimated, because it is not guaranteed that residual systematic errors cancel out over multiple observations as if they were random errors. Note that the joint fit approach is theoretically valid under the same assumptions for which the weighted mean is valid, and the two approaches are expected to lead to very similar results (we checked that this happens in this case). In conclusion, our reanalysis confirms the results reported in Wong et al. (2014) for the eclipse depth, and relative interepoch variability, but not the 4 detection of a nonflat phase curve’s slope during the eclipse, as larger error bars are needed to account for the possible residual systematics.
Eclipse depth  This work  Wong et al. 2014 

Best estimate  (1.57 0.03)10  1.58010 
Dispersion  7.210  8.410 
0.42  0.86  
Phase constant (days)  This work  Wong et al. 2014 
Best estimate  (4.4 2.0)10  6.010 
Dispersion  7.010  11.310 
1.0  4.3 
4 Conclusions
We have applied a blind signalsource separation method to analyze twelve photometric observations of the eclipse of the exoplanet XO3b obtained with Warm Spitzer/IRAC at 4.5 m. The method is an evolution of the pixelICA method proposed and successfully used by our team to analyze real and synthetic transit observations. By adding a wavelet transform of the time series, we extend the applicability of pixelICA to detrend lowS/N observations with instrumental systematics stronger than the astrophysical signal. Wavelet pixelICA results are consistent within 1 with results reported in the literature. They also have smaller dispersions in the eclipse parameters measurements, even including the most recent results that appeared on the arXiv while this paper was under review. While the individual error bars are usually more conservative, as they fully accounts for the possible uncertainties, the final error bar on the eclipse depth is equal or smaller than the ones obtained with other methods discussed in the literature.
No significant interepoch variations are detected over twelve repeated observations in 3 years interval. This is convincing evidence that, with appropriate data detrending methods, transit and eclipse measurements based on Spitzer/IRAC observations can achieve this level of precision and reproducibility, and therefore are useful to characterize the atmospheres of exoplanets. Also, the lack of any detectable astrophysical variability, for the XO3b system, allows to combine multiple observations to increase the accuracy in stellar and planetary parameters.
Appendix A Testing the robustness of results
a.1 The effect of binning
Given the large number of frames (14,912) for each observation, binning data is very useful to decrease the computational time needed to run the MCMCs for the eclipse parameters and scaling coefficients for the independent components (see Section 2.4.7). Some authors suggest that an optimal choice of the binning size can be useful to reduce the noise on the timescale of interest (Deming et al., 2015; Kammer et al., 2015), provided that the theoretical curve is similarly binned as necessary to eliminate bias, and the bin size is not too long to cause significant loss of astrophysical information (Kipping, 2010).
An additional choice for the pixelICA algorithm is whether to bin the pixel time series prior the ICA separation, or to bin the independent components extracted from unbinned pixel time series. We found that the two options are almost equivalent, as the eclipse signals obtained after removing the systematic components from the raw lightcurve are identical (discrepancies 12 order of magnitudes smaller than the fitting residuals), except in cases for which the unbinned ICA separation fails to retrieve an eclipse component. It is worth to note that for the unbinned case, the amplitude of the total noise plus systematics is higher than the the eclipse depth. Thus we decided to bin the individual pixellightcurves prior ICA retrieval.
We compare the results obtained for all the observations with bin sizes of 32 and 64 frames, i.e. 64 and 128 s, respectively. First, we test the gaussianity of fitting residuals by calculating their root mean square (rms) as a function of the bin size, . If fitting residuals are white noise, the rms would scale as . Figure 5 shows that, in both cases, the rms of fitting residuals slightly deviates from the expected behaviour of white noise. Those deviations are smaller for the analysis with bin size of 32 frames.
The parameter results obtained with the two binning choices are all consistent within 0.5 , and in average within 0.16 (see also Figures 6 and 12). Both the error bars and the overall dispersions are smaller for the cases with bin size of 32 frames by factors of 1.4. We consider the results obtained with bin size of 32 frames as our best results.
a.2 Specific wavelet ICA options
The DWT of a time series is specified by the choice of a mother wavelet and the number of levels (see Section 2.3.2). We adopted a Daubechies4 mother wavelet, and onelevel decomposition. For a few observations we tested multiple choices of the mother wavelet, among the Daubechies, Biorthogonal, Symlets and Coiflets families, and number of decomposition levels. We refer to Daubechies (1992); Percival & Walden (2000) for details about the wavelet properties. We found that different choices of the mother wavelet are not significant, e.g. discrepancies in the eclipse signals are 12 orders of magnitudes smaller than the fitting residuals, while level decompositions higher than 1 usually appear to make impossible for ICA to retrieve the eclipse. The difficulties with higherlevel DWTs may arise from subsampling the average coefficients, and the fact that some of the lowfrequency nongaussian components may be smeared over higher levels of detail coefficients.
a.3 About background subtraction
Uncorrected background may bias the normalized amplitude of the eclipse depth, as well as the phase curve’s slope, if background is not constant over time. The typical morphology of a background time series is either a constant function or a slow monotonic drift. The lack of a distinct temporal structure and nongaussianity makes it difficult to disentangle with ICA, as well as other statistical methods. For this reason, we performed ad hoc background subtractions before ICA detrending (see Section 2.4.3). In this Section, we discuss the impact of this step in the analyses, by comparing results obtained with and without background subtraction. These are also used to infer the maximum parameter errors that can be caused by an inappropriate background correction.
Figure 7 shows an example of background time series estimated for one observation. The measured mean background level slightly varies from one observation to the other, but it is always less than 0.6 of the total flux from the system, hence potentially affecting the eclipse depth by less than 10, well below the error bars. Background is also not constant during one observation: it has a small ramp for the first 20 minutes (except for the eclipses extracted from longer observations), then continues to slowly increase.
Figure 8 shows the parameter results obtained with and without background subtraction. The greatest discrepancies are observed for the phase constant, which is systematically higher by 22010 days for the cases without background subtraction, suggesting the presence of uncorrected systematics. It is quite remarkable that our individual error bars automatically account for those systematics, but, given their nonrandom nature, the weighted mean error bars for the phase constant cannot be compared. It is difficult to prove the superiority of the results obtained with background subtraction over the others, as the residuals between the full models and the relevant raw lightcurves (see Section 2.4.7) are very similar for the two cases, e.g. similar amplitudes and time correlations, leading to similar error bars. Slightly smaller parameter dispersions are obtained with background subtraction rather than without, in particular 710 vs 910 for the eclipse depth, and 710 vs 810 days for the phase constant. The higher mean value of the phase constant obtained without background subtraction, i.e. 1310 days, appears to be less likely, as it would require a higher than expected increase in the atmospheric temperature due to stellar irradiation during the eclipse, and/or strong horizontal disomogeneities either in temperature, chemical composition and/or clouds (Cowan & Agol, 2011; Kataria et al., 2013; Agúndez et al., 2014). If we assume that systematics are removed with background subtraction, we can take the weigthed mean as best estimate for the phase constant.
The discrepancies between eclipse depth measurements with and without background subtraction are in the range 1010, and, in average, the eclipse depth is smaller by 410 for the case without background subtraction. Although this is more than the 10 difference expected from the mean background level relative to the mean stellar flux (see discussion in the previous paragraph), in this case, the two weighted means are consistent within 0.5 . We found that the main effect of background on the eclipse depth is due to correlations between the measured eclipse depth at the eclipse center and the phase constant, in terms of Pearson correlation coefficients between the relevant MCMCs. The details of this study are beyond the scope of this paper.
a.4 Using different arrays of pixels
In previous studies, we found that, for highS/N transit observations, pixelICA performances are only slightly dependent on the choice of the array (Morello et al., 2014, 2015). The 33 array is too narrow compared to the Spitzer/IRAC Point Spread Function, and for this reason, it is not photometrically stable (although it usually leads to consistent results). Larger arrays are more photometrically stable, but they include more noise. Also, noisier pixellightcurves, appear to increase the uncertainties in the ICA decomposition (, see Section 2.4.8), so that the smallest error bars were usually obtained using the 55 array.
In this work, we analyzed all datasets with two different choices of the pixelarray, i.e. 55 and 77. For the analyses with the 77 array, we only adopted the 64 frames bin size, for which the MCMC fitting is faster. Figure 9 compares the results obtained with the two different arrays and the same bin size. The two sets of results are consistent well within 1 , but error bars obtained with the 77 array are 11.5 times larger than the ones obtained with the 55 array, despite the residuals in the fits are similar and often smaller for the 77 cases. This confirms the conclusions obtained from previous analyses, in particular:

the 55 array leads to smaller error bars than other squared arrays of pixels;

parameter results from different arrays are consistent well within 1 .
In this case, the choice of a less optimal array increases the error bars more significantly than in our previous analyses, most likely because of the lower S/N of the observations analyzed here.
Figure 10 compares the results obtained with and without background subtraction for the case of 77 array. The same considerations discussed in Appendix A.3 for the 55 array are valid for the 77 array.
a.5 Breaking degeneracy with the phase constant
At the end of Appendix A.3 we revealed the existence of a correlation between the eclipse depth and phase constant parameter, which is stronger for data with a higher slope, such as the cases without background subtraction. This may affect the eclipse depth estimate in case of residual systematics with a slope, e.g. uncorrected background. In our analyses, the maximum bias on the eclipse depth due to correlations with this kind of systematics is 410, provided the systematic slope is not larger than our individual error bars.
Here we test the consequences of adopting zero phase constant while fitting for the eclipse depth. This is equivalent to the assumption that the phase curve is flat before and after the eclipse, and the slope is entirely due to instrumental effects. Compared to the cases with the phase constant as free parameter, fitting residuals with zero phase constant are not significantly larger. Figure 11 compares results for the eclipse depth with free and zero phase constant, with and without background subtraction. The weighted mean results are reported in Figure 12.
We note that differences in eclipse depth measurements with and without background subtraction are smaller for the zero phase constant models, suggesting that eclipse depth measured with zero phase constant models are less affected by residual systematics with a slope. Free phase constant models are valid in a more general context, as, differently from the zero phase constant models, they approximate planet’s flux variability during the observation. The consistency between the results obtained with the two classes of models indicates that, in this case, planet’s flux variability in the proximity of the eclipse is smaller than the error bars.
a.6 Robustness of the eclipse depth measurement
Figure 12 reports the weighted mean eclipse depth estimated for all the tests discussed in the Appendices. Note that they are mutually consistent at the 0.5 level, and the range is 610.
Footnotes
 slugcomment: Not to appear in Nonlearned J., 45.
References
 Agúndez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, A&A, 564, A73
 Ballard, S., Charbonneau, D., Deming, D., Knutson, H. A., Christiansen, J. L., Holman, M. J., Fabrycky, D., Seager, S., & A’Hearn, M. F., 2010, PASP, 122, 1341
 Ballerini, P., Micela, G., Lanza, A. F., & Pagano, I. 2012, A&A, 539, A140
 Beaulieu, J. P., Carey, S., Ribas, I., & Tinetti, G. 2008, ApJ, 677, 1343
 Beaulieu, J. P., Tinetti, G., Kipping, D. M., Ribas, I., Barber, R. J., Cho, J. Y. K., Polichtchouk, I., Tennyson, J., Yurchenko, S. N., Griffith, C. A., Batista, V., Waldmann, I. P., Miller, S., Carey, S., Mousis, O., Fossey, S. J., & Aylward, A. 2011, ApJ, 731, 16
 Berta, Z. K., Charbonneau, D., Bean, J., Irwin, J., Burke, C. J., Désert, J. M., Nutzman, P., & Falco, E. E. 2011, ApJ, 736, 12
 Brown, T. M. 2001, ApJ, 553, 1006
 Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377
 Charbonneau, D., Allen, L. E., Megeath, S. T., Torres, G., Alonso, R., Brown, T. M., Gilliland, R. L., Latham, D. W., Mandushev, G., O’Donovan, F. T., & Sozzetti, A. 2005, ApJ, 626, 523
 Cowan, N. B., & Agol, E. 2011, ApJ, 726, 82
 Cowan, N. B., Machalek, P., Croll, B., Shekhtman, L. M., Burrows, A., Deming, D., Green, T., & Hora, J. L. 2012, ApJ, 747, 82
 Daubechies, I. 1992, SIAM, Ten Lectures on wavelets, ISBN: 9780898712742
 Deming, D., Wilkins, A., McCullough, P., Burrows, A., Fortney, J. J., Agol, E., DobbsDixon, I., Madhusudhan, N., Crouzet, N., Désert, J. M., Gilliland, R. L., Haynes, K., Knutson, H. A., Line, M., Magic, Z., Mandell, A. M., Ranjan, S., Charbonneau, D., Clampin, M., Seager, S., & Showman, A. P. 2013, ApJ, 774, 95
 Deming, D., Knutson, H. A., Kammer, J., Fulton, B. J., Ingalls, J., Carey, S., Burrows, A., Fortney, J. J., Todorov, K., Agol, E., Cowan, N. B., Désert, J. M., Fraine, J., Langton, J., Morley, C., & Showman, A. P. 2015, ApJ, 805, 132
 Désert, J. M., Lecavelier des Etangs, A., Hébrard, G., Sing, D. K., Ehrenreich, D., Ferlet, R., & VidalMadjar, A. 2009, ApJ, 699, 478
 Donoho, D. L. 1995, ITIT, 41, 613
 Evans, T. M., Aigrain, S., Gibson, N., Barstow, J. K., Amundsen, D. S., Tremblin, P., & Mourier, P. 2015, arXiv:1504.05942v1
 Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
 Haario, H., Laine, M., Mira, A., & Saksman, E. 2006, Statistics and Computing, 16, 339
 Hirano, R., Narita, N., Sato, B., Harakawa, H., Fukui, A., Aoki, W., & Tamura, M. 2011, PASJ, 63, L57
 Ingalls, J. G., Krick, J. E., Carey, S. J., Stauffer, J. R., Lowrance, P. J., Grillmair, C. J., Buzasi, D., Deming, D., DiamondLowe, H., Evans, T. M., Morello, G., Stevenson, K. B., Wong, I., Capak, P., Glaccum, W., Laine, S., Surace, J., StorrieLombardi, L. 2016, arXiv:1601.05101v1
 Inuso, G., La Foresta, F., Mammone, M., & Morabito, F. C. 2007, Int. J. Conf. Neural Netw., 20, 1524
 JohnsKrull, C. M., McCullough, P. R., Burke, C. J., Valenti, J. A., Janes, K. A., Heasley, J. N., Prato, L., Bissinger, R., Fleenor, M., Foote, C. N., GarciaMelendo, E., Gary, B. L., Howell, P. J., Mallia, F., Masi, G., & Vanmunster, T. 2008, 677, 657
 Kammer, J. A., Knutson, H. A., Line, M. R., Fortney, J. J., Deming, D., Burrows, A., Cowan, N. B., Triaud, A. H. M. J., Agol, E., Désert, J. M., Fulton, B. J., Howard, A. W., Laughlin, G. P., Lewis, N. K., Morley, C. V., Moses, J. I., Showman, A. P., & Todorov, K. O. 2015, ApJ, 810, 118
 Karande, K. J., & Talbar, S. 2014, Springer India, Independent Component Analysis of Edge Information for Face Recognition, ISBN:9788132215110
 Kataria, T., Showman, A. P., Lewis, N. K., Fortney, J. J., Marley, M. S., & Freedman, R. S. 2013, ApJ, 767, 76
 Kipping, D. K. 2010, MNRAS, 408, 1758
 Knutson, H. A., Madhusudhan, N., Cowan, N. B., Christiansen, J. L., Agol, E., Deming, D., Désert, J. M., Charbonneau, D., Henry, G. W., Homeier, D., Langton, J., Laughlin, G., & Seager, S. 2011, ApJ, 735, 27
 Knutson, H. A., Lewis, N., Fortney, J. J., Burrows, A., Showman, A. P., Cowan, N. B., Agol, E., Aigrain, S., Charbonneau, D., Deming, D., Désert, J. M., Henry, G. W., Langton, J., & Laughlin, G. 2012, ApJ, 754, 22
 Knutson, H. A., Fulton, B. J., Montet, B. T., Kao, M., Ngo, H., Howard, A. W., Crepp, J. R., Hinkley, S., Bakos, G. Á., Batygin, K., Johnson, J. A., Morton, T. D., & Muirhead, P. S. 2014, ApJ, 785, 126
 La Foresta, F., Mammone, N., & Morabito, F. C. 2012, in Neural Nets, Vol. 3931, ed. B. Apolloni (Berlin: Springer), 78
 Lagarias, J. C., Reeds, J. A., Wright, M. H., & Wright, P. E. 1998, SIAM Journal of Optimization, 9, 112
 Lewis, N. K., Knutson, H. A., Showman, A. P., Cowan, N. B., Laughlin, G., Burrows, A., Deming, D., Crepp, J. R., Mighell, K. J., Agol, E., Bakos, G. A., Charbonneau, D., Désert, J. M., Fischer, D. A., Fortney, J. J., Hartman, J. D., Hinkley, S., Howard, A. W., Johnson, J. A., Kao, M., Langton, J., & Marcy, G. W. 2013, ApJ, 766, 95
 Lin, J., & Zhang, A. 2005, NDTI, 38, 421
 Mammone, N., La Foresta, F., & Morabito, F. C. 2012, IEEE Sensors Journal, 12, 533
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171
 MoralesCaldéron, M., Stauffer, J. R., Davy Kirkpatrick, J., Carey, S., Gelino, C. R., Barrado y Navascués, D., Rebull, L., Lowrance, P., Marley, M. S., Charbonneau, D., Patten, B. M., Megeath, S. T., & Buzasi, D. 2006, ApJ, 653, 1454
 Morello, G., Waldmann, I. P., Tinetti, G., Peres, G., Micela, G., & Howarth, I. D. 2014, ApJ, 786, 22
 Morello, G., Waldmann, I. P., Tinetti, Howarth, I. D., G., Micela, G., & Allard, F. 2015, ApJ, 802, 117
 Morello, G. 2015, ApJ, 808, 56
 Moussaoui, R., Rouat, J., & Lefebvre, R. 2006, In: IEEE International Conference on Acoustic, Speech and Signal Processing, vol. 5, pp. V645âV648 (2006)
 Percival, D. B., & Walden, A. T. 2000, Cambridge University Press, Wavelet Methods for Time Series Analysis, ISBN: 9780521685085
 Tessenyi, M., Tinetti, G., Savini, G., & Pascale, E. 2013, Icarus, 226, 1654
 Tinetti, G., Deroo, P., Swain, M. R., Griffith, C. A., Vasisht, G., Brown, L. R., Burke, C., & McCullough, P. 2010, ApJ, 712, L139
 Stein, C. M. 1981, AnSta, 9, 1135
 Stevenson, K. B., Harrington, J., Nymeyer, S., Madhusudhan, N., Seager, S., Bowman, W. C., Hardy, R. A., Deming, D., Rauscher, E., & Lust, N. B. 2010, Nature, 464, 1161
 Stevenson, K. B., Harrington, J., Fortney, J. J., Loredo, T. J., Hardy, R. A., Nymeyer, S., Bowman, W. C., Cubillos, P., Bowman, M. O., & Hardin, M. 2012, ApJ, 754, 136
 Swain, M. R., Tinetti, G., Vasisht, G., Bouwman, J., Chen, P., Yung, Y., Deming, D., & Deroo, P. 2009, ApJ, 690, L114
 Torres, G., Fischer, D. A., Sozzetti, A., Buchhave, L. A., Winn, J. N., Holman, M. J., & Carter, J. A. 2012, ApJ, 757, 161
 Waldmann, I. P. 2012, ApJ, 747, 12
 Waldmann, I. P., Tinetti, G., Deroo, P., Hollis, M. D. J., Yurchenko, S. N., & Tennyson, J. 2013, ApJ, 766, 7
 Waldmann, I. P. 2014, ApJ, 23, 780
 Waldmann, I. P., Tinetti, G., Rocchetto, M., Barton, E. J., Yurchenko, S. N., & Tennyson, J. 2015, ApJ, 802, 107
 Waldmann, I. P., Rocchetto, M., Tinetti, G., Barton, E. J., Yurchenko, S. N., & Tennyson, J. 2015, ApJ, 813, 13
 Winn, J. N., Holman, M. J., Torres, G., McCullough, P., JohnsKrull, C., Latham, D. W., Shporer, A., Mazeh, T., GarciaMelendo, E., Foote, C., Esquerdo, G., & Everett, M. 2008, ApJ, 683, 1076
 Wong, I., Knutson, H. A., Cowan, N. B., Lewis, N. K., Agol, E., Burrows, A., Deming, D., Fortney, J. J., Fulton, B. J., Langton, J., Laughlin, G., & Showman, A. P. 2014
 Zellem, R. T., Lewis, N. K., Knutson, H. A., Griffith, C. A., Showman, A. P., Fortney, J. J., Cowan, N. B., Agol, E., Burrows, A., Charbonneau, D., Deming, D., Laughlin, G., & Langton, J. 2014, ApJ, 790, 53
 Zhao, C.H., Liu, J., & Sun, J.D. 2006, Joumal of Electronics & Infornation Technology, 28, 1565