Free Access
Issue
A&A
Volume 526, February 2011
Article Number A12
Number of page(s) 13
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/200913093
Published online 14 December 2010

© ESO, 2010

1. Introduction

The ransiting hot jupiter HD 189733b orbits a small, bright, main-sequence K2V star (K = 5.5) and produces deep transits of  ~2.5% (Bouchy et al. 2005). The planet has a mass of Mp = 1.13 Jupiter mass (MJ), a radius of Rp = 1.16 Jupiter radius (RJ) in the visible (Bakos et al. 2006b; Winn et al. 2007), and a brightness temperature ranging between 950 and 1220 K (Knutson et al. 2007a,b). This implies a large atmospheric scale height (200 km), allowing transits to characterize the chemical composition and vertical structure of the atmosphere.

Transmission spectroscopy during transit allows the upper part of the atmosphere to be probed down to altitudes where it becomes optically thick. Therefore, one can infer the composition of the atmosphere by measuring the wavelength-dependent planetary radius (Seager & Sasselov 2000; Burrows et al. 2001; Hubbard et al. 2001). Pioneering observational work using this method on HD 209458b with the Hubble Space Telescope (HST) has led to the first detection of an exoplanetary atmosphere (Charbonneau et al. 2002) and its escaping exosphere (Vidal-Madjar et al. 2003, 2004, 2008; Ballester et al. 2007). The optical transmission spectra of this planet shows evidence of several different layers of Na (Sing et al. 2008a,b) as well as Rayleigh scattering by molecular hydrogen (Lecavelier des Etangs et al. 2008b) and the possible presence of TiO/VO (Désert et al. 2008). Sodium has been detected in the atmosphere of HD 189733b with ground-based observations (Redfield et al. 2008). Using the Advanced Camera for Survey (ACS) onboard HST, Pont et al. (2007, 2008) detected atmospheric haze at high altitude, interpreted as Rayleigh scattering, possibly by small particles (Lecavelier des Etangs et al. 2008a). The detection of water and methane reported by Swain et al. (2008) in spectroscopic observations between 1.5 and 2.5 μm has been challenged by Sing et al. (2009) with new HST observations at 1.66 and 1.8 μm. More recently, Lecavelier des Etangs et al. (2010) reported the presence of an extended exosphere of atomic hydrogen surrounding the planet using far-UV observations.

The SpitzerSpace Telescope (Werner et al. 2004) has demonstrated its potential to probe exoplanetary atmospheres through emission and transmission spectra. Multi-wavelength observations of the eclipses (secondary transits) of HD 189733b behind the star have made possible the study of the atmosphere’s emission spectrum (Charbonneau et al. 2005; Deming et al. 2005, 2006, 2007). This planet has been extensively observed during secondary transits and throughout the phases of its orbit (Knutson et al. 2007a,b, 2009). Strong water (H2O) absorption in the dayside emission spectrum has been detected along with possible variability (Grillmair et al. 2007, 2008). The likely presence of carbon monoxide (CO, Charbonneau et al. 2008), and carbon dioxide (CO2, Swain et al. 2009) have also been reported.

Knutson et al. (2007a,b) obtained the first accurate near infrared (NIR) transit measurements for this planet at 8.0 μm using the Infrared Array Camera (IRAC, Fazio et al. 2004) onboard Spitzer. The search for molecular spectroscopic signatures by comparing two photometric bands with Spitzer has been attempted by Ehrenreich et al. (2007). This study concluded that uncertainties in the measurements were too large to draw firm conclusions on the detection of water at high altitudes. Yet, using the same data set but a different analysis (Beaulieu et al. 2008), Tinetti et al. (2007) claimed to detect the presence of atmospheric water vapour from comparison of the absorption at 3.6 μm and 5.8 μm.

We recently performed a consistent and complete study of these observations obtained with the four IRAC channels (Désert et al. 2009), including a detailed assessment of systematics. We concluded that there is no excess absorption at 5.8 μm compared to 3.6 μm. We also derived a slightly larger radius at 4.5 μm that cannot be explained by H2O or Rayleigh scattering. We interpreted this small absorption excess as caused by to the presence of CO molecules (Désert et al. 2009). This is consistent with the low level of emission measured from the planetary eclipse at the same wavelength (Charbonneau et al. 2008).

Here we present new Spitzer/IRAC observations of a primary transit of HD 189733b at 3.6 μm. This work is part of our ongoing efforts to characterize the transmission spectrum of HD 189733 using space-based observatories (Désert et al. 2009; Sing et al. 2009). The present work focuses on new Spitzer/IRAC observations of the planetary transit, where the photons were accumulated in subarray mode. The data presented here consist of a high-cadence time series, which is described in Sect. 2 and analyzed in Sect. 3. We present a method to correct for the stellar variability in Sect. 4. The photometric precision achieved with this observation allows us to derive accurate measurements of the planetary radius and orbital parameters. Our results together with a comparison with previous studies and theoretical predictions are also given in Sect. 4.

2. Observations and data reduction

2.1. Observations

We obtained Spitzer General Observer’s time in Cycle 3 and 4 (PI: Vidal-Madjar; program IDs 30590 and 40732), in the four Spitzer/IRAC bandpasses. In total, three primary transits of HD 189733b were observed. Our primary scientific objectives were to detect the main gaseous constituents (H2O and CO) of the atmosphere of this hot-Jupiter.

Our first observations of the HD 189733 system (hereafter visit 1) were acquired simultaneously at 3.6 and 5.8 μm (channels 1 and 3) on 2006 October 31. Visit 2 was completed 2007 November 23 with measurements at 4.5 and 8 μm (channels 2 and 4). For both visits, the system was observed in IRAC’s stellar mode for 4.5 h, of which 1.8 h was during planetary transit. The observations were binned into consecutive subexposures with a cadence of 0.4 s for channels 1 and 2 and 2.0 s for channels 3 and 4. We obtained a total of 1936 frames for channels 1 and 3 and 1920 frames for channels 2 and 4. Details of the analysis of these data (Visit 1 and 2) can be found in Désert et al. (2009).

Here we focus on the data obtained during Visit 3 on 2007 November 25 (cycle 4, program 40732) at 3.6 μm (channel 1); its analysis is presented in this paper. We used IRAC’s 32  ×  32 pixel subarray mode (1.2′′ × 1.2′′ pixels) for visit 3, which lasted 4.5 h with 1.8 h during transit. The primary goal of this last visit was to obtain an accurate measurement of the planetary radius at 3.6 μm. The subarray mode has a greater efficiency than the stellar mode in collecting photons.

HD 189733 is brighter in channels 1 and 2 than in channels 3 and 4. For visit 1, the flux was 1700 mJy at 3.6 μm, which is close to the saturation limit with 0.4 s exposure time as used in visit 1. IRAC’s subarray mode can solve this issue by observing at high speed cadence of 0.1 s with sets of 64 subframes taken back-to-back with no gaps between subframes. During visit 3, we obtained 1920 consecutive exposures, yielding 122 880 frames with an effective exposure time of 0.08 s and a frame interval of 0.1 s. A delay of 8.3 s occurs between consecutive exposures. While planning the observation in the Astronomical Observing Request (AOR) format, we carefully avoided placing the star near dead pixels. We also did not dither the pointing so the target remained on a given pixel throughout the observation. This minimizes errors from imperfect flat-field corrections and therefore increases the relative photometric accuracy. This strategy has been used successfully for several previous Spitzer observing runs on HD 189733 (Knutson et al. 2007a,b; Ehrenreich et al. 2007; Agol et al. 2009; Désert et al. 2009), on HD 149026 (Nutzman et al. 2009), on GJ 436 (Deming et al. 2007; Gillon et al. 2007; Demory et al. 2007), and more recently with Warm-Spitzer on HD 80606 (Hébrard et al. 2010).

2.2. Data reduction

For visit 3, we used basic calibrated data (BCD) frames, produced by the standard IRAC calibration pipeline, to correct for dark current, flat-fielding, and detector nonlinearity and to convert the observations into flux units. We are able to measure the centroid position of a stellar image (computed by Gaussian fitting) to a precision of 0.01 pixel, using the DAOPHOT-type Photometry Procedures, GCNTRD, from the IDL Astronomy Library1. We find that the centroid position varies by less than 10% of a pixel during a complete observation. We used the APER routine to perform aperture photometry with a radius of 5 pixels to minimize the contribution of HD 189733B, which stands at an angular separation of (Bakos et al. 2006a).

The background level for each frame was measured with APER as the median value of the pixels inside an annulus centered on the star with inner and outer radii of 16 and 18 pixels, respectively. Typical background values are 10 ± 0.2 electrons per pixel compared to  ~110 000 electrons of total stellar signal. Therefore, photometric errors are not dominated by fluctuations in the background. Short exposures yield a typical signal-to-noise ratio (S/N) of  ~300 per individual observation.

After producing the photometric time series, we iteratively selected and trimmed outliers greater than 4σ by comparing the 122  880 photometric measurements to the best-fit transit light curve model. This removed any remaining measurements affected by transient hot pixels. In this way, we discarded 138 frames, or approximately 0.1% of the observations. We also discarded the first 2542 measurements (planetary phase before  −0.0472), which are affected by a small ramp effect. Thus, our cleaned light curve consists of 120 200 photometric measurements. We then binned the transit light curve by a factor of four to increase computing efficiency without losing information on the pixel-phase effect (see below). Our final cleaned, binned transit light curve contains 30 050 data points.

3. Analysis

3.1. Fitting the transit light curve

Owing to instrumental effects, our measured out-of-transit flux values are not constant. Telescope pointing jitter results in fluctuations of the stellar centroid position (Fig. 1), which, in combination with intra-pixel sensitivity variations, produces systematic noise in the raw light curves (Fig. 2, upper panel). A description of this effect, known as the pixel-phase effect, is given in the Spitzer/IRAC data handbook (Reach et al. 2006, p. 50; see also Charbonneau et al. 2005). To correct the light curve, we define a baseline function that is the sum of a linear function of time and a quadratic function of the X and Y centroid positions. This function, with five parameters, Ki, is described in Désert et al. (2009).

We modeled the transit light curve with four parameters: the planet-star radius ratio, Rp/R, the orbital semi-major axis to stellar radius ratio (system scale), a/R, the impact parameter, b, and the time of mid transit, Tc. We used the IDL transit routine OCCULTNL, developed by Mandel & Agol (2002), to model the light curve. Our limb-darkening corrections consist of three non linear, limb-darkening coefficients, as defined in Sect. 4.3 and presented in Table 1. We also used a linear function of time as the baseline (Aj, j = 1, 2) as presented in Désert et al. (2009).

We used the MPFIT package2 to perform a Levenberg-Marquardt least-squares fit of the transit model to our cleaned and binned observations. The best-fit model was computed over the whole parameter space (Rp/R, a/Rb, Tc, Aj, Ki). The baseline function described above is combined with the transit light curve function so the fit is constrained by eleven parameters (four for the transit model, two for the linear baseline, and five for the pixel phase effect); this gives 30  039 degrees of freedom.

3.2. Mean values and errors determination

We used the prayer bead method (Moutou et al. 2004; Gillon et al. 2007) to determine the mean value as well as the statistical and systematical errors for the measured parameters.

As shown by Pont et al. (2006), the existence of low-frequency correlated noise (dubbed red noise) between different exposures must be considered to obtain a realistic estimation of the uncertainties. To obtain an estimate of the systematic errors in our observations we followed the method described by Gillon et al. (2006) and derived the covariance from the residuals of the light curve. In the prayer bead method, the residuals of the initial fit are shifted systematically and sequentially by one frame, and then added to the transit light curve model before fitting again. The error on each photometric point is the same and is set to the rms of the residuals of the first fit obtained.

Totally, 30  050 shifts and fits of transit light curves were produced to derive a set of parameters and to extract their medians and their corresponding errors. We set our uncertainties equal to the range of values containing 68% of the points in the distribution in a symmetric range about the median for a given parameter (Fig. 3). Error bars estimated using the prayer bead method are generally larger than the ones obtained using other methods like “delta khi square” (Fig. 4) or Markov Chain Monte Carlo (MCMC) because they include the effect of the residual red noise. Therefore these error bars are likely the most conservative.

thumbnail Fig. 1

X (black) and Y (red) position of the star centroid in pixels as function of time (binned by 4). These positions are used in the decorrelation function that is used for the fit of the transit light curve.

Open with DEXTER

4. Results and discussion

thumbnail Fig. 2

Top panel: raw transit light curve obtained in subarray mode binned by four with its fit overplotted in red. Pixel-phase effect is clearly visible. Middle panel: normalized transit light curve. Bottom panel: residuals obtained from the fits of the transit light curve.

Open with DEXTER

The raw aperture photometry of the transit light curve obtained in subarray mode binned by four is plotted in Fig. 2 together with the lightcurve decorrelated from the dependence on instrumental parameters. We derived the system parameters and their error bars in a consistent way using the distribution obtained from the prayer bead method and present the results in Table 2. From our fits, we evaluated the radius ratio Rp/R, the impact parameter b, the system scale a/R and the central time of the transit Tc (Table 2). We obtain a mean χ2 of 30  024 for n = 30 039 degrees of freedom. However, we found that the parameters distributions are far from being Gaussian (Fig. 3). This is explained by the presence of a non negligible red-noise left after corrections as revealed by comparing the rms of the unbinned and binned residuals (Fig. 5) following the method described by Gillon et al. (2006). This strenghtens the suggestion that the error bars would have been smaller if they had been estimated with methods other than the prayer bead.

4.1. Signal-to-noise ratios

The standard deviation of the difference between the best-fit model and the observed transit light curve is found to be 1.8 × 10-3 (rms) on individual data points. This value corresponds to a signal-to-noise of 560:1 for a frame binned by four; as expected, this is nearly twice the theoretical (photon noise) signal-to-noise ratio for individual frames. We achieve an average of 91% of the theoretical signal-to-noise on individual frames.

Table 1

Limb-darkening coefficients used for IRAC/channel 1 bandpass centered at 3.6 μm.

We also estimated the red-noise by comparing the rms of residuals with different bin sizes (Fig. 5) and we find that the correlated systematics left after decorrelation of the pixel-phase effect have an amplitude of 1.4 × 10-4 in a bin of a 1000 points. The amplitude of the white (σw) and red (σr) noise are obtained by solving the following system of equations: where σ4 and σ1000 are the standard deviation taken over a sliding average of four and of one thousand points respectively. We estimated the white noise amplitude, σw = 0.003589, and the low-frequency red noise amplitude, σr = 0.0001397.

4.2. Transit timing

thumbnail Fig. 3

Distributions obtained from the prayer bead method for each fitted parameter. Vertical continuous lines correspond to the median of the distribution. Vertical dashed lines correspond to plus or minus 68% of the distribution.

Open with DEXTER

The measured central time of the transit can be compared to the expected transit time from a known ephemeris. We measure timing residuals for the transit according to the ephemeris of Knutson et al. (2009) (P = 2.21857578 ± 0.00000080 days, Tc = 2 454 399.23990 ± 0.00017 HJD). We find that the center of our transit is Tc = 2 454 430.310594 (HJD), which corresponds to an observed minus calculated (O−C) transit time of 47 s ± 4(±16). The uncertainties are set to the uncertainty in the observed transit time, while the values in parenthesis give the uncertainty in the predicted time. The total uncertainty in the O−C values is the sum of these two values. This corresponds to a time shifted of nearly 3σ after its expected value. However, it is difficult to draw conclusion from this result because the determination of mid-transit times are extremely sensitive to fitting the ingress and egress of the eclipse profile. Consequently, correlated noise or a nonuniform stellar brightness distribution could change the ingress/egress shapes of the transit and thus the central time away from what is expected for a uniform stellar disk. A recent study uses seven transits and seven eclipses of this exoplanet to derive a precise set of transit time measurements with an average accuracy of 3 s, and reveal a lack of transit-timing variations (Algol et al. 2010).

4.3. Limb darkening

Southworth (2008) shows that the determination of the light curve parameters, especially the radius of the planet, can be sensitive to the applied limb darkening model and its coefficients, yielding a possibility of systematic errors. This is particularly important for transits observed in the visible wavelength, but less important for transits observed in the infrared. Nevertheless, we investigated the importance of stellar limb darkening on the light-curve solutions and parameter uncertainties. We compared our results obtained from fits with various limb darkening laws and different numbers of fit coefficients.

We first note that the fit is clearly improved when including limb darkening (Fig. 6). We found that accounting for the effects of limb-darkening at 3.6 μm decreased the resulting best-fit transit depth by (ΔRp/R)3.6   μm = 0.0005, which corresponds to 5σ in the present result. We used a non linear limb darkening law with three fixed limb darkening coefficients (Table 1), which is slightly different from the one used in Désert et al. (2009). In the present work, we set c1 to zero and calculate new sets of coefficients c2, c3 and c4, resulting in a three-parameter non linear limb darkening law, (3)The c1(1 − μ1/2) term was removed because it reflects the intensity distribution at small μ values and is not needed when the intensity at the limb is desired to vary linearly at small μ values (for further details see Sing et al. 2009, 20103). Compared to the quadratic law, the added μ3/2 term provides the flexibility needed to more accurately reproduce the stellar model atmospheric intensity distribution at near-infrared and infrared wavelengths. The limb-darkening coefficients for the three-parameter non linear law were computed using a Kurucz ATLAS stellar model4 with Teff = 5000 K, log g = 4.5, and [Fe/H]  = 0.0 in conjunction with the transmission through the IRAC filters and fitted the calculated intensities between μ = 0.05 and μ = 1 (Table 1). This limb-darkening law is slightly different from the one used in the previous analysis (Désert et al. 2009) particularly at the limb (μ < 0.05). Nevertheless, the signal-to-noise ratio of our data is not high enough to distinguish the changes between the two sets of limb-darkening coefficients because the difference in the normalized transit light curves is less than 10-6.

The photometric precision of our transit light curve allowed us to perform fits using a non-linear limb-darkening model with the linear coefficient (c2) as a free parameter. We find c2 = 1.12994 ± 0.0243 in agreement with the value quoted in Table 1. As a test, we also fitted the transit light curve using a linear limb darkening law with one single coefficient c1 as free parameter. We derived c1 = 0.27821 ± 0.0291 and find that Rp/R is 2.5σ above the solution obtained with a the non linear limb darkening law, but with a less good χ2, showing that this solution can be excluded and the non linear law is preferred.

thumbnail Fig. 4

Delta-Khi-square distribution of parameters. Vertical dashed red lines correspond to the values for the minimal Khi-square and vertical continuous red lines correspond to plus or minus 1σ extracted with the prayer bead method. Vertical dotted green lines correspond to plus or minus 1σ extracted with the Khi-square method.

Open with DEXTER

Table 2

Fitted parameters of the transit light curves at 3.6 μm for visit 1 and 3.

thumbnail Fig. 5

Root-mean-square of binned residuals versus bin size. The solid continuous red line is proportional to N −1/2 and is normalized to match the value for bin size N = 4. The top continuous green curve corresponds to bin residuals without pixel-phase decorrelation. Black crosses correspond to bin residuals with pixel-phase effect decorrelation. Remaining red-noise becomes apparent for bins larger than 50 data points (6.25 s).

Open with DEXTER

thumbnail Fig. 6

Residuals from fits with or without limb-darkening corrections (bottom and top panel respectively). Data (diamonds) are binned by 120 (15 s) and a smooth of these binned residuals is overplot in red. Vertical dashed lines mark the beginning and the end of the transit (first and last contact). The fit is better with limb darkening. No signature remains in the residuals.

Open with DEXTER

thumbnail Fig. 7

Top panel: normalized transit light curves obtained in stellar mode (diamonds) and in subarray mode (crosses). The subarray mode transit light curve is shifted vertically for display purpose. The subarray mode transit light curve is binned by 64 to obtain the same number of data points (~1900 points) as the stellar mode transit light curve. Bottom panel: residuals obtained from the fits of the transit light curve. Residuals from the binned subarray mode transit light curve are shifted vertically (for display) and are four times smaller.

Open with DEXTER

thumbnail Fig. 8

Joint distributions of parameters obtained from the prayer bead method using the subarray mode observations. The x-axis corresponds to radius ratios for all the panels. The top panel is the fit of the central time, the middle and bottom panel, correspond to the impact parameter and the scale of the system respectively. The diamonds with error bars correspond to the results obtained from the observations of visit 1 (Désert et al. 2009). The dashed red lines correspond to linear fits of the points from visit 3, and reveal the correlations between Rp/Rb and a/R.

Open with DEXTER

thumbnail Fig. 9

Orbital parameters derived from several studies. Top panel: impact parameter b; bottom panel: system scale a/R. Diamonds are values from several other studies: data are from Pont et al. (2007) in the visible, from Sing et al. (2009) in the near infrared and from Knutson et al. (2007a,b) at 8.0 μm also measured by Carter & Winn (2010). Blue squares are from our visit 1 and 2 observations in the four IRAC channels (Désert et al. 2009). The red square are from the observations of visit 3 (this work). Continuous and dashed straight black lines indicate the average and error values obtained from data excluding outliers. These values are considered as best estimates for b and a/R. Parameters derived from visit 1 at 3.6 and 5.8 μm respectively are significantly different from the best estimates. The empty squares correspond to the derived parameters using to the observation of visit 1 and the transit light curve models which take into account the occulted spots during ingress or egress.

Open with DEXTER

thumbnail Fig. 10

(Top) a): relative flux variation for HD 189733 observed with ground-based follow-up encompassing the time of our IRAC observations (connected crosses). The black circular spot indicates the time of our Spitzer stellar mode observations at 3.6 and 5.8 μm (visit 1). The variations are caused by rotational modulation in the visibility of star spots with a rotation period of 11.953 days (Henry & Winn 2008). The continuous curve overplotted in red corresponds to a sinusoid with a period of 11.95 days. (Bottom) b): the same as above for the stellar mode observations at 4.5 and 8.0 μm (visit 2) and subarray mode observations at 3.6 μm (visit 3).

Open with DEXTER

4.4. Comparison with previous observations

We compare the parameters derived here with the results obtained from the observations of visit 1 and 2 (Désert et al. 2009). Results from fits to the transit light curve obtained in the observations of visit 1 gathered at 3.6 μm are summarized in Table 2. We measured a dispersion of 2.1 × 10-3 (rms) on the individual  ~1900 data points with an exposure time of 0.4 s. This corresponds to a S/N of 400:1 per image. To compare with our present data set, we binned the transit light curve of visit 3 by 64 and obtained  ~1900 points with a rms = 4.9 × 10-4. The two transit light curves, obtained for visit 1 and 3 at 3.6 μm, are overplotted with 1900 points each in Fig. 7. The S/N by subarray frame is four times higher than by stellar mode exposure, in agreement with expected photon noise for a total exposure time sixteen times longer.

We are particularly interested in the comparison of the ratios of the radius of the planet over the radius of the star at 3.6 μm in the two observation sets. The radius ratio found in visit 3 stands 4σ above the value derived from visit 1 at 3.6 μm (Table 2). The impact parameters b and system scales a/R measured in the transit light curves provide information for interpreting these results. We found that the b and a/R values derived from visit 1 observations at 3.6 μm, are respectively 3σ below and above the values obtained here at the same wavelength with visit 3 observations (Figs. 8 and 9). The same results are obtained for visit 1 observations at 5.8 μm (both observations at 3.6 μm and 5.8 μm were accumulated simultaneously during visit 1). As already noticed in Désert et al. (2009), the b and a/R measured with visit 1 and 2 stellar mode observations are consistent between channels observed simultaneously, but disagree between channels observed at two different epochs. This suggests that unrecognized systematics remains between the observations obtained at the two epochs. This systematic effect could be either instrumental or astrophysical. Below, we estimate the impact of those systematics on the radius determination, and possible astrophysical origins of these discrepancies.

The value of the impact parameter b has a direct effect on the derived planetary radius. This can be seen theoretically from the analytic approximation (Carter et al. 2008) (4)where τ is the ingress or egress duration(τ = tII − tI), tI and tII are the time of first and second contact, and T is the total transit duration. Consequently, underestimating b from 0.66 to 0.63, without changing the timing of ingress and egress would lead to overestimate Rp/R by 0.001, which corresponds to 10σ in the result from visit 3.

We estimated the influence of systematic errors on the estimate of b on the determination of the radius.

We fixed the b value of visit 1 to the one obtained at higher S/N in visit 3 (b = 0.661) and fitted visit 1 transit light curves at 3.6 μm and 5.8 μm. In that case we obtained (Rp/R)3.6  μm = 0.15532 for visit 1; these values agree with the results obtained in visit 3 (Fig. 9).

In conclusion, we find consistent values for Rp/R and a/R in all data sets if the impact parameter b is set to the same value. Our final value agree with the values b = 0.671 ± 0.008 and a/R = 8.92 ± 0.09 derived by Pont et al. (2008) with HST/ACS, with the values obtained by Sing et al. (2009) with HST/NICMOS, and with a/R = 8.924 ± 0.022 from Carter & Winn (2010) using several Spitzer/IRAC transit light curves at 8 μm. This suggests that b and a/R derived from visit 1 are affected by unknown systematics, which could be either instrumental or astrophysical such as starspots (see Sect. 4.5).

4.5. Impact of starspots on the derived parameters

The source HD 189733 is an active K star, which has been observed to vary photometrically by  ±1.5% at visible wavelengths (Henry & Winn 2008; Croll et al. 2007; Miller-Ricci et al. 2008). Those variations are caused by rotational modulations in the visibility of star spots on a rotation period of 11.953 ± 0.009 days (Hébrard & Lecavelier des Etangs 2006; Henry & Winn 2008).

In this section, we consider the possibility that transit light curves can be modified by the presence of spots on the surface of the star, introducing systematic errors on the parameter estimates. In particular, we investigate the possible contribution of stellar spots and brightness inhomogeneities to explain the discrepancy between the planet radii at 3.6 μm measured in visits 1 and 3 data (Sect. 4.5.2). We also investigate the effect of spots to explain the bias in the orbital parameters estimated using visit 1 data (Sect. 4.5.3).

4.5.1. Ground based photometric follow-up

We use ground based photometric follow-up to estimate and to correct for the stellar activity at the epochs of our three visit observations.

The star has been followed-up in V-band over several years with the Tennessee State University 0.8 m automated photometric telescopes (APT) at Fairborn Observatory (Henry et al. 2008). These data show that the stellar flux variation reach a peak-to-peak maximum of 3% in the visible (Fig. 10).

The observed variations in the visible have to be translated into a corresponding variation in the infrared. Let us define fλ, the relative variation in flux raising to stellar activity: (5)where F is the measured flux and Fquiet is the reference stellar flux (e.g., without stellar activity). fλ is negative when the star is fainter and positive when it is brighter than the reference brightness.

Previous HST ACS observations of HD 189733 have established that its stellar spots have an effective temperature approximately 1000 K cooler than that of the stellar photosphere (Pont et al. 2008) with a spot size of  ~2.8% of the stellar disk area. Following Knutson et al. (2009), we considered that the amplitude of the stellar flux variations scales approximately as the ratio of the flux of blackbodies at 5000 and 4000 K. Thus, we used the normalized difference of Planck functions to estimate the maximum variation of flux at 3.6μm f(3.6   μm). The 3% maximum variation in the V-band caused by the presence of spots on the stellar surface scales to f(3.6   μm) ~ 0.8% in the infrared.

The ground-based photometry suggests that visit 1 happened during the maximum brightness of the star (Fig. 10). Unfortunately, no ground-based observations were acquired during visit 2 and 3. Thus, we have no direct measurement of the stellar activity at those precise epochs. Nevertheless, the activity period P = 11.953 ± 0.009 (Henry & Winn 2008) can be used to interpolate the star brightness at the epochs of visit 2 and 3. We found that the star should be close to its maximal brightness during visit 2, and  ~2% below the maximal brightness during visit 3, corresponding to f(3.6   μm) ~ −0.5% at this epoch.

4.5.2. Stellar spots outside the zone occulted by the planet

thumbnail Fig. 11

Census of the measured values of Rp/R ratios in the four Spitzer/IRAC channels for HD 189733b obtained using three transit observations during visit 1 and 2 (blue circles, Désert et al. 2009), and visit 3 (red square; this work) with their 1σ error bars. The data point obtained for visit 3 has to be corrected from stellar spots (see text 4.5.3). The IRAC bandpasses for each channel are shown at the top with dotted line. The theoretical planet-to-star radius ratios with absorption by only water molecules or by only Rayleigh scattering are plotted in light blue and grey respectively.

Open with DEXTER

thumbnail Fig. 12

Measured values of Rp/R ratios at 3.6 μm as fonction of the relative decrease of stellar flux due to spots (in V-band). The blue and red squares are the Rp/R derived from the observations of visit 1 and 3 respectively. The dotted line corresponds to the theoretical correction of the relative difference between the measured radius ratio and the expected value as function of the relative decrease of stellar flux in the V-band fλ. Using ground based photometric follow-up (Fig. 10), we set f(3.6   μm) to zero for visit 1 (see text 4.5.1). The Rp/R we derive from the observations of visit 3 (no correction) is above the one obtained in visit 1 (4σ). To compare both visits, the Rp/R derived from visit 3 has to be corrected for spots (Sect. 4.5.3). Simultaneous observation in V-band shows that fV − band =  −2% during visit 3 (Fig. 10). (Rp/R)Corrected for visit 3 agree with Rp/R derived in visit 1 for a slope α =  −3.

Open with DEXTER

Spots and bright faculae cause variations in the star brightness (F). Nonetheless, if the surface brightness of the zone occulted by the planet is not modified, the absolute depth of the transit light curve (ΔF) remains the same. Therefore, at epochs of low stellar brightness like during visit 3, the relative transit depth (ΔF/F) can be significantly overestimated, leading to an overestimate of the planet-to-star radius ratio.

As seen in Corot-2b (Czesla et al. 2009), the transit depth constraining the planet radius ratio is correlated with the star brightness. We can empirically define α by: (6)where (Rp/R)True is the true radius ratio that could be obtained from the transit light curve when the star presents no spots and is at the reference brightness. The relative error on the measured radius ratio is therefore given by (7)In the Corot-2b data (Czesla et al. 2009), we found α ~ 1.7. With the hypothesis that the stellar surface brightness outside the spot areas is not modified by the stellar activity, we obtained α =  −1. With the photometric variations obtained from ground-based observation at epochs close to the Spitzer observations, we found f3.6   μm =  −0.5% in visit 3 compared to visit 1. Therefore, the radius ratio Rp/R derived from visit 3, which is  ~0.75 ± 0.18% above the radius derived from visit 1 (Table 2), can be explained using Eq. (7) if α ~ −3 ± 1 (Fig. 12).

In short, the discrepancy between the visit 1 and 3 radius ratios at 3.6 μm can be explained by the effect of stellar activity and the presence of spots during the visit 3.

4.5.3. Occulted starspots

The planet occulting dark or bright areas on the surface of the star during transit produces a rise or decrease in the flux. Because the ingress and egress parts of the transit light curve strongly constrain the orbital parameters of the planet, dark or bright areas occulted by the planet during the ingress or egress change the apparent time of the transit contacts, introducing systematic errors in the derived b and a/R parameters (Eq. (4)). Consequently, the low b value obtained in channel 1 and 3 using visit 1 data could be caused by the presence of cold spots or bright faculae occulted by the planet.

To estimate the influence of occulted star cold or bright regions during the ingress and egress on the determination of our parameters, we fitted the visit 1 and 2 observations with a transit light curve model including occulted star spots (cold or hot) during ingress and egress. These fits have a goodness-of-fit equivalent to the fits using the model without spots. For visit 2, the model with occulted spots gives similar results as Désert et al. (2009). Most importantly, for visit 1, the model well fitsl the data with a bright spot occulted during the ingress part of the light curve. The best-fit is found assuming a spot size similar to the size of the planet and a spot surface brightness 8% brighter than the star. Compared to the results given in Désert et al. (2009), this new model applied to visit 1 data at 3.6 μm gives a larger b by 0.013 (2σ), a smaller a/R by 0.06 (1σ), and error bars larger by a factor of about 1.5. Interestingly, the Rp/R found with this new model change by less than 0.3σ. Similar conclusions are obtained with data at 5.8 μm (Fig. 9).

As a conclusion, the surprising results for b and a/R found by Désert et al. (2009) in visit 1 data can simply be explained by astrophysical systematics in the transit light curve due to stellar activity. However, the difference in the measured values for the radius ratios Rp/R between visit 1 and 3 cannot be explained by these dark and bright areas occulted during ingress or egress. It must be caused by the fluctuation of the whole stellar disk brightness (see Sect. 4.5.2).

4.6. Opacities in the IRAC bandpasses

The Rp/R measured at 3.6 μm from visit 1 is 4σ above the value measured in visit 3. This difference can be reconciled by correcting for the 2% stellar variability monitored with the ground-based observations between the two visits. However, at least two colors measured simultaneously are required to ensure that the measured variations of the wavelength-dependent radius of the planet are caused by atmospheric chemical absorbers. In this framework, the Rp/R measured at 3.6 μm from visit 1 is 4σ above the opacity due to water if we consider that the radius simultaneously measured at 5.8 μm is due to water opacity. Therefore, other species beyond H2O must be present in the atmosphere of the planet and absorb at 3.6 μm. One possible physical explanation for the large radius at 3.6 μm could be that it is due to Rayleigh scattering (Lecavelier des Etangs et al. 2008b). Effectively, the radius mesured at this wavelength agrees at one sigma level with its expected value if we consider Rayleigh scattering opacity (Fig. 11). We note that the possible presence of a strong methane (CH4) band at 3.3 μm could also contribute to the opacity in this bandpass (Sharp & Burrows 2007). The radius we derived at 4.5 μm could be interpreted as due to COmolecules that exhibit large opacities in this bandpass (Désert et al. 2009). The radii we derived at 5.8 μm and 8.0 μm are at more than 3σ above the absorption due to Rayleigh scattering, and could be interpreted by the presence of water molecules whose opacity overcomes Rayleigh scattering at those wavelengths (Fig. 11). These two bandpasses exhibit approximatly the same radius value, consistent with a model which includes water.

Therefore water vapour could be present in the atmosphere of this hot-Jupiter, but it cannot be proved simply by comparing the planetary radii measured at 3.6 and 5.8 μm, because the radius measured at 3.6 μm is affected by some other species absorption.

5. Conclusion

We obtained a new high-S/N photometric transit light curve at 3.6 μm using the IRAC/Spitzer subarray mode observations. This observing mode offers a higher photometric cadence, therefore a better signal-to-noise ratio than the previous measurements (stellar mode) obtained at this wavelength. We derive a value of Rp/R = 0.15566 ± 0.00011 at 3.6 μm, which is 4σ above the value derived from previous observations in the same bandpass (Désert et al. 2009). This difference is explained by the stellar variability between the two observational epochs. Simultaneous ground-based observation allowed us to correct for these variations. A photometric variation of  ~2% in V-band was monitored between the two measurements and provides a correction of  −0.75% for the new 3.6 μm measurement. Note that the first observation was obtained during the maximum brightness of the star, which led to the minimal possible value of the measured Rp/R = 0.1545 ± 0.0003 at 3.6 μm for this planet (Désert et al. 2009). This final planet-to-star radius ratio at 3.6 μm is above the expected value if only water molecules were contributing to the absorption at both 3.6 and 5.8 μm.

We showed the importance of a long term ground-based monitoring program in order to take into account the stellar activity and to compare results obtained from observations carried out at different epochs. Any estimate of the planet-to-star radius ratio should be associated with a corresponding stellar activity level.


Acknowledgments

D.K.S. is supported by CNES. 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. D.E. is supported by the Centre National d’Études Spatiales (CNES). D.E. acknowledges financial support from the French Agence Nationale pour la Recherche through ANR project NT05-4_44463. We thank our anonymous referee for his comments that strengthen the presentation of our results.

References

All Tables

Table 1

Limb-darkening coefficients used for IRAC/channel 1 bandpass centered at 3.6 μm.

Table 2

Fitted parameters of the transit light curves at 3.6 μm for visit 1 and 3.

All Figures

thumbnail Fig. 1

X (black) and Y (red) position of the star centroid in pixels as function of time (binned by 4). These positions are used in the decorrelation function that is used for the fit of the transit light curve.

Open with DEXTER
In the text
thumbnail Fig. 2

Top panel: raw transit light curve obtained in subarray mode binned by four with its fit overplotted in red. Pixel-phase effect is clearly visible. Middle panel: normalized transit light curve. Bottom panel: residuals obtained from the fits of the transit light curve.

Open with DEXTER
In the text
thumbnail Fig. 3

Distributions obtained from the prayer bead method for each fitted parameter. Vertical continuous lines correspond to the median of the distribution. Vertical dashed lines correspond to plus or minus 68% of the distribution.

Open with DEXTER
In the text
thumbnail Fig. 4

Delta-Khi-square distribution of parameters. Vertical dashed red lines correspond to the values for the minimal Khi-square and vertical continuous red lines correspond to plus or minus 1σ extracted with the prayer bead method. Vertical dotted green lines correspond to plus or minus 1σ extracted with the Khi-square method.

Open with DEXTER
In the text
thumbnail Fig. 5

Root-mean-square of binned residuals versus bin size. The solid continuous red line is proportional to N −1/2 and is normalized to match the value for bin size N = 4. The top continuous green curve corresponds to bin residuals without pixel-phase decorrelation. Black crosses correspond to bin residuals with pixel-phase effect decorrelation. Remaining red-noise becomes apparent for bins larger than 50 data points (6.25 s).

Open with DEXTER
In the text
thumbnail Fig. 6

Residuals from fits with or without limb-darkening corrections (bottom and top panel respectively). Data (diamonds) are binned by 120 (15 s) and a smooth of these binned residuals is overplot in red. Vertical dashed lines mark the beginning and the end of the transit (first and last contact). The fit is better with limb darkening. No signature remains in the residuals.

Open with DEXTER
In the text
thumbnail Fig. 7

Top panel: normalized transit light curves obtained in stellar mode (diamonds) and in subarray mode (crosses). The subarray mode transit light curve is shifted vertically for display purpose. The subarray mode transit light curve is binned by 64 to obtain the same number of data points (~1900 points) as the stellar mode transit light curve. Bottom panel: residuals obtained from the fits of the transit light curve. Residuals from the binned subarray mode transit light curve are shifted vertically (for display) and are four times smaller.

Open with DEXTER
In the text
thumbnail Fig. 8

Joint distributions of parameters obtained from the prayer bead method using the subarray mode observations. The x-axis corresponds to radius ratios for all the panels. The top panel is the fit of the central time, the middle and bottom panel, correspond to the impact parameter and the scale of the system respectively. The diamonds with error bars correspond to the results obtained from the observations of visit 1 (Désert et al. 2009). The dashed red lines correspond to linear fits of the points from visit 3, and reveal the correlations between Rp/Rb and a/R.

Open with DEXTER
In the text
thumbnail Fig. 9

Orbital parameters derived from several studies. Top panel: impact parameter b; bottom panel: system scale a/R. Diamonds are values from several other studies: data are from Pont et al. (2007) in the visible, from Sing et al. (2009) in the near infrared and from Knutson et al. (2007a,b) at 8.0 μm also measured by Carter & Winn (2010). Blue squares are from our visit 1 and 2 observations in the four IRAC channels (Désert et al. 2009). The red square are from the observations of visit 3 (this work). Continuous and dashed straight black lines indicate the average and error values obtained from data excluding outliers. These values are considered as best estimates for b and a/R. Parameters derived from visit 1 at 3.6 and 5.8 μm respectively are significantly different from the best estimates. The empty squares correspond to the derived parameters using to the observation of visit 1 and the transit light curve models which take into account the occulted spots during ingress or egress.

Open with DEXTER
In the text
thumbnail Fig. 10

(Top) a): relative flux variation for HD 189733 observed with ground-based follow-up encompassing the time of our IRAC observations (connected crosses). The black circular spot indicates the time of our Spitzer stellar mode observations at 3.6 and 5.8 μm (visit 1). The variations are caused by rotational modulation in the visibility of star spots with a rotation period of 11.953 days (Henry & Winn 2008). The continuous curve overplotted in red corresponds to a sinusoid with a period of 11.95 days. (Bottom) b): the same as above for the stellar mode observations at 4.5 and 8.0 μm (visit 2) and subarray mode observations at 3.6 μm (visit 3).

Open with DEXTER
In the text
thumbnail Fig. 11

Census of the measured values of Rp/R ratios in the four Spitzer/IRAC channels for HD 189733b obtained using three transit observations during visit 1 and 2 (blue circles, Désert et al. 2009), and visit 3 (red square; this work) with their 1σ error bars. The data point obtained for visit 3 has to be corrected from stellar spots (see text 4.5.3). The IRAC bandpasses for each channel are shown at the top with dotted line. The theoretical planet-to-star radius ratios with absorption by only water molecules or by only Rayleigh scattering are plotted in light blue and grey respectively.

Open with DEXTER
In the text
thumbnail Fig. 12

Measured values of Rp/R ratios at 3.6 μm as fonction of the relative decrease of stellar flux due to spots (in V-band). The blue and red squares are the Rp/R derived from the observations of visit 1 and 3 respectively. The dotted line corresponds to the theoretical correction of the relative difference between the measured radius ratio and the expected value as function of the relative decrease of stellar flux in the V-band fλ. Using ground based photometric follow-up (Fig. 10), we set f(3.6   μm) to zero for visit 1 (see text 4.5.1). The Rp/R we derive from the observations of visit 3 (no correction) is above the one obtained in visit 1 (4σ). To compare both visits, the Rp/R derived from visit 3 has to be corrected for spots (Sect. 4.5.3). Simultaneous observation in V-band shows that fV − band =  −2% during visit 3 (Fig. 10). (Rp/R)Corrected for visit 3 agree with Rp/R derived in visit 1 for a slope α =  −3.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.