THE CLIMATE OF HD 189733b FROM FOURTEEN TRANSITS AND ECLIPSES MEASURED BY SPITZER

After conversion to counts, we flagged and cleaned the images of outliers, such as cosmic rays. We cleaned the images at the pixel level by rejecting outliers in the time series for each pixel. This worked well since the telescope pointed at nearly the same location for the entirety of each of our observations, so the flux in each pixel stayed relatively constant making outliers easy to flag. The pixel flux is affected by pointing variations, discussed in Section 3.2, so we did the outlier rejection by taking the difference between the pixel time series and a five-exposure (2 s) running median of the time series. The duration of the median was chosen to be shorter than the shortest of the pointing excursions, one of which occurred in the middle of the transit of the phase-function observation with a 1 pixel pointing change over 4 s.

For each pixel, the median-subtracted time series was sorted. Then, the standard deviation was computed from the 68.3% confidence limits on the median-subtracted time series. Finally, outliers were flagged in the median-subtracted time series that differed by more than 4σ from zero. The flagged pixels were then replaced by the five-exposure median.

This procedure performed well in eliminating cosmic rays and rogue pixels. After it was carried out, the photometry showed no significant outliers and, as we show below, was very close to the photon-noise limit.

Spitzer is affected by pointing variations that cause the star to change position slightly on the detector; this requires accurate centroiding to perform precise photometry. There are four varieties of pointing variations we observe in the data: small amplitude, short-timescale "jitter" which appears by eye as a damped random walk; a periodic pointing fluctuation which occurs on a timescale of ∼1 hr with an amplitude of about 0.1 pixel; gradual drifts over long timescales; and occasional short-timescale sharp excursions, as mentioned in Section 3.1.

In the initial stages of our data reduction, we found that the photometry was extremely sensitive to the pointing drifts. Since the 8 μm camera is undersampled, variations in pointing at the ∼0.1 pixel level can lead to ∼10% variations in the flux of the central pixels. For large photometric apertures, this is not a problem since small changes in the centroid do not change the enclosed flux much. However, large apertures contain pixels with lower illumination which have a longer-lasting detector ramp (see Section 4.2), so a smaller aperture containing higher illumination pixels is more desirable since these pixels have a ramp that saturates more quickly. Thus, we realized that a very accurate centroiding algorithm would be necessary for using smaller apertures, so we set out to test a wide range of different centroiding algorithms to see which performed best. We tried several centroiding algorithms: flux-weighted centroiding (e.g., Knutson et al. 2007), parabolic fitting (e.g., Todorov et al. 2010), point response function (PRF) fitting (e.g., Laughlin et al. 2009), and two-dimensional Gaussian fitting (e.g., Désert et al. 2009). The simplest algorithms to implement are the first two since these involve no optimization; the last two involve iterative nonlinear optimization, but the PRF fitting turned out to be too slow and problematic to implement.

We first tested the centroiding algorithms by creating simulated jitter using the IRAC Channel 4 PRFs. From these tests, we found that the two-dimensional Gaussian fitting gave the least scatter in the derived centroid relative to the input centroid; we found that keeping the x and y standard deviations the same gave as good a fit to the centroid as allowing the two to vary independently. We next ran the algorithms on the two stars (the target star and the M-dwarf companion; Bakos et al. 2006a) in the phase-function data from Knutson et al. (2007) and on the two stars in the observations of HD 80606 (Laughlin et al. 2009). These tests were critical since these were long time series so the stars had time to drift a significant fraction of a pixel, and the data contain noise. We first computed the centroid of both stars in each image (x₁, y₁ and x₂, y₂), then subtracted the x and y coordinates for each pair of stars (x₁ − x₂ and y₁ − y₂), and finally binned the x and y differences until the standard deviation of the x and y differences reached a minimum. A perfect centroiding algorithm ought to have perfect tracking between the two stars, resulting in a standard deviation only due to the photon shot noise and finite spatial resolution of the instrument. Various choices can be made for each of these algorithms, such as what portion of the array to fit or whether to smooth the data first before centroiding, so we spent some time experimenting with these and other choices.

In short, we found that the two-dimensional Gaussian performed the best out of all the centroiding algorithms. The algorithm selects a 7 × 7 sub-array from the image centered on the brightest pixel of a star. It then fits a two-dimensional Gaussian to this sub-array, allowing the center (centroid), amplitude, and width to vary, four free parameters in all for each image. We used the mpcurvefit.pro routine which implements a nonlinear Levenberg–Marquardt algorithm to optimize these parameters (Markwardt 2009). For the HD 189733 phase-function observations, the scatter in the data was 0.0018 pixels in x and 0.0051 pixels in y when binned by 512 exposures (205 s), while for HD 80606 the scatter was 0.0015 and 0.0021 in x and y when binned by four exposures (56 s); further binning resulted in minimal decrease of the scatter. The second best centroiding algorithm was the flux-weighted algorithm which had standard deviations of 0.016 and 0.032 in x and y for HD 189733, and 0.0019 and 0.011 in x and y for HD 80606. Thus, the two-dimensional Gaussian centroid performed better by a factor of ∼5 than the flux-weighted centroid; not only that, but the scatter in the flux-weighted centroid is due to a systematic error, while the scatter in the two-dimensional Gaussian centroid is almost completely random. This can be seen in Figure 1 which shows that the two-dimensional Gaussian centroid has weaker correlation of centroid difference of the two stars versus the centroid of one of the stars; on the other hand the other centroiding techniques show significant correlation between the offsets of the two stars and the pixel position of one of the stars, indicating a systematic error in the centroid determination.

Zoom In Zoom Out Reset image size

Applying these centroiding algorithms to our 12 transits and eclipses and the transit and eclipse from Knutson et al. (2007), we find that the scatter in the difference in centroids of the two stars (adding in quadrature the x and y components) ranges from 0.0035 to 0.0042 pixels for the Gaussian centroid after binning by 128 exposures (51 s); this is a factor of 3–7 times smaller than the flux-weighted centroid, and also no systematic trend in x, and a weak systematic trend in y. Since the Gaussian centroids of the two stars track one another well and the difference in their positions is nearly uncorrelated with their position on the detector, we are confident that the Gaussian centroid is giving the correct absolute position of these stars. When the data are fit with a 3.5 pixel radius aperture (as described in more detail below), the two-dimensional Gaussian centroid yields a χ² which is smaller for 9 of 14 transits/eclipses than the flux-weighted centroid, while the total χ² is smaller by 130 (after discarding the first 55 minutes of data for each eclipse/transit which has the steepest portion of the ramp).

In conclusion, we recommend the two-dimensional Gaussian for Spitzer IRAC Channel 4 sub-array centroiding of bright targets as it appears to behave in a near-optimal manner. As we will show, this results in very small scatter in the resulting photometry.

To subtract the background, we used a similar procedure as that described in Knutson et al. (2007); namely, we fit a Gaussian to a histogram of the counts from a subset of pixels located in the corners of the image (excluding the M-dwarf companion, and excluding the top row). This background contributes about 1.9%–2.6% of the total flux in our 4.5 pixel radius aperture, which we subtract from the time series frame by frame. As discussed by Harrington et al. (2007) and in the IRAC Instrument Handbook, the flux and the background of the 1st–5th and 58th frame of every set of 64 exposures is systematically lower than the other exposures in a set. However, after we carried out background subtraction frame by frame, the offset in these exposures does not appear in our total time series. Consequently, we believe it is due to a bias offset that affects the entire frame uniformly, and thus is easily removed.

Spitzer was not envisioned as an instrument for carrying out high-precision photometry on bright targets; consequently it was designed without sub-millimagnitude exoplanet photometry in mind. An instrumental artifact that appears at the ∼ mmag level in photometry, but can be up to 10% for low-illumination pixels over 33 hr, is the so-called detector ramp (Deming 2009): a pixel which is illuminated uniformly in time shows a gradual increase in the detected flux (see Figure 3). This is an important effect to correct for in fitting photometric time series; unfortunately there has not been a full understanding of this effect for the Spitzer IRAC detector. Here, we derive a toy model which qualitatively matches the behavior of the ramp (Section 4.2.1). We show that prior functions used for ramp corrections in other analyses of IRAC data (e.g., Deming et al. 2006; Désert et al. 2009) do not have the correct functional form to describe the observed ramp (Section 4.2.2).

Understanding how to remove the detector ramp has evolved with time. Early work (Deming et al. 2006) ratioed the transit star to other sources, and modeled the baseline in the ratio as linear or quadratic in time. Fitting functions to the ramp directly have used either exponentials in time (Harrington et al. 2007) or polynominals in the log of time (Knutson et al. 2009; Désert et al. 2009). These approaches have been adequate at lower signal-to-noise ratios (S/Ns), but for the present high S/N data we are motivated to find an improved functional form. We propose a new functional form, motivated by the toy model, which has the correct behavior and matches the observed ramps better (Section 4.2.3). We apply a range of tests to this ramp model, and show that it performs better than other ramp models (Section 4.2.4).

4.2.1. Toy Model for the Detector Ramp

The ramp effect is hypothesized to be due to trapping of electrons in detector defects ("charge trapping"). When a pixel is first illuminated, the charge traps are effectively empty, and some fraction of the electrons generated by the incident flux are retained by the traps instead of being read out by the array. As these charge traps fill, the effective gain of the detector goes up, until eventually the effect disappears. Thus bright, non-variable sources should have a detected flux that asymptotes to a constant value. When a pixel is not illuminated (or illuminated at very low intensity), the trapped charge gradually releases with time, causing the charge traps to become empty; this leads to ghost images after exposure of bright sources. A consequence of this model is that higher illumination pixels fill their charge traps more quickly, thus showing a much shorter detector ramp timescale. Although it is not clear that this model is correct, its predicted behavior agrees qualitatively with the observed IRAC photometric properties: for a bright source, the central pixels have a short ramp which saturates quickly, while the pixels with lower illumination in the wings of the point-spread function (PSF) show a more gradual ramp. This behavior, though, is difficult to model quantitatively as the pixel illumination varies with time due to Spitzer pointing variations (see Section 3.2).

A simple toy model can be developed for charge trapping as follows. Let γ_m be the fraction of volume of a detector pixel filled with charge traps, γ(t) be the fraction of volume of a pixel with empty charge traps at time t, and β be the total well depth (electrons). A pixel is illuminated below saturation with an intensity causing I(t) electrons to be released per second, while the measured intensity is I'(t) (e⁻ s⁻¹). As the pixel is illuminated, the charge traps fill up at a rate proportional to the intensity times the fraction of volume of empty charge traps; however, there is also a timescale τ at which electrons in charge traps are released, causing γ to increase. This gives

$$\begin{equation} {d\gamma \over dt} = -{I(t)\over \beta } \gamma -{\gamma _{\rm m}-\gamma (t)\over \tau }. \end{equation} \tag{ 1 }$$

The measured intensity is then

$$\begin{equation} I^\prime (t) = (1-\gamma) I(t) + \beta {\gamma _{\rm m}-\gamma \over \tau }. \end{equation} \tag{ 2 }$$

These equations have no closed-form solution for an arbitrary I(t); however, we can solve for their behavior in certain limits. For instance, if the illuminating intensity is constant, I(t) = I₀, for times t ⩾ t₀, then

$$\begin{equation} \gamma (t) = \gamma _{\rm m}(1-I\tau /\beta)^{-1} + \left(\gamma (t_{\rm 0})-\gamma _{\rm m}(1-I\tau /\beta)\right)e^{(\tau ^{-1}-I/\beta)(t-t_{\rm 0})}. \end{equation} \tag{ 3 }$$

As we are observing a bright star, we can simplify this equation by assuming that τ⁻¹ ≪ I₀/β, but we are still below saturation and in the linear regime (I₀t_exp ≲ 0.9β, where t_exp is the exposure time). In this limit

$$\begin{equation} I^\prime (t) \approx I_{\rm 0} (1-\gamma) = I_{\rm 0} (1-\gamma (t_{\rm 0}) e^{-{I_{\rm 0}\over \beta }(t-t_{\rm 0})}). \end{equation} \tag{ 4 }$$

This gives the expected ramp behavior: more strongly illuminated pixels have an apparent intensity, I', that asymptotes to a constant more quickly on a timescale β/I₀. At modest illumination, this becomes

$$\begin{equation} I^\prime (t) = I_{\rm 0} (1-\gamma _{\rm 0} + \gamma _{\rm 0} I_{\rm 0} \beta ^{-1} (t-t_{\rm 0})); \end{equation} \tag{ 5 }$$

a gradual linear ramp.

In the limit of zero illumination,

$$\begin{equation} \gamma (t) = \gamma _{\rm m}-(\gamma _{\rm m}-\gamma _{\rm 0}) e^{-(t-t_{\rm 0})/\tau }. \end{equation} \tag{ 6 }$$

Thus, the apparent intensity is

$$\begin{equation} I^\prime (t) = {\beta \over \tau } (\gamma _{\rm m}-\gamma _{\rm 0})e^{-(t-t_{\rm 0})/\tau }. \end{equation} \tag{ 7 }$$

This leads to persistent or ghost images that decay exponentially with time after observation of a bright target when the illumination is so strong that γ₀ ≪ γ_m.

4.2.2. Prior Models for the Detector Ramp

4.2.3. New Model for the Detector Ramp

Motivated by the toy model in Section 4.2.1, we decided to try an exponential ramp function. As this model predicts, the time constant is a function of pixel illumination. However, due to the pointing variations, we were not successful in correcting for the ramp on the pixel level. Instead we tried a ramp correction function that is simply the sum of two exponential terms:

$$\begin{equation} F^\prime /F= a_{\rm 0}-a_{\rm 1} e^{-t/\tau _{\rm 1}}-a_{\rm 2} e^{-t/\tau _{\rm 2}}, \end{equation} \tag{ 8 }$$

where F' is the flux affected by the ramp and F is the flux corrected for the ramp. Although this does not have exactly the correct behavior for the sum of pixels with different illuminations (assuming the toy model is correct), it does have the correct asymptotic behavior, and qualitatively represents the correction from higher and lower illumination pixels.

Figure 3 shows the ramp function overplotted with our data for the 14 transits and eclipses.

4.2.4. Performance of Double-exponential Ramp

We carried out photometry with apertures ranging in radius from 1 pixel to 7 pixels. We fit each transit and eclipse separately, and then computed the standard deviation of the data divided by the best-fit model (this is essentially the residuals in magnitudes). We find that the residuals are minimized at 4.39 mmag for an aperture radius of 3.5 pixels; this is the same aperture chosen in Knutson et al. (2007). For aperture radii of 3.0 pixels, the scatter is 4.5 mmag, while for 4.0 pixels it is 4.40 mmag, and for 4.5 pixels is 4.47 mmag, indicating that there is a shallow dependence on scatter with aperture size.

More importantly, we wish to minimize the presence of red noise in the residuals of the data. Consequently, we looked at the scatter in the residuals binned over a range of bin sizes from one exposure to 1920 exposures as a function of aperture size. We then took the mean of the scatter of the binned data over all 14 transits and eclipses, and computed the product of this mean scatter divided by the unbinned scatter over all bin sizes. The minimum occurs for an aperture of 4.5 pixels; although this has a slightly larger residual scatter without binning, the binned residuals are smaller relative to the unbinned residuals than for the 3.5 pixel radius aperture case. We have also measured the power spectrum of the residuals, and we find that the 4.5 pixel radius aperture minimizes the long-period power. Thus, we feel that this aperture size represents an appropriate compromise between small scatter in the unbinned data (which varies weakly with aperture size) versus minimization of the red noise component.

Figure 4 shows the scatter in the binned residuals, averaged over the 14 transits, as a function of bin size. Even up to bin sizes of 8832 exposures (1 hr bin), the scatter in the data does not deviate significantly from the inverse square root of the bin size; this indicates that the residuals are uncorrelated, and thus there is little (if any) red noise present. Remarkably the scatter in the 1 hr bins reaches 30 μmag; however, this is after subtracting off the double-exponential ramp model.

Zoom In Zoom Out Reset image size

For an aperture radius of 4.5 pixels, the median counts per exposure is 66,792. If photon counting errors dominate, then the expected noise level is 3.87 mmag per exposure. Including read noise (4.5 e⁻ per pixel) and sky noise (∼30 counts per pixel), the expected uncertainty is 3.99 mmag (we did not use the BCD uncertainties since these overpredict the noise properties according to the Spitzer Observer's Manual). The standard deviation in the residuals is 4.47 mmag, which is only 15.5% greater than the photon counting error and 12.0% greater than the expected errors including read noise and sky noise. Thus, the noise properties after correcting for the detector ramp are very close to the expected photon noise. The 3.5 pixel radius aperture has a residual scatter which is only 11% above the photon noise; however, this aperture size appeared to have more significant red noise, so we opted for the larger 4.5 pixel radius aperture.

We use the JPL Horizons ephemeris for the Spitzer orbit to convert the satellite time (keyword DATE_OBS) to BJD in Coordinated Universal Time (UTC).⁸ This correction is important since heliocentric and BJD can differ by as much as a few seconds, and different time systems can vary by seconds depending on the number of leap seconds included, which is close to the timing precision we can achieve with these data (Eastman et al. 2010).

We compute the errors on model parameters by calculating the residuals from each fit, shifting these by a random number of observations, and then adding the shifted residuals back to the best-fit model, then re-fitting; the so-called prayer-bead analysis (e.g., Agol & Steffen 2007). For each transit and eclipse, we carried out 2000 iterations of the prayer-bead shifts. This has the advantage of preserving correlations in the noise of the data that might still be present. For instance, if the ramp model is incorrect, there may be systematic deviation due to using the wrong ramp model, and these deviations are preserved within the residuals. This approach has some disadvantages; for instance, if the noise behaves differently within eclipse than outside eclipse, this might exaggerate the noise outside eclipse. Another disadvantage of this technique is that the number of independent trials is limited by the size of the data set since point order has to be preserved; consequently, we also randomly chose to reverse the residuals or change their sign to give more independent noise realizations. There is also the possibility that the effects of correlated noise may be removed in the fit. Even with these disadvantages we expect that this technique gives a fairly conservative estimate for the uncertainties on model parameters.

We allowed the model parameters to vary independently for each transit/eclipse. These fits were necessary since a simultaneous fit to the entire data set is computationally intensive due to the large number of data points; we avoided pre-binning the data to preserve as much information as possible about the noise in the final results. Figure 6 shows the average of all seven transits and all seven eclipses, corrected for the detector ramp and folded to the same orbital phase. We have binned the data by 8960 exposures to 69 data points for clarity.

Zoom In Zoom Out Reset image size

For the transits, we found that the sky velocity, v, and impact parameter, b, have no evidence for variation. Figure 7 shows each of these parameters plotted versus the transit number. The sky velocity has an average fractional uncertainty of 0.62%; the scatter in the measured values is 0.57%, and thus is consistent with being constant. Combining our data together, we find the average sky velocity is 25.125 ± 0.064 R_* day⁻¹, while a limit on the variation of the sky velocity is dv/dt = (−5.5 ± 6.6) × 10⁻⁴R_* day⁻². The impact parameter has an average measured value of 0.6631 ± 0.0023 R_*, with an average fractional uncertainty for each observation of 0.93% and a fractional scatter for the seven observations of 0.67%, also consistent with being constant. We constrain the change in impact parameter to be db/dt = (−0.02 ± 2.67) × 10⁻⁵R_* day⁻¹. Thus, our data indicate that the impact parameter and sky velocity of the transits remain constant to <1% over a duration of 590 days.

Zoom In Zoom Out Reset image size

We measured the transit and eclipse times for the seven transits and eclipses, shown in Figure 8, as well as in Tables 1 and 2. The errors on the transit times range from 2.4 to 3.6 s and are some of the most precise transit times ever measured, comparable to, or better than, the three HST transit times reported in Pont et al. (2007). We fit separate ephemerides to the transits and eclipses; the results are shown in Table 3. If we instead fit the transit times with the quadratic function: $t_{\rm n} = t_{\rm 0} + P n + \frac{1}{2} \dot{P} P n^2$, where t_n is the time of the nth transit, and $\dot{P}$ is the change in period of the orbit, we find $\dot{P} = -0.06 \,{\pm}\, 0.02$ s yr⁻¹. Since this is primarily due to the last data point, which may be an outlier, we do not view this as a significant detection.

Zoom In Zoom Out Reset image size

Table 1. Transit Parameters

Phase	F*	tc	$\sigma _{t_{\rm c}}$	Δtc	b	v	(Rp/R*)2	u1	−a1, − a2	τ1, τ2
	(counts)	(BJD-2454000 days)	(s)	(s)	(R*)	(R* day−1)	(%)		(counts)	(10−2 days)
−109.0	67222.6	37.611919 ± 0.000034	3.0	−4.2	0.6694 ± 0.0062	25.02 ± 0.15	2.4088 ± 0.0037	0.08 ± 0.03	75, 550	0.6010, 6.6616
1.0	66782.7	281.655291 ± 0.000040	3.4	−0.0	0.6586 ± 0.0066	25.28 ± 0.16	2.4022 ± 0.0047	0.11 ± 0.03	281, 107	0.4735, 5.4525
2.0	66722.4	283.873884 ± 0.000042	3.6	1.4	0.6592 ± 0.0066	25.24 ± 0.18	2.4253 ± 0.0063	0.11 ± 0.03	756, 752	0.4929, 5.5313
52.0	66758.8	394.802711 ± 0.000028	2.4	5.2	0.6636 ± 0.0067	25.05 ± 0.17	2.4333 ± 0.0051	0.13 ± 0.03	756, 712	0.4891, 5.6423
63.0	66995.3	419.207003 ± 0.000036	3.1	1.7	0.6594 ± 0.0049	25.18 ± 0.11	2.4225 ± 0.0049	0.12 ± 0.02	773, 722	0.5322, 6.2047
158.0	67385.6	629.971694 ± 0.000033	2.9	1.8	0.6641 ± 0.0062	24.90 ± 0.16	2.3984 ± 0.0062	0.16 ± 0.02	639, 570	0.4560, 5.5904
159.0	67242.2	632.190128 ± 0.000039	3.4	−10.4	0.6685 ± 0.0060	24.99 ± 0.16	2.3965 ± 0.0074	0.08 ± 0.03	715, 669	0.4890, 5.8697

Download table as:  ASCIITypeset image

Table 2. Eclipse Parameters

Phase	F*	tc	$\sigma _{t_{\rm c}}$	Δtc	Fp/F*	−a1, − a2	τ1, τ2
	(counts)	(BJD-2454000 days)	(s)	(s)	(%)	(counts)	(10−2 days)
−108.5	67466.9	38.722278 ± 0.000265	22.9	9.1	0.3345 ± 0.0057	0, 0	0.0000, 0.0000
0.5	66870.3	280.546423 ± 0.000263	22.7	−32.5	0.3469 ± 0.0060	609, 535	0.3793, 4.5632
1.5	66808.5	282.765713 ± 0.000350	30.2	29.3	0.3420 ± 0.0068	255, 118	0.1134, 4.8741
51.5	66904.4	393.693798 ± 0.000414	35.8	−26.2	0.3368 ± 0.0042	750, 727	0.6350, 6.1896
63.5	67097.2	420.317184 ± 0.000287	24.8	16.2	0.3623 ± 0.0064	759, 647	0.5089, 6.0830
157.5	67582.4	628.862649 ± 0.000390	33.7	−30.8	0.3378 ± 0.0091	330, 239	0.0310, 5.6520
158.5	67318.3	631.081875 ± 0.000313	27.0	25.5	0.3477 ± 0.0075	737, 685	0.5505, 5.6994

Download table as:  ASCIITypeset image

The uncertainty on the transit times and eclipse times, as well as the derived ephemerides, is inversely proportional to depth of the transits and eclipses. This is due to the fact that the timing precision is proportional to the inverse of the flux gradient with time during ingress and egress. The ingress and egress durations are the same for the transit and eclipse (assuming a circular orbit), so the ratio of the flux gradient scales with the ratio of their depths. The ratio of the depths is proportional to the ratio of the surface brightness of the star to the surface brightness of the planet (limb darkening is weak for this star at 8 μm), so the transit time precisions are smaller than the eclipse precisions by the ratio of the surface brightness of the planet to the star, which is about 14.3%, or a factor of 7.0.

Figure 8 shows the deviations from the ephemeris for both the transits and the eclipses. The transits have a scatter of 5.1 s, which is very close the observational errors; there is no evidence in our data for TTVs over a period of 590 days. The eclipses also appear to be precisely periodic—their scatter with respect to the best-fit ephemeris is 27 s, which is comparable to the errors on each data point. The period derived from the transits differs from that derived from the eclipses by only 0.1 s!

The eclipses appear 69 ± 11 s later than 1/2 of an orbital period after the transits. As discussed in Knutson et al. (2007), this is partly due to the light travel across the system, 30.8 ± 0.6 s (this uncertainty is due to the uncertainty in stellar mass, M_* = 0.806 ± 0.048 M_☉; Torres et al. 2008), while the remaining 38 ± 11 s can be mostly accounted for by the hot spot on the planet causing an offset in the time of eclipse when the planet is modeled as a uniform disk, as shown in Section 6.5.

These transit data show no significant timing variations, but from these we can constrain the maximum mass allowed of additional planets in the system. Transit timings are a particularly sensitive probe for planets in or near mean-motion resonance (MMR) and previous studies have ruled out Earth mass or super-Earth mass planets in low-order MMR for several systems. Prior TTV analyses of the HD189733 system (Hrudková et al. 2010; Miller-Ricci et al. 2008) used data with timing precision of order 30 s and were sensitive to planetary masses of near (and below) 1 Earth mass in favorable MMRs. Our Spitzer observations of HD189733 have nearly a factor of 10 better timing precision and consequently have improved sensitivity to secondary planets by that same factor. Here, we calculate the maximum mass that an additional planet could have based upon these transit data. To do so we note that the χ² per degree of freedom of the timing residuals is slightly more than unity. We therefore scale the timing uncertainty by a factor of 1.5 and then multiply by two to achieve our 2σ (≃ 95% confidence level) upper bound on the timing variations.

Figure 9 shows 95% confidence-level constraints on secondary planets with near circular orbits in this system based upon these data and the RV measurements from Boisse et al. (2009). These limits are derived from the analytic formula given in Agol et al. (2005) and Steffen & Agol (2005). We do not attempt an in-depth numerical analysis of these transit times here—the robustness of limits derived from analytic formulae was demonstrated in Agol & Steffen (2007), Nesvorný (2009), and Nesvorný & Beaugé (2010). These data exclude planets above 2 Earth masses for any orbit that lies closer to the known planet than either the interior or exterior 2:1 MMR. The transit-timing mass exclusion is superior than the exclusion from RV data for periods from 1 to 5 days, excluding all planets with masses greater than 3 Earth masses within this range. In addition they exclude planets with masses well below the mass of Mars—approximately 0.2 Mars masses or 2 Moon masses at 95% percent confidence—in circularly orbiting 2:1 or 3:2 MMRs (interior or exterior). For non-circular orbits, the sensitivity generally increases. However, in low-order MMR the mass sensitivity can decline as much as a factor of 10 for eccentric orbits (see, for example, Agol & Steffen 2007).

Zoom In Zoom Out Reset image size

As discussed in Knutson et al. (2007), the 8 μm phase function indicates that the hottest point on the planet is offset from the sub-stellar point. This was predicted by Cooper & Showman (2005), and is attributed to the advection of energy by a super-rotating wind encircling the equator of the planet. This offset hot spot means that the ingress and egress of the secondary eclipse will have a shape that differs from our model, which utilizes a uniform disk. In particular, this means that the steepest portion of ingress and egress will be offset from the uniform disk case; since the hot spot is on the trailing side of the planet with respect to the direction of motion, this causes a delay in the eclipse time when fit with a uniform disk model (Charbonneau et al. 2005; Williams et al. 2006). In Knutson et al. (2007), we estimated that the hot spot would cause a delay of at most 20 s; however, the fit to the phase function in that paper did not correct for stellar variation which caused the location of the hot spot to be underestimated, leading to an underestimated uniform time offset.

To estimate the magnitude of this effect, we used a simplified model of the longitudinal planet brightness which is discussed in Cowan & Agol (2010b). Briefly, each position on the planet is treated as a parcel of gas which moves eastward, absorbing star light as it passes across the day side, all the while radiating with a time constant τ_rad. This model can be parameterized by a single parameter, = τ_rad/τ_adv, where τ_adv is the time it takes a parcel of gas to circle the planet. Small values of ("instant" re-radiation) lead to darker night sides and dayside temperatures which are in equilibrium with the incident stellar flux. Large values of lead to nearly uniform temperatures at each latitude. Thus, in the small or large limits we expect no timing offset since the day sides are symmetric. Only with ∼ 1 is there an offset hot spot causing a phase function which peaks before secondary eclipse, as well as a slight offset in the times of eclipse ingress and egress if fit with a uniform planet.

We computed the effect of on the time of secondary eclipse by solving for the planet dayside longitudinal surface brightness at the equator in the Rayleigh–Jeans limit and assuming a constant temperature with latitude. We computed the eclipse ingress and egress from this model for the planet surface brightness, we fit this simulated eclipse light curve with a model for the eclipse of a uniform planet, and from this best fit we determined the offset in the time of eclipse, the so-called uniform time offset defined by Williams et al. (2006). Figure 10 shows this time offset as a function ; a positive offset means that the secondary eclipse occurs later than expected for a uniform planet. The maximum offset predicted by this model is 43 s, which agrees with the observed eclipse time offset. Our measured eclipse time is plotted as a dashed line in this figure, with the uncertainty indicated by the horizontal shaded rectangle.

Zoom In Zoom Out Reset image size

The location of the peak in the planet phase function from Knutson et al. (2007) provides another constraint on the location of this hot spot, or equivalently on the value of (see Figure 11). We used the relation between the infrared and optical stellar variability derived in Section 5 to derive the change in stellar flux at 8 μm during the phase-function measurement, about 0.6 mmag over 26.6 hr. We then fit the last 2/3 of the measured 8 μm phase function to estimate (Figure 11); we discarded the first 1/3 of the phase-function data since it is strongly affected by the ramp correction. We find a best-fit value of = 0.74 ± 0.07, which we have also plotted as a vertical shaded region in Figure 10. This value of predicts a timing offset of 33 s, which is consistent with the measured offset of 38 ± 11 s. Figure 12 shows a direct comparison of the secondary eclipse to the average of our seven secondary eclipses. The top panel shows the binned data as well as the best-fit secondary eclipse at 1/2 orbital period after transit plus the 30.8 s light-travel time delay (solid line), as well as the = 0.74 model with light-travel time delay (dashed line). The bottom panel shows the residuals binned into eleven bins: pre-ingress, post-egress, eclipse, and four bins each in ingress and egress; the residuals are plotted for the uniform planet model (diamonds) and = 0.74 model (filled circles with error bars). The uniform planet model shows points which are on average higher in ingress and lower in egress, which is a sign of the shifted hot spot. The hot-spot model provides a better fit to the data, although there is still scatter in the residuals which just reflects the low significance of the eclipse hot-spot detection (the uniform time offset is only 3.5σ: 38 ± 11 s). Note that we have not optimized the hot-spot model, but only computed the light curve from the best fit to the phase function (which gives = 0.74).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Consequently, there is no evidence for non-zero ecos ω in this system. For the estimated value of , the remaining time offset is 6 ± 11 s, which yields ecos ω = 0.00005 ± 0.00009.

Another prediction of this model is the day–night brightness difference. For = 0.74, we predict a nightside brightness which is 57% of the dayside brightness (day defined at mid-eclipse; night at mid-transit). This corresponds to a decrease in brightness from the day to night side which is 0.15% of the stellar flux, which is very close to the value of 1.2 ± 0.2 mmag derived in Section 5. The minimum planet brightness (for the visible hemisphere) divided by the maximum planet brightness for this best-fit model is 50%, while the peak in observed planet brightness is 23°, 0.065 orbital periods, or 3.5 hr before the secondary eclipse. On the planet, the hottest longitude is 13° from the sub-stellar point; note that this differs from the hemispherically integrated peak due to asymmetry in the longitudinal dayside intensity. Figure 11 shows the phase-function data after correction for stellar variation and binned by 1000 data points, versus planet orbital phase with the best-fit model for the planet variability, = 0.74, overplotted. Also plotted is our estimate of the nightside flux based on Equation (9).

In the fits to the transits, we allowed the depth of each transit to vary independently through the ratio of the planet to stellar radius, p = R_p/R_*. This ratio also affects the duration of ingress and egress, but the amount of time spent in ingress or egress is small, so the primary effect is on the depth of transit. Figure 13 shows the derived transit depths, p², for the seven observed transits. There is evidence for variability in the transit depth—the uncertainty on the individual transit depths ranges from 37 to 74 μmag, while the scatter is 145 μmag. This corresponds to a scatter in the fractional variation of transit depth of 0.6%, while the ratio of maximum to minimum is 1.5%. Fitting the transit depths with a single value of p gives a Δχ² of 40.8 for 6 degrees of freedom. Thus, transit depth variation is detected with high significance.

Zoom In Zoom Out Reset image size

We allowed the limb-darkening parameter to vary for each transit, which ranged from 0.08 to 0.13 for the seven transits, with an overall mean of u₁ = 0.12 ± 0.01. Even though limb darkening is weaker in the infrared, an LTE Kurucz model atmosphere (Kurucz 1992) with parameters close to the values inferred for HD 189733, T_eff = 5000 K, $[\rm Fe/H] = 0.0$, and $\log [g\rm {(cm \,s^{-2})}] = 4.5$, predicts a linear limb-darkening coefficient of 0.136, close to what we measure. We checked that the variation in transit depth is not due to variations in the best-fit limb-darkening parameter by holding the limb-darkening coefficient fixed at the mean value; this did not affect our measured values of transit depth.

In fact, the variation in transit depth is not necessarily due to a change in p. The average depth of each transit is given by δ = 〈I_path〉πR²_p/(F_* + F_p) (Mandel & Agol 2002), where 〈I_path〉 is the average surface brightness within the path of the planet across the star and F_p, F_* are the planet and stellar flux during transit. Although our model assumes a linear limb-darkening law, if the path of the planet passes over a region of the star with brighter than average surface brightness, then a larger depth will be inferred.

So, it is possible that the change in transit depth is due to changes in 〈I_path〉/F_*, R_p, R_*, or F_p/F_*. Variations in R_p or R_* seem unlikely to be responsible for the transit depth variation as this would require changes in radius of 0.3% for either body: either 220 km for the planet (∼1/3 of a thermal scale height) or 1600 km for the star. Both of these variations seem too large to be physically plausible. Variations in 〈I_path〉/F_* due to fluctuations in F_* are ruled out as the transit depth variations are uncorrelated with the optical stellar flux variations. Variations in F_p (the nightside planet flux) are less than 0.35 mmag relative to F_*, which is too small by a factor of 17 to be responsible for the transit depth variations. Thus, the most likely possibility is that the transit depth variations are due to a variation in the occulted stellar intensity, 〈I_path〉. This requires only a variation of 0.6% in the surface brightness of the path of the planet relative to the average stellar surface brightness, which is much smaller than the 12% change in surface brightness from center to limb inferred for the best-fit limb darkening. We have checked that individual star spots are not responsible for this variation by computing the standard deviation of the residuals in transit divided by the square root of the counts, which is 1.153, versus the same quantity computed out of transit, which is 1.154, so there is no evidence for individual star spots causing this difference. Similar inferences have been drawn for Channels 1 and 3 by Beaulieu et al. (2008), while star spots are easier to detect in the optical where the contrast between star spots and the stellar disk is larger, and have already been detected for this star by Pont et al. (2007), who also resolved the shape and color dependence of the spots.

Figure 14 shows the measured eclipse depth in units of the stellar flux at mid-eclipse. The weighted mean eclipse depth is 0.344% ± 0.0036% and the χ² fit to the eclipse depths with this mean value is 14.6 for 6 degrees of freedom. The errors on each eclipse depth vary between 0.004% and 0.009%, while the scatter in the depths is 0.01%, so there is no detection of significant eclipse depth variability. This scatter corresponds to 2.7% variation of the mid-eclipse planet brightness; this can be taken as an upper limit on the planet variability at 68% confidence.

Zoom In Zoom Out Reset image size

Although our eclipse depths do not exhibit significant variability, they only probe variability on timescales shorter than the baseline of the observations (roughly 2 years). To verify that our data do not show any time structure, we plot in Figure 15 the change in eclipse depth against the time between observations, for each pair of observations (N × (N − 1)/2 = 21 for N = 7). We add uncertainties in quadrature to estimate the uncertainty on the flux differences. The resulting locus is flat, showing that measured eclipse depth is uncorrelated with the time of the observation. Note that this sets a limit not only astrophysical variability scenarios, but also systematic errors.

Zoom In Zoom Out Reset image size

This limit on the variation in eclipse depths is sufficient to rule out the predicted variation computed for HD 189733b by Rauscher et al. (2008). The most extreme prediction they make is for their η = 0.05, ${\bar{U}} =$ 800 m s⁻¹ model which has a standard deviation of ∼8% in the dayside brightness, with the largest excursions of 20%. The largest difference in brightness we see is 8%, while the scatter in planet brightness is 2.7%, so by both measures the observed variation is a factor of ∼3 smaller than the predictions of this particular model.

Other models predict smaller variations in the dayside brightness, such as Showman et al. (2009b) who compute the 8 μm brightness variation for HD 189733b should be less than 1%. Our upper limits are consistent with this model, but unfortunately do not constrain the model due to our uncertainties that are larger than the predicted variations.

With upper limits on both day and nightside variability, it is worth asking which of these puts stronger constraints on the planet's physical properties. Consider the simple model of Cowan & Agol (2010a), which parameterizes the planet's day and nightside brightnesses in terms of the planet's Bond albedo and recirculation efficiency. One can treat brightness variability—of the day or night—as being due to changes in albedo and/or changes in recirculation efficiency, and compute how these affect the day and nightside brightness. We find that the dayside variability upper limit provides a better constraint on variation in the Bond albedo or recirculation efficiency than does the nightside variability limit, which is ∼6× larger than the dayside limit.

Another possible origin of planet dayside variability are transient local variations in the surface brightness, for example, due to large-scale "storms." We use a toy model where the planet's day side has a uniform temperature, T_d, except for a storm with covering fraction 0 < f < 1 (the y-axis) and temperature difference ΔT from mean dayside temperature (the x-axis). Figure 16 shows exclusion limits on the largest putative storm that could form or dissipate without appearing in our data. The increasingly dark shades of gray denote areas of parameter space excluded at 1σ through 12σ. According to Showman et al. (2009a), the radius of deformation for HD 189733b is 0.3 of the planetary radius, an order of magnitude larger than for Jupiter. The covering fraction for a typical storm on such a planet would be f = 0.1, for which we cannot rule out storms differing by 324 K (68% confidence) from the average dayside temperature. For comparison, Jupiter's Great Red Spot has a filling fraction of roughly 0.03, and a temperature 12 K cooler than the rest of the planet. The bottom line is that our data rule out only the most extreme weather fluctuations on HD 189733b.

Zoom In Zoom Out Reset image size

Due to the high precision of our data and weak limb darkening in the infrared, we can considerably improve the determination of certain stellar and planet parameters for this system from our data. Since both the small inferred value of ecos ω and theoretical predictions indicate that e should be close to zero, we set e = 0 in deriving the system parameters. The uncertainties on the stellar and planet parameters are computed for each transit or eclipse by computing the system parameters from the model parameters from each simulation (using the relations in Winn 2010), computing the standard deviation of the results from the simulations as an estimate of the errors on each parameter, and then taking a weighted mean of all transits/eclipses to obtain the final mean value of the best-fit parameters.

Table 4 presents the system parameters determined from all 14 transits and eclipses. We have focused on parameters that are most directly constrained from the photometry, which are either dimensionless, or have units of density. Compared to the values derived in Torres et al. (2008) and Pont et al. (2007), the uncertainties on our values are smaller by a factor of 2–10. For the planet surface gravity, g_p, we use the velocity semi-amplitude K = 200.56 ± 0.88 m s⁻¹, derived by Boisse et al. (2009), and for the planet density, ρ_p, we use the stellar mass M_* = 0.806 ± 0.048 M_☉ given in Torres et al. (2008), propagating the uncertainties assuming they are uncorrelated and Gaussian.