Evaluating Climate Variability Of The Canonical Hot Jupiters Hd 189733b & Hd 209458b Through Multi-epoch Eclipse Observations

TL;DR

This study analyzes multi-epoch eclipse data of two hot Jupiters to assess their atmospheric variability, finding no significant changes and setting upper limits on possible IR variability over several years.

Contribution

It provides the first constraints on temporal IR variability of HD 189733b and HD 209458b using multi-epoch observations spanning years.

Findings

No statistically significant variability detected in IR emission.

Upper limits on variability amplitudes are established for both planets.

Observed eclipse timing offsets are consistent with hotspot models.

Abstract

Here we present the analysis of multi-epoch secondary eclipse observations of HD 189733b and HD 209458b as a probe of temporal variability in the planetary climate using both Spitzer channels 1 and 2 (3.6 and 4.5 um). Constraining temporal variability will inform models and identify physical processes occurring at either length scales too small to directly observe or at pressure levels that are inaccessible to transit observations. We do not detect statistically significant variability and are able to place useful upper limits on the IR variability amplitudes in these atmospheres. There are very few planets with multi-epoch observations at the required precision to probe variability in dayside emission. The observations considered in this study span several years, providing insight into temporal variability at multiple timescales. In the case of HD 189733b, the best fit eclipse depths for the channel 2 observations exhibit a scatter of 102 ppm about a median depth of 1827 ppm and in channel 1 exhibit a scatter of 88 ppm about a median depth of 1481 ppm. For HD 209458b, the best fit eclipse depths for the channel 2 observations exhibit a scatter of 22 ppm about a median depth of 1406 ppm and in channel 1 exhibit a scatter of 131 ppm about a median depth of 1092 ppm. The precision and scatter in these observations allow us to constrain variability to less than (5.6% and 6.0%) and (12% and 1.6%) for channels (1,2) of HD 189733b and HD 209458b respectively. There is a difference in the best fit eclipse timing compared to the predicted time consistent with an offset hotspot as predicted by GCMs and confirmed in previous phase curve observations.

Figures

Tables

Table 1. Table 3 Ephemerides

Parameter	HD 189733b	HD 209458b
$T_{\star }$(K)	5040$\pm$50	6065$\pm$50
$M_{\star }$ (M${}_{\mathrm{Sun}}{}^{})$	0.806$\pm$0.048	1.131$\pm$0.026
$R_{\star }$ (R${}_{\mathrm{Sun}}{}^{})$	0.756$\pm$0.018	1.155$\pm$0.015
$M_p$ (M${}_{\mathrm{Jup}}{}^{})$	1.144$\pm$0.056	0.690$\pm$0.024
$R_p$ (R${}_{\mathrm{Jup}}{}^{})$	1.138$\pm$0.027	1.359$\pm$0.015
$R_p/R_{\star }$	0.024122$\pm$5.8$\times$10-5	0.014607$\pm$2.4$\times$10-5
log(g) (log(c/s2))	3.339$\pm$0.03	2.969$\pm$0.0187
Period (days)	2.21857567$\pm$1.5$\times$10-7	3.52474859$\pm$ 3.8$\times$10-7
$i$ (∘)	85.7100$\pm$0.0023	86.710$\pm$ 0.05
m$\sin i$ (M${}_{\mathrm{Jup}}{}^{})$	1.140 $\pm$0.056	0.689$\pm$0.024
$a/R_{\star }$	8.84$\pm$0.27	8.81$\pm$0.186

Table 2. Table 4 Results of Best Fit Eclipse Depth and Time

AOR	Eclipse Depth	Eclipse Time	SDNR	${\beta }_{\mathrm{red}}$
	(ppm)	O-C (min)
HD 209458b Ch 2
38703616	1370 $\pm$ 37	2.22 $\pm$ 1.04	0.003240	1.43
48014336	1401 $\pm$ 42	1.44 $\pm$ 0.83	0.003085	1.17
50489088	1424 $\pm$ 60	3.03 $\pm$ 0.90	0.003085	1.16
50490880	1380 $\pm$ 58	1.62 $\pm$ 0.79	0.003053	1.27
50490624	1443 $\pm$ 56	3.14 $\pm$ 0.82	0.003060	1.27
50490368	1417 $\pm$ 45	2.60 $\pm$ 0.96	0.003151	1.05
50490112	1420 $\pm$ 52	2.81 $\pm$ 0.72	0.003125	1.18
50489856	1424 $\pm$ 47	1.57 $\pm$ 0.65	0.003116	1.09
HD 209458b Ch 1
48013824	1050 $\pm$ 73	4.37 $\pm$ 1.61	0.004995	2.59
50496512	909 $\pm$ 60	1.48 $\pm$ 1.27	0.004873	1.99
50496000	1162 $\pm$ 77	3.38 $\pm$ 1.39	0.004834	1.96
50496256	1209 $\pm$ 67	2.81 $\pm$ 1.13	0.004805	1.95
50495744	845 $\pm$ 92	4.21 $\pm$ 1.84	0.005055	2.29
50495488	1008 $\pm$ 76	1.81 $\pm$ 1.37	0.004946	1.78
50494976	1049 $\pm$ 78	2.48 $\pm$ 1.30	0.004898	2.03
50494208	1209 $\pm$ 72	1.61 $\pm$ 1.36	0.0049606	2.03
50493440	1301 $\pm$ 76	2.20 $\pm$ 1.15	0.004983	1.87
50492928	994 $\pm$ 105	6.82 $\pm$ 2.61	0.004939	2.74
41629440	1216 $\pm$ 77	1.28 $\pm$ 2.32	0.006171	4.14
41628416	1128 $\pm$ 72	1.99 $\pm$ 2.72	0.006168	3.68
HD 189733b Ch 2
38390784	1812 $\pm$ 50	1.52 $\pm$ 0.69	0.004596	1.65
38390016	1806 $\pm$ 64	0.31 $\pm$ 0.97	0.004698	1.74
50489344	1716 $\pm$ 55	0.96 $\pm$ 0.45	0.004631	1.06
50488832	1815 $\pm$ 61	-0.37 $\pm$ 0.61	0.004606	1.41
50488576	1984 $\pm$ 75	0.31 $\pm$ 0.64	0.004630	1.31
50488320	1774 $\pm$ 50	0.24 $\pm$ 0.55	0.004572	1.12
50488064	2016 $\pm$ 65	-0.02 $\pm$ 0.48	0.004654	1.00
HD 189733b Ch 1
41592320	1574 $\pm$ 305	1.40 $\pm$ 6.0	0.004228	4.55
41591296	1431 $\pm$ 86	0.29 $\pm$ 1.74	0.003382	2.13
40152064	1271 $\pm$ 310	0.16 $\pm$ 4.08	0.001708	3.76
50495232	1514 $\pm$ 105	-0.37 $\pm$ 1.35	0.003370	1.51
50494464	1473 $\pm$ 73	0.81 $\pm$ 0.93	0.003298	1.50
50493696	1498 $\pm$ 53	0.80 $\pm$ 0.67	0.003266	1.54
50493184	1491 $\pm$ 115	-0.74 $\pm$ 1.37	0.003306	1.85

Equations

$$x_{\mathrm{cen}}=\frac{{\sum }_{j,k}(I_{jk}j)}{{\sum }_{j,k}{I}_{jk}};y_{\mathrm{cen}}=\frac{{\sum }_{j,k}(I_{jk}k)}{{\sum }_{j,k}{I}_{jk}}.$$

$$\tilde{\beta }=\frac{{\left(\sum P_i\right)}^2}{\sum \left(P_i^2\right)},$$

$$R_{model}=1-a_1\times e^{(-\frac{t}{a_2})}.$$

$$\Delta S^t=\sum _{i=1}^n{c}_i{\hat{P}}_i^t+DE(t)+ft+g{t}^2+h,$$

Full text

††slugcomment: Accepted for Publication in AJ December 17, 2019

1. Introduction

Studying exoplanet atmospheres is challenging since we often are not able to spatially resolve them. They are typically too faint to disentangle their light from that of their much brighter host star. Transiting exoplanets provide a unique opportunity to study exoplanet atmospheres in spite of this challenge. Disc integrated secondary eclipse observations provide valuable information on temperature, albedo, and chemical composition averaged over the entire hemisphere. Primary transits give a limb averaged atmospheric molecular spectrum and pressure - temperature profile. However, we must employ other novel techniques to begin to probe physical and chemical processes and structures at smaller scales. Temporal variability is one way to probe weather features and their their movements without spatially resolving them. Variability is also a tool to study portions of the atmosphere inaccessible via other observational techniques. Transmission spectroscopy observations in the near infrared typically probe millibar pressure levels high in the atmosphere. Processes deeper in the atmosphere are obscured due to the long path lengths of the transit geometry. Time varying changes in emission may be indicative of processes occurring deeper in the atmosphere that are, otherwise, unobservable with current methods and insturments.

Close-in giant planets (hot Jupiters) are interesting objects with which to study atmospheric dynamics. They are assumed to be tidally locked based on their short orbital periods (¡5 days) and minimal separation from their host stars. As a result, they experience a constant radiative forcing on their permanent dayside. The combination of the rotation rate and radiative forcing of these planets are predicted to produce large scale weather structures unlike anything in our solar system.

There have been many efforts to model the circulation and temperature structure of hot Jupiters with many of these efforts based on the properties of the two most well studied targets, HD 189733b and HD 209458b (HD 189733b& HD 209). Simulations by Showman & Guillot (2002); Cho et al. (2003) and Menou et al. (2003) predicted that the Rossby deformation radius and Rhines scale of hot Jupiters should be comparable to the planetary radius resulting in atmospheric dynamics comprised of a few jets and large scale polar vortices. They predicted that vortices would be large enough that their migration could have an affect on the observed eclipse depth. 2007RauscherVar demonstrated that, based on those model predictions, changes in eclipse depth as great as 20% could be observed. Later simulations by Showman et al. (2009) predicted far more stable atmospheric structures. Variability was predicted to be less than several percent for HD 189733b at 8 $\mu$m and the result of wave dynamics deep in the atmosphere rather than migrating vortices.

Agol et al. (2010a) observationally probed model predictions with a multi - epoch set of transit and eclipse observations of HD 189733b at 8 µm. The results of that work placed an upper limit on eclipse depth variability at 2.7%, effectively ruling out larger predictions. There have been few other attempts to study temporal variability due to the lack of multiple observations over a sufficient temporal baseline.

We choose as our targets for this study HD 189733b and HD 209458b both because of their suitability for observations with the Spitzer observatory and because they are the only targets with a large number of observations with a single instrument spanning a number of years. These two planets are the most well studied and characterized exoplanets to date. Both orbit bright stars, Kmag of 5.5 and 6.3 for HD 189733b and HD 209458b respectively, making them ideal targets for characterization and have previously observed secondary eclipse depths ¿1000 ppm in both warm Spitzer channels 1 and 2.

HD 189733b is arguably the most thoroughly studied exoplanet. It has been observed with photometric (e.g. Ehrenreich et al., 2007; Désert et al., 2009, 2011b) and spectroscopic (e.g. Gibson et al., 2012) transits as well as secondary eclipse (e.g. Charbonneau et al., 2008a; Agol et al., 2010a) and phase curve (Knutson et al., 2007, 2009, 2012) observations across a multitude of wavelengths. In addition to full orbit phase curves at multiple wavelengths, it is the only planet to have been mapped via the eclipse mapping technique (Majeau et al., 2012; de Wit et al., 2012; Rauscher et al., 2018). HD 209458b was the first transiting exoplanet discovered (Charbonneau et al., 2000). It has also been observed in both transit and eclipse geometries photometrically (e.g. Charbonneau et al., 2008b; Diamond-Lowe et al., 2014) and spectroscopically (e.g. Deming et al., 2013; Line et al., 2016) , as well as full orbit phase curve observations (Zellem et al., 2014).

The large amount of observational and theoretical effort invested into these two canonical hot Jupiters make them perfect choices for a study of this nature. There are eclipse observations from multiple programs spanning years from which to build the necessary temporal baseline. Additionally, their orbital properties are well studied and well constrained allowing for the detection of small signals in their lightcurves due to planetary thermal structure that can be clearly distinguished from other potential sources such as orbital eccentricity.

Understanding variability will be essential in current and future attempts to leverage multiple epoch observations to build precision for high resolution techniques. When stacking observations, whether it be transmission spectroscopy or emission photometry, one makes the assumption that each observation is an independent measurement of a constant signal. Showing that any temporal variability is below the level of precision of your observations allows them to be combined in analysis. If this set of Spitzer observations are not sensitive to temporal variability they can be combined to achieve the precision necessary to perform eclipse mapping.

Eclipse mapping (Majeau et al., 2012; de Wit et al., 2012) uses the deviations to the shape of ingress and egress caused by a non uniform dayside temperature distribution to create two dimensional thermal maps. In general, one needs $\sim$10 points over the ingress/egress with a precision of at least a tenth of the eclipse depth to achieve the necessary resolution to map hot Jupiters. The only way to achieve that precision with current observatories is to stack multiple observations. In the James Webb Space Telescope (JWST) era, observing time will be expensive and full orbit phase curves may prove too costly. However, mapping may be done much more efficiently with as little as two eclipse observations. Understanding temporal variability then is crucial in knowing if stacking multi-epoch observations is justifiable in that context. The observations considered here give the best current insight into orbit to orbit changes in hot Jupiter atmospheres.

2. Observations

Here we consider all existing secondary eclipse observations of HD 189733b and HD 209458b in channels 1 and 2 (3.6 or 4.5$\,\mu{\rm m}$ bandpass) . The observations analyzed here are part of Programs 60021 (PI: H. Knutson), 10103 (PI: N. Lewis), 90186 (PI: K. Todorov), and 70100 (PI: M. Swain). The details of each Astronomical Observing Request (AOR) are displayed in Table 3.1 and 3.1. All of the observations were carried out in sub-array mode ($32\times 32$ pixels, $39{\tt"}\times 39{\tt"}$). Observations in Program 10103 utilize a 30 minute peak-up observation preceding them to stabilize the image on the detector ‘sweet spot’ and decreases the likelihood of a ramp in the data (Ingalls et al., 2012). Program 60021 observations are full orbit phase curves from which we extract the AORs containing the eclipse from the larger data set. Both HD 189733 and HD 209458 are bright targets with Kmag of 5.5 and 6.3 respectively. As a result, frame times for all of the observations are 0.1 seconds. The two exceptions are the observations of HD 209458b in channel 2 that were conducted with a 0.4 second frame time and the observations of HD 189733b from program 70100 that were conducted such that the stellar centroid was near the corner of the pixel so that a longer (0.4 s) exposure time could be used. This technique proved problematic in that the systematics were not easily corrected with existing techniques and, as a result, several of these observations were discarded from this study.

3. Data Analysis Methods

3.1. Photometric Extraction

For each AOR we began with Basic Calibrated Data (BCD) available on the Spitzer Heritage Archive. Each BCD file contains a cube of 64 frames of $64\times 64$ pixels. Each set of 64 images comes as a single FITS file with a time stamp corresponding to the start of the first image. We determine the time of each frame in the set by adding the appropriate multiple of the frame time to the time stamp of the first image. Each frame was corrected for bad pixels or NaN values by masking the invalid pixels. Each frame is background subtracted based on the median pixel value outside of a 10 x 10 pixel box around the stellar centroid. Pixel values in the background region that deviate from the mean by more than 5 $\sigma$ are clipped and the median value remaining is taken to be the background. Stellar centroid positions for each frame are determined by a first moment center of mass calculation:

[TABLE]

We perform two-dimensional Gaussian Centroiding as well but find the first moment calculation to be more stable and provides less scatter in the extracted photometry.

The noise pixel parameter ($\tilde{\beta}$) (Lewis et al., 2013; Kilpatrick et al., 2017) is calculated for each frame as

[TABLE]

where each $P_{i}$ is the response measured in each pixel across the frame. We also save the 5 x 5 array of background subtracted pixel values about the stellar centroid for each frame. Each array is normalized such that they sum to one. We then perform aperture photometry about the stellar centroid using the aperture_photometry function from the Astropy package Photutils. Circular apertures of fixed radii ranging from 1.8 - 2.8 pixels and variable radii apertures ranging from $\beta_{\rm pix}\times$ {0.7…1.2} and $\beta_{\rm pix}+$ {-0.6…1.4}. The resultant time series photometry is then filtered by removing any points that deviate by more than 3 $\sigma$ from the median values of flux or x,y position.

3.2. Systematics Correction and Model Fitting

The intrapixel sensitivity variation (Ingalls et al., 2012), the change in measured flux as a function of stellar centroid position and methods of correction, are well documented (e.g. Ingalls et al., 2016; Kilpatrick et al., 2017). Here, we employ two independent methods of correction for this systematic variation; the Nearest Neighbors method (NNBR), otherwise known as Gaussian Kernel Regression with data (Lewis et al., 2013) and Pixel Level Decorrelation (PLD) (Deming et al., 2015).

Each eclipse fit was based on the model of Mandel & Agol (2002) for a uniform occultation implemented in python by the BATMAN package (Kreidberg, 2015). The orbital, stellar, and planetary parameters listed in table 3 were used as input fixed parameters to the model. Spitzer IRAC data is known to have an exponential ramp in flux over the first 30–60 minutes of observing; however the peak-up technique has alleviated this problem to some extent (Kilpatrick et al., 2017; Ingalls et al., 2012). As a precaution, each data set was fit in several ways: without alteration, trimmed at 20 minutes from the beginning of the observation, and with an exponential ramp model of the form

[TABLE]

PLD has a quadratic visit long temporal variation included in all cases rather than the exponential ramp and combines the systematics and astrophysical models by a linear expression

[TABLE]

where the $\hat{P}^{t}_{i}$’s are the response of the $i^{th}$ pixel in a 3 $\times$ 3 grid around the centroid of the image and $DE(t)$ is the eclipse model and $\Delta S^{t}$ is the total flux, including the astrophysical signal and all systematics, at each time $t$.

For each AOR the eclipse model and systematics model were combined and best fit values for all free parameters were determined using a non-linear, least squares fitting algorithm. The standard deviation of the normalized residuals (SDNR) was used as a metric for selecting the best fit out of the 30 different apertures for each AOR. The Bayesian Information Criterion (BIC) was used to determine the appropriate treatment between various trimming and ramp model options. We find that no trimming or exponential ramp are necessary in all cases of Channel 2 observations and the quadratic ramp with PLD is favored in all Channel 1 observations. The results from the best fit aperture were passed to a Markov Chain Monte Carlo implemented by emcee (Foreman-Mackey et al., 2013) to derive uncertainties of each free parameter. We use a number of walkers at least twice the number of free parameters and run for $10^{5}$ steps per walker before testing for convergence using Gelman Rubin statistics with a threshold for acceptance of 1.01 (Gelman & Rubin, 1992).

4. Results

Here we present the best fit eclipse depths and center of eclipse times, along with their corresponding uncertainty, for each of the observations. The results presented in Table 4 are derived using the PLD method with a quadratic temporal term. No additional trimming or ramps are modeled. The best fit aperture was consistently a fixed radius aperture between 2.2 and 2.4 pixels. We include in Table 4 the standard deviation of the normalized residuals (SDNR) and the $\beta_{\rm red}$ factor, as defined in Gillon et al. (2010), as a measure of correlated noise remaining in the data after systematic corrections. Included in Figure 1 are representative fits for each of the two targets in each of the two channels with the systematics removed and binned at two minute intervals. In the case of HD 189733b, the best fit eclipse depths for the channel 2 observations exhibit a scatter of 102 ppm about a median depth of 1827 ppm and in channel 1 exhibit a scatter of 88 ppm about a median depth of 1481 ppm. For HD 209458b, the best fit eclipse depths for the channel 2 observations exhibit a scatter of 22 ppm about a median depth of 1406 ppm and in channel 1 exhibit a scatter of 131 ppm about a median depth of 1092 ppm.

As shown in Figures 2 and 3, the scatter in the data is consistent with the uncertainty in each measurement in all cases except for the channel 1 HD209 data. Given the number of observations at the achieved precision, we would expect the resultant eclipse depths to represent a sampling of a distribution represented by the shaded areas in Figures 2 and 3. The channel 1 observations of HD209, which exhibit more scatter than indicated by the error bars, is thought to be due to unresolved systematics. This is further supported by the larger values of $\beta_{\rm red}$ in channel 1 observations in comparison to channel 2.

4.1. Constraints on Temporal Variability

As a probe of periodic variability in the eclipse depths we calculate the absolute difference in eclipse depths for each pair of observations. The absolute difference is plotted against the time between observations. Small numbers of observations do not lend themselves to more thorough methods of detecting power at certain time scales such as Lomb-Scargle periodograms, however, this simple approach would show a spike in absolute difference at any relevant timescale. The apparent random scatter of the data points in Figure 4 does not indicate any periodic structure in any of the observations. Given that we find no evidence for variability in our observations we compare the standard deviation and the magnitude of the eclipse depths to constrain variability to less than (5.6% and 6.0%) and (12% and 1.6%) for channels (1,2) of HD 189733b and HD 209458b respectively. We synthesized a periodic signal over a grid of varying periods and amplitudes to demonstrate that our observations would be sensitive to variability above the levels of our constraint. We then simulated measurements using the relevant uncertainty at each of the observation times. We evaluate the delta BIC for a constant signal versus a periodic signal at each grid point averaged over 100 trials. We find that the detection threshold is not sensitive to period over the timescales probed by these observations; only amplitude. The delta BIC is noisy and insignificant until amplitudes equivalent to our previously stated constraints after which it increases rapidly with signal amplitude. This metric provides evidence that if a signal were present our observations would be sensitive to it.

4.2. Stellar Variability

Stellar variability in the form of spots or plages can affect the measured planetary transit signal, potentially altering one’s interpretation of its atmospheric composition, particularly for bright targets with large transit signals (Pont et al., 2008; Silva-Valio, 2008; Czesla et al., 2009; Wolter et al., 2009; Agol et al., 2010b; Berta et al., 2011; Carter et al., 2011; Désert et al., 2011a; Sing et al., 2011; Fraine et al., 2014; McCullough et al., 2014; Oshagh et al., 2014; Damasso et al., 2015; Barstow et al., 2015; Zellem et al., 2015, 2017; Cauley et al., 2017, 2018; Rackham et al., 2017, 2018, 2019; Morris et al., 2018). HD 189733 is an active K0 star which has been observed to vary by as much as $\pm$ 1.5% at visible wavelengths with a rotation period of 11.95 days (Knutson et al., 2012, and references within) whose variability potentially impacts the observed transit depth of the planet (McCullough et al., 2014; Knutson et al., 2012). Although the amplitudes of these variations are reduced as our Spitzer/IRAC observations are in the infrared (Oshagh et al., 2014; Rackham et al., 2017; Morris et al., 2018) and of eclipses (Zellem et al., 2017), we examine the amplitude of these changes in comparison to the scatter and uncertainty in our eclipse depth measurements.

We employ ground-based monitoring spanning the 2009 to 2014 observing seasons with the Tennessee State University 0.8 m Automated Photoelectric Telescope (APT) at Fairborn Observatory in southern Arizona (e.g., Henry, 1999; Eaton et al., 2003, Fig. 5). However, most of our Spitzer /IRAC observations occur outside the APT observing season (Fig. 5). While previous studies have interpolated APT monitoring to their Spitzer observations (Knutson et al., 2012) using an activity model (Aigrain et al., 2012) to fit the APT data and interpolated it to their Spitzer observations, we conservatively take the peak-to-trough Strömgren $b$+$y$ variability observed by APT from 2009 to 2014 (4.7%) and interpolate it to the infrared using the scaling presented in Knutson et al. (2012) (1.6%). Using Equation 7 in Zellem et al. (2017), we estimate the effect of HD 189733’s variability on the observed eclipse depths and find that variability does not statistically impact HD 189733b’s eclipse depths as the changes in its measured eclipse depths induced by variability are less than the measurement uncertainties of the eclipses themselves (Fig. 5).

5. Discussion

5.1. Model Predictions

Given that our assessment of atmospheric variability is in direct relation to the planet’s atmospheric dynamics, we compare our Spitzer variability estimates to predictions from three-dimensional general circulation models (GCMs). In particular, we use GCM simulations of HD 189733b from Showman et al. (2009) and HD 209458b from Kataria et al. (2016), respectively, which utilize the Substellar and Planetary Radiation and Circulation (SPARC) model (Showman et al., 2009). The SPARC model couples the general circulation model maintained at the Massachusetts Institute of Technology (the MITgcm, Adcroft et al. (2004)) with a plane-parallel, two-stream version of the multi-stream radiation code developed by Marley & McKay (1999). Further details on these models are provided in Showman et al. (2009)), Kataria et al. (2016) and references therein.

Using the prescription described in Fortney et al. (2006)) and Showman et al. (2009), we compute simulated eclipse depths derived from our GCM results for each planet over a period of 1-2 Earth years (Figure 6). At each Spitzer bandpass, the simulations, on average, exhibit a periodic variation in eclipse depth of 1-1.5 % over a period of $\sim$43 days for HD 189733b and a variation of 0.5-1 % over a period of $\sim$40 days for HD 209458b. This predicted variability is the result of a global sloshing mode or wave dynamics deep in the planetary atmosphere that leave only a small measurable effect at observable pressures Showman et al. (2009). Therefore, hot Jupiters such as these should exhibit low ($\sim$1%) variability over timescales much longer than the planet’s orbital period and smaller than the typical uncertainty in eclipse depth with current instruments.

The low observed atmospheric variability of HD 189733b and HD 209458b with Spitzer is in contrast to observations of brown dwarfs and solar system gas giants where temporal variability at infrared wavelengths has been observed at the 10-50% level (see review by Artigau, 2018). The low level of observable infrared variability for hot Jupiters like HD 189733b and HD 209458b compared with brown dwarfs and solar system gas giants likely arises from differences in the relative strengths in radiative and advective processes taking place in their atmospheres. Because brown dwarfs are self-luminous and comparatively fast rotators ($T_{orb}\sim$hrs), their circulation is dominated by multiple bands of jets and vortices (Showman et al. 2019). In contrast, the high stellar insolation and synchronous (and hence slower) rotation of hot Jupiters results in the emergence of strong day-night forcing that produce fast ($\sim$1 km/s) planetary-scale jets. The strong radiative forcing and global-scale circulation patterns in the observable portion of hot Jupiter atmospheres suppress small scale variations in the planet’s thermochemical structure that would contribute to large-amplitude variations in the dayside flux from the planet.

5.2. Further Observational Probes of Variability

Our ability to probe variability in exoplanet atmospheres is currently limited by the available targets and observational facilities. It is important to remember that facilities like Spitzer were not originally designed for high-precision time-series observations. With hot Jupiter targets that orbit bright nearby host stars like HD 209458b and HD 189733b; per eclipse precisions on the order of 100 ppm can be achieved with Spitzer. Future observational facilities, such as JWST , will both improve the achievable precision on eclipse depth measurements and greatly expand the wavelengths over which they can be obtained (Beichman et al., 2014). Hot Jupiter models predict significant variations in temperature, chemistry, and circulation patterns as a function of pressure level in the atmosphere (e.g. Showman et al., 2009), which translates into a strong wavelength dependence in the predicted levels of observable variability (e.g. Lewis et al., 2014). The observations presented here at 3.6 and 4.5 $\mu$m for HD 209458b and HD 189733b are predicted to probe a fairly limited pressure range in the $\sim$1-100 mbar level of these planets atmospheres (e.g. Showman et al., 2009). Future observations of hot Jupiters at higher precision spanning a larger range of wavelengths as well as observations of cooler exoplanets will be critical for expanding our understanding of the physical processes driving variability, or the lack thereof, in exoplanet atmospheres.

5.3. Non-Uniform Disc

There is prior evidence of non-uniform thermal distributions in both of the targets studied here. The analysis of the full orbit phase curves of HD 189733b (Knutson et al., 2012) provide evidence for an eastward shift of the hottest portion of the planetary atmosphere with respect to the substellar point of 5.29 $\pm$ 0.59 and 2.98 $\pm$ 0.82 hours in channels 1 and 2 respectively. These offsets correspond to longitudinal offsets of 35∘ $\pm$ 4∘ and 20∘ $\pm$ 6∘. Similarly, phase curve analysis of HD 209458b at 4.5 $\mu$m by Zellem et al. (2014) indicates an offset of 41∘ $\pm$ 6∘.

Williams et al. (2006) and Rauscher et al. (2007) predicted that a non-uniform dayside temperature distribution will change the shape of ingress/egress in comparison to a uniformly bright model. Also, as Williams et al. (2006) predicted and Agol et al. (2010a) observed, fitting the eclipse of a planet with a non-uniform dayside temperature distribution with a uniform occultation model will result in an apparent delay in the eclipse time as a result of these changes to the shape of ingress/egress when the timing of the eclipse is a free parameter. In Figures 2 and 3 we plot the difference between the center of eclipse time and what would be expected by assuming that the eclipse would occur at exactly one half of a period from the transit. HD 209458b consistently shows an offset in the eclipse time of 155 $\pm$ 25 seconds in Channel 1 and 129 $\pm$ 16 seconds in Channel 2. HD 189733b shows an offset of 28 $\pm$ 27 seconds in Channel 1 and 35 $\pm$ 13 seconds in Channel 2.

Figure 7 shows the magnitude of the effect of longitudinal offset on the observed center of eclipse time. We create a model planet with a hotspot of radius (0.5 Rplanet) utilizing the python package SPIDERMAN (Louden & Kreidberg, 2017). The global (out of hot spot) temperature is set to 1200 K and the hot spot temperature to 1700 K. Starting with the hot spot centered at the substellar point, we move the hotspot longitudinally in 1∘ increments and fit at each iteration with a uniform brightness model with center of eclipse time as a free parameter. We plot in Figure 7 the delay in eclipse time relative to the predicted time exactly one half period from transit. Note that some portion of apparent delay is present even in the 0∘ offset case. This is due to light travel time across the solar system ($\sim$ 45 sec. for HD 209). The offset of the maximum brightness observed in the phase curves would suggest that eclipse observations fit with a uniform model should exhibit an offset in eclipse timing of 40 and 60 seconds (Ch. 1 and 2) for HD 189733b and 90 seconds for HD 209458b channel 2 observations. The secondary eclipse fits presented here show strong agreement with the phase curve observations of HD 209458b. We are unable to confirm the phase curve offsets for HD 189733b with the precision attained with this set of observations.

As a further probe of non-uniformity, we refit the data while fixing the eclipse time to the expected value based on previous observations of epoch of transit and assuming a circular orbit. We combine the residuals of all the fits of each target in each channel to probe for any coherent deviations from the uniform model. Fig. 8 shows the stacked residuals from all of the channel 2 fits of HD 209. Here we see strong evidence for deviations in the shape of ingress/egress that could be caused by an eastward shifted hotspot. There are other factors that could also mimic this effect. The effect of eccentricity and shape of planet are explored in detail by de Wit et al. (2012). However, the precision with which the orbital properties of both of these canonical hot-Jupiters have been characterized allow many of factors to be ruled out or tight enough constraints to be placed on them to make a significant detection of non-uniformity likely in future eclipse mapping studies (Rauscher et al., 2018).

6. Conclusion

This analysis of multi - epoch secondary eclipse observations of the canonical hot Jupiters, HD 189733b and HD 209458b, finds no evidence for temporal variability greater than the precision of our observations. We expect variability in secondary eclipse depth to be driven by the dynamics of the atmosphere thus motivating comparison of observations to GCM predictions. The simulations, on average, exhibit a periodic variation in eclipse depth of 1-1.5 % for both planets at both 3.6 and 4.5 $\mu$m. Based on the uncertainty in each measurement of secondary eclipse depth, we are able to constrain variability to less than (5.6% and 6.0%) and (12% and 1.6%) for channels (1,2) of HD 189733b and HD 209458b respectively. The lack of evidence of variability provides motivation and justification to combine multi - epoch observations to achieve the precision necessary to perform high spatial resolution techniques such as eclipse mapping. The apparent offset of the center of eclipse time evident in these observations provides further evidence for a non-uniform dayside temperature distribution for the planet HD 209458b. We compare the evidence for non - uniformity to previous analyses of full phase curve observations and find agreement with HD 209458b observations at 4.5 µm. The precision achieved by these observations when phase folding and stacking all like observations is enough to resolve deviations to ingress/egress caused by a non - uniform temperature distribution in at least the case of HD 209458b at 4.5 $\mu$m. This level of precision meets the requirements to utilize the eclipse mapping technique to produce spatially resolved, two dimensional maps of the planetary dayside thermal distribution.

7. Acknowledgements

This work was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program under Grant Number 80NSSC17K0484. This work is based on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. We also acknowledge that part of this work was completed at the Space Telescope Science Institute (STScI) operated by AURA, Inc. This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org.