A Search for a Sub-Earth-Sized Companion to GJ 436 and a Novel Method to Calibrate Warm Spitzer IRAC Observations

We acquired observations of GJ 436 nearly continuously during 2008 May 5–29, interrupted for several hours at approximately 2 day intervals for data downloads. A complete description of the EPOXI photometric extraction pipeline is given in Ballard et al. (2010) and summarized here. We used the existing Deep Impact data-reduction pipeline to perform bias and dark subtractions, as well as preliminary flat-fielding (Klaasen et al. 2005). We first determined the position of the star on the CCD using PSF-fitting, by maximizing the goodness of fit (with the χ² statistic as an estimator) between an image and a model PSF (oversampled by a factor of 100) with variable position, additive sky background, and multiplicative brightness scale factor. We then processed the images to remove sources of systematic error due to the CCD readout electronics. We first scaled down the two central rows by a constant value, then we scaled down the central columns by a separate constant value, and finally we scaled the entire image by a multiplicative factor determined by the size of the subarray. We performed aperture photometry on the corrected images, using an aperture radius of 10 pixels, corresponding to twice the half-width at half-maximum of the PSF. To remove the remaining correlated noise due to the interpixel sensitivity variations on the CCD, we fit a 2D spline surface, with the same resolution as the CCD, to the brightness variations on the array as follows. We randomly drew a subset of several thousand out-of-transit and out-of-eclipse points from the light curve (from a data set of ∼29,988 points) and found a robust mean of the brightness of the 30 nearest neighbors for each. We fit a spline surface to these samples and corrected each data point individually by linearly interpolating on this best-fit surface. We used only a small fraction of the observations to create the spline surface in order to minimize the potential transit signal suppression introduced by flat-fielding the data with itself. To produce the final time series, we iterated the preceding steps, fitting for the row and column multiplicative factors, the subarray size scaling factor, and the 2D spline surface that minimized the out-of-transit white noise of the photometric time series. We included one additional step to create the final 2D spline, which was to iteratively remove an overall modulation from the GJ 436 light curve that we attributed to star spots. After we took these steps to address the systematics associated with the observations, the red noise was largely removed. Figure 1 shows the GJ 436 time series before and after the 2D spline correction. After the correction is applied, the precision of the light curve is 56% above the photon limit. We note that the version of the EPOXI GJ 436 light curve presented in Ballard et al. (2010) is very slightly different from the version used in this analysis. Because of the possibility of suppression of additional transits by the 2D spline correction method, we produced a version of the light curve that masked the points that occurred during transit times of the putative GJ 436c from contributing to the CCD sensitivity map. This set of additional masked points is 16 hr total in duration over the entire data set, consisting of 2 hr intervals centered at each of 10 candidate transit events—two of which partially coincide with transits of GJ 436b and were already masked from the 2D spline correction.

Zoom In Zoom Out Reset image size

Under Spitzer GO Program 50056 (PI: H. Knutson), GJ 436 was monitored for 70 hr continuously from 2008 July 12–15 by the Spitzer Space Telescope (Werner et al. 2004) with the IRAC subarray (Fazio et al. 2004) in the 8 μm channel, to obtain a full phase curve of the hot Neptune GJ 436b. This observing sequence consisted of 0.4 s integration exposures in 9195 blocks of 64 images each. These data were preflashed using the same technique as the HD 149026 observations (described in Knutson et al. 2009b). Preflashing effectively removes most of the "detector ramp" effect, which is characterized by an initial upward asymptote in the measured flux, followed by a gradual downward slope (Charbonneau et al. 2008; Knutson et al. 2008, Knutson et al. 2009a. In this case we observed a region in M17 centered at R.A. 18^h20^m28^s and decl -16°12'20'' with fluxes ranging between 3500–8000 MJy Sr^-1 for 30 minutes prior to the start of the observations. This significantly reduced the amplitude of the detector ramp in these data, although there is still a small rise in flux of approximately 0.05% visible in the first 5 hr of observations. The data after this point have a flux variation of less than 0.05% over the remaining 65 hr. The photometric extraction was performed by a method similar to the one described in detail in Knutson et al. (2009b), which we summarize here. We determined the position of the star on the array by centroid, taking a position-weighted sum of the flux within a circular aperture with a radius of 5 pixels. We find that using this aperture size for the position determination returned the photometric lowest rms in the resulting time series. To perform background subtraction, we first created a median image from all the observations and identified four regions in the corners of the array where the light from the point-spread function (PSF) was minimized, in order to minimize contamination by diffraction spikes in our background estimates. We created a histogram of the values of these pixels and calculated the background from the mean of this distribution. We then performed aperture photometry with apertures ranging from 2.5 pixels in radius to 7.0 pixels in radius (the radius used for the aperture photometry affects the significance of the putative additional transit, which is explained in § 3.2). We discarded points that lie more than 3σ from the median value within each set of 64 images (in order to remove images affected by transient hot pixels). For the purposes of addressing the additional planet hypothesis, we were concerned with only a 5 hr portion of the light curve that fell in the latter half of the 70 hr phase curve, where we no longer observe any evidence for a detector ramp. Nonetheless, we fit and divided a quadratic function in time during the window of interest in order to remove any remaining trends in the data, which could be due to spots on the star, position-dependent aperture losses, and the phase variation from the planet itself.

We observed GJ 436 in subarray mode, using a 0.1 s integration time in the 4.5 μm channel. These observations span 18 hr over UT 2010 January 28–29. The observing sequence consisted of 7640 blocks of 64 images each. We experimented with two methods of locating the position of the star on the array. We first found a flux-weighted sum of the position within a circular aperture with a radius of 3 pixels. We also experimented with determining the position of the star by fitting a PSF within a 3.5 pixel aperture (we found similar results with larger apertures). Using the 100-times-oversampled PSF provided by the Spitzer Science Center for the 4.5 μm channel, we performed a χ² minimization in which we allowed the X and Y positions of the PSF, a multiplicative brightness factor, and an additive sky background value to vary. We compared a histogram of the positions to get a sense for the precision of each measurement technique, first subtracting a running median calculated individually for each point from the nearest 20 points in time (to account for positional drift). We determined that these histograms had similar widths with the two methods: 6.0 × 10^-3 pixels in Y with PSF-fitting versus 6.4 × 10^-3 pixels in Y with centroid. However, the precision of the final time series was not improved by using the PSF-fitting measured positions; rather, the precision was degraded by 20%. We attribute this degradation to the decreased positional precision of the PSF positions: although their bulk scatter is less than the scatter from centroiding, the precision with which we measure an individual position is set by the resolution of the PSF at 0.01 pixels, rather than the 1σ error from the scatter in the centroided positions of 6.4 × 10^-3 pixels. We considered the possibility of improving our positional accuracy by producing a more highly oversampled PSF with an interpolation of the Spitzer Science Center PSF, but found the computing time of PSF-fitting with such a large PSF to be prohibitively expensive.

We then measured the stellar brightness in each image by performing aperture photometry on the basic calibrated data products across a range of apertures from 2.1 to 6 pixels. We calculate the sky background, which is almost negligible, from a 3σ-clipped mean of the pixels inside a ring with width of 10 pixels from a radius of 7 pixels to 17 pixels. We found that the precision of the final time series was optimized at an aperture of 2.1 pixels in radius.

We then corrected for the well-known intrapixel sensitivity variation observed in IRAC channels 1 and 2 (Reach et al. 2005; Charbonneau et al. 2005; Knutson et al. 2008, Knutson et al. 2009b). In lieu of fitting a polynomial in X and Y to the brightness variations, we instead implemented a point-by-point correction. We first binned the light curve into 20 s bins (approximately 145 points bin^-1). For each binned point, we evaluate a weighted-sensitivity function using all unbinned points (excluding those points that occurred inside the bin in question and outliers more than 3σ from the mean of nearby points) given by

where f_j is the flux value of the j th observation, which is assigned a weight based on its distance in X and Y from the i th point being corrected. The function S(t_j) is a boxcar function in time, which is 1 for all points that are permitted to contribute to the sensitivity map, and 0 for all points that are not. The position of this function defines a "mask" interval during which time observations do not contribute to the intrapixel sensitivity map. The purpose of this mask is to exclude points that do not reflect an accurate representation of the pixel sensitivity, such as those during transit; if we were to include these points in the creation of the intrapixel sensitivity map, we would both suppress the transit and introduce additional correlated noise outside of transit. The sensitivity function is then normalized by dividing by the Gaussian function multiplied by the boxcar S(t_j). We find the best results using σ_x = 0.017 pixels and σ_y = 0.0043 pixels (the width of the Gaussian weighting function is much smaller in Y, because the dependence of the flux on the Y position is much stronger). We then divide the i th binned flux value by W(x_i,y_i) to remove the effects of intrapixel sensitivity. Using a point-by-point sensitivity function, we do not need to assume a functional form for the intrapixel sensitivity; however, there is the very important caveat that flat-fielding the data by itself has the effect of suppressing the depths of additional transits, if they are present. We discuss how we avoid this effect in § 3.3 with the use of the mask function. In Figure 2 we show three-dimensional views of the weighted-sensitivity function W(x_i,y_i). The large-scale features of W(x_i,y_i) are well approximated by polynomials in X and Y, per previous techniques, but we also find a smaller-scale "corrugation" effect in the Y direction, where the sensitivity exhibits low-level sinusoidal-like variations with a separation of approximately 5/100 of a pixel between peaks.

Zoom In Zoom Out Reset image size

We investigated the authenticity of these features using several tests. For Gaussian noise, we expect the weighted-sensitivity function to exhibit only smoothed random noise, with features equal in size to the smoothing kernel, and so the corrugation features in Y should have a size scale near 9/1000 of a pixel (because σ_y is 0.0043 pixels). Futhermore, we should be able to predict the χ² improvement by comparing the data to W(x_i,y_i) as opposed to the null hypothesis. A Gaussian time series will be better fit by a smoothed version of itself than by a flat line. The χ² should improve by the number of smoothing kernels contained in the interval in question: for a time series with N_y values ranging over 0.25 of a pixel and a smoothing kernel of σ_y = 0.0043, then the number of smoothing kernels (defining 3σ from the center as the extent of the kernel) contained in the interval is approximately 10, and we should expect the χ² of the fake time series to improve by 10 when compared to the weighted-sensitivity function instead of a flat line. Conversely, in order for authentic features attributable to the intrapixel sensitivity to be believable, the features must be significantly larger than the smoothing kernel, and the weighted-sensitivity function must provide a much better fit to the data than predicted from a smoothed version of the data itself. To test this hypothesis, we created a random Gaussian data set sampled at the same X and Y positions on the detector, with a standard deviation equal to that of the actual time series. We then created a sensitivity function using the same kernel in X and Y as for the actual observations and compared the results; this comparison is shown in Figure 3.

Zoom In Zoom Out Reset image size

We find that the weighted-sensitivity function from the random data set looks as we expected it to look, with corrugation features near the size of the full width at half-maximum (FWHM) of the smoothing kernel. The improvement in χ², using a flat line with value 0 as compared to the binned weighted-sensitivity function, is 11.4 for the fake data set, very close to the predicted improvement of 10. For comparison, the improvement in χ² for the real data between a flat line and the weighted-sensitivity function is 24.3, and the amplitudes of the features are larger. Furthermore, the peak-to-peak size scale of features is 0.05 pixels, more than six times the FWHM of the smoothing kernel, which argues against their being smoothing artifacts. The difference in width of the weighted-sensitivity function between the real and fake data sets (the thickness of the black points in the vertical direction) can be attributed to the weak X dependence of the pixel sensitivity. In the lowest right-hand panel of Figure 2, there is evidence for pixel sensitivity variation in the X direction over a size scale comparable to the width of the 0.04 pixel cross section used to generate Figure 3. For these reasons, we are confident that the high-frequency features, with a period of 0.05 pixels, in the real weighted-sensitivity function are authentic. The amplitude of this effect is 100 ppm: 40% of the amplitude of a 0.75 R_⊕ transit in front of a 0.438 R_⊙ star such as GJ 436 (Ballard et al. 2010). Therefore, removing this correlation with position was crucial to our ability to rule out the putative transit of depth 250 ppm.

With this photometric-reduction procedure, we achieve a precision of 0.0053 per 0.1 s exposure on the unbinned time series. Compared to the photon noise-limited precision of 0.0042, we are 26% above the photon noise limit, although the presence of remaining correlated noise means that the scatter with larger bin sizes deviates more from the ideal Gaussian limit. We achieve a sensitivity of 71 ppm per 20 minute bin, compared to the shot-noise limit of 41 ppm at that bin size.

For comparison, we also perform a reduction using a polynomial fit in X and Y to remove the intrapixel sensitivity variations in flux (Reach et al. 2005; Charbonneau et al. 2005; Knutson et al. 2008, Knutson et al. 2009b). We express the measured flux f^' in terms of the incident flux f and express the X and Y positions of the star on the detector with the following expression:

We find that the precision we achieve using a polynomial intrapixel sensitivity function for the same bin size of 20 minutes is 230 ppm, as compared to a precision of 71 ppm using the weighted-sensitivity map W(x_i,y_i). If we divide the time series into three portions of 5 hr duration and fit separate polynomial coefficients for each portion, we achieve a precision of 91 ppm for a bin size of 20 minutes—still 1.3 times larger than the precision using W(x_i,y_i). However, although we can improve the overall precision of the time series by fitting the polynomial coefficients independently for increasingly short durations, we are never able to reliably recover transits of a 0.75 R_⊕ planet in a time series reduced with a polynomial sensitivity function. We discuss this analysis in § 3.3.

In Ballard et al. (2010), we conducted a search for additional transiting planets in the EPOXI light curve of GJ 436. In that work, we demonstrate our sensitivity to additional transits by injecting light curves corresponding to additional planets with varying planetary radius, period, and phase, and then attempting to recover them by maximizing the χ² goodness of fit. We found, when we carefully accounted for the signal suppression introduced by reducing the data with the 2D spatial spline method, that we were sensitive to Earth-sized planets with good (≥50%) probability for periods less than only 0.5 days. However, we discovered weak evidence for an additional transiting planet, which fell well below the criterion we established for a detection. In Ballard et al. (2010), we empirically established the criterion for detection that used the improvement of χ² corresponding to the best-fit transit signal compared to the χ² of the null hypothesis; we could reliably recover the correct period of any signal which produced an improvement of Δχ²≥250. The transit signal corresponding to a 0.75 R_⊕ sized planet is well below this threshold. However, we find that the largest deviations occurred with a regular period near the 4∶5 resonance with the Neptune-sized planet. Five of these events produced a combined improvement to the χ² of 140. However, these five comprise only half of the transits that we would expect to see with this ephemeris: Of these remaining five events, two coincide with transits of GJ 436b, one occurs during a gap in the phase coverage, and two show no evidence of a transit. If we include all predicted candidate transit events in our χ² calculation (including the two events that coincide with a transit of GJ 436b), the null hypothesis gives a better solution than any transit model. However, if points that coincide closely in time with transits of the GJ 436b are excluded from the calculation, the improvement over the null χ² is 70. There are two motivations to exclude these in-transit points: First, we fit a slope with time to the points immediately outside of the transit of GJ 436b (from 3 minutes to 30 minutes before the start of transit and after the end of transit) and divide this slope in order to normalize each transit before we fit for the system parameters (see Ballard et al. 2010), and this procedure may suppress other signals. Second, there is also the small possibility of an occultation of GJ 436c by GJ 436b.

We observed comparable Δχ² improvement (within 2σ) for periods ranging between 2.1074 days (0.797 the period of GJ 436b) and 2.1145 days (0.800 the period of GJ 436b). Figure 4 shows the 10 events separately. We also plot a transit model corresponding to a planet with radius 0.75 R_⊕ and period 2.1076 days (this period was selected by combining the EPOXI and Spitzer 8 μm observations described in the next section).

Zoom In Zoom Out Reset image size

The constraints on the ephemeris of the putative GJ 436c from the EPOXI data alone meant that the accuracy with which we could predict the times of transits 1.5 years out from the EPOXI observations was poor. However, the extant Spitzer 8 μm phase curve was gathered only two months after the EPOXI observations took place, such that our accuracy on the predicted time was 5 hr (defined by the duration of the 2σ confidence interval from EPOXI). We performed a boxcar search of this light curve, allowing both the time of transit and the depth of transit to vary, since the EPOXI observations provided only weak constraints on the planetary radius. We discovered a transitlike signal in these data within the time window predicted from EPOXI, but the signal was present (that is, produced a Δχ² improvement larger than any other feature during the time window of interest) only when apertures smaller than 4 pixels in radius were used to perform the photometry. The top panel of Figure 5 shows a portion of the 8 μm light curve, with 3.5 pixel aperture used to extract the photometry. The solid black curve shows the best-fit transit light-curve solution, and the event at a time of BJD 2454660.4 is the secondary eclipse of the hot Neptune GJ 436b. The second panel shows the significance of the best χ² as a function of time. The bottom two panels show the results using a 7.0 pixel aperture to extract the photometry. While the significance of the secondary eclipse remains constant, the significance of the putative additional transit depends on aperture size.

Zoom In Zoom Out Reset image size

The dependence of the putative transit signal on the size of the aperture argued against the planet hypothesis; rather, a signal that disappears at larger radii is more likely due to position-dependent flux losses. We concluded that the Spitzer 8 μm observations neither definitely confirmed nor refuted the planet hypothesis, so we gathered additional Spitzer observations at 4.5 μm, where we could obtain a higher-precision light curve, to definitely resolve the question. The single candidate transit from Spitzer greatly decreased the possible parameter space of the planet's ephemeris, and so we were able to predict a transit to occur 1.5 years after the EPOXI observations within an 15 hr window. We also predicted the radius of the putative planet, from a χ² minimization of the combined EPOXI and Spitzer 8 μm observations, to be 0.75 R_⊕.

The method we use to correct the 4.5 μm data has the possible effect of suppressing transit signals, since the point-by-point correction relies critically on the assumptions that the flux variations are due only to the pixel sensitivity variations and that the stellar flux is constant. If a transit occurs during the observations, then the derived value of the pixel sensitivity W(x_i,y_i) in the location of the detector where the transiting points occur will be lower by a fraction that depends on the pixel sampling. If the points in transit comprise 10% of the observations in that location on the detector, then the value for the pixel sensitivity in that location will be low by 10%. This will have the effect of suppressing the transit depth by 10%, since we will incorrectly attribute a fraction of the decrement at that time to the pixel performance, rather than to an astrophysical variation. Therefore, to correctly recover the true transit signal, we must iteratively identify and then mask the points that occur in transit. This presents a challenge, since we do not know exactly when the transit should occur. We tested the procedure of simply masking points that occur within a transit duration of each point being corrected (so that the sensitivity function for each point is calculated using points more than half a transit duration removed in time), but found that this did not produce the highest-quality time series. Although this masking procedure prevents suppression of the transit signal when the mask is corrected located over the transit, the transiting points are allowed to contribute to the sensitivity function at other times, which introduces correlated noise to the remainder of the time series. Although the depth of the transit signal is preserved, its significance relative to other features in the time series is diminished, reducing our ability to correctly identify it. We therefore concluded that the correct procedure is to locate the position of the transit in time and to mask the set of transiting points for each individual flux correction. The challenge is to locate the true position of the transit in order to correctly place the mask. We addressed this question by producing n time series, where n is the number of tested mask positions (in this case, we tested mask positions in 30 minute intervals over a roughly 15 hr time series, resulting in 33 mask positions). We hypothesize that when the mask correctly coincides with a transit, both the depth of the transit and its significance relative to the next most significant feature in the time series should increase, thus enabling us to identify the time of transit.

In order to establish that we could reliably detect the 0.75 R_⊕ radius planet, we injected a transit of this size into the Spitzer 4.5 μm observations. We attempted to blindly recover the injected transit time by varying the mask position (which was always 1 hr in duration) by 30 minute intervals, producing a separate time series for each mask position. Then, for each of these n time series, we evaluated the significance of a boxcar light-curve function with the predicted transit depth and duration, allowing only the time of the transit to vary. We identified the time of transit from the time of most significant improvement to the χ² from the boxcar search among all the time series, which indeed occurred when the mask was correctly located over the injected transit. We then repeated this procedure, shifting the injected transit time in hour increments, to ensure that we could detect the transit at any time during the observations. Figure 6 shows the 4.5 μm light curve, with 0.75 R_⊕ transit signal injected at a location denoted by the solid line. The dashed lines mark the beginning and end of the window in time during which data are masked from the sensitivity-function calculation. In this case, the significance of the detection using the improvement over the null χ² increases by a factor of 2 when the mask is correctly located over the transit, which enabled us to blindly locate the time of the injected transit.

Zoom In Zoom Out Reset image size

We assess our sensitivity to 0.75 R_⊕ planet transits by noting both the absolute significance of the detection and the ratio of the detection significance to the significance of the next best solution for all injected signals. We find that for all transits occurring after a time of BJD 2455224.9 (corresponding to a period of 2.1071 days, which is safely before any ephemeris, consistent with the solution from EPOXI), we recover the correct transit time with significance Δχ² ≳ 20, and in all cases the significance of the transit signal (once the mask is centered over the transit time, so the suppression is minimized) is at least 60% higher than the significance of the next-highest solution. Figure 7 shows the recovery statistics for all injected transit times, starting at BJD 2455224.90 and ending at BJD 245525.44, in increments of 1 hr.

Zoom In Zoom Out Reset image size

We note also that the average detection significance is Δχ² of 45, as compared to a predicted significance of 37 (using the scatter of 71 ppm for 20 minute bins, compared to the transit depth of 250 ppm, and assuming an hour-long transit). Having demonstrated our sensitivity to transits as small as the putative GJ 436c, we then repeated the preceding analysis on the actual time series—generating different versions of the time series for each mask position, while keeping the duration of the mask constant at 1 hr. The solution with the best improvement over the null hypothesis is shown in Figure 8, with the beginning and end of the interval during which data are masked from the weighted-sensitivity function shown by dashed lines. The best solution is actually an antitransit in this case; when the mask is located over these points, the solution with highest significance gives an improvement over the null χ² of 15. The most significant solution with a transit decrement solution has a significance of 7. We do not find any signal with the significance at which we detected injected signals of the expected depth.

Zoom In Zoom Out Reset image size

We repeated this analysis on a version of the time series that was reduced using a polynomial fit to the intrapixel sensitivity variation. The standard deviation of the time series, as discussed in § 2.3, is 30% higher than the standard deviation of the time series reduced using the weighted-sensitivity function, even after we fit coefficients independently to 5 hr pieces of the time series. We find that we are able to recover the correct injected transit time only 20% of the time when the time series is reduced using a polynomial fit to the sensitivity function.