SPITZER OBSERVATIONS CONFIRM AND RESCUE THE HABITABLE-ZONE SUPER-EARTH K2-18b FOR FUTURE CHARACTERIZATION

Our photometric extraction routine is outlined in Crossfield et al. (2015) and Petigura et al. (2015), following the approach introduced by Vanderburg & Johnson (2014). In brief, during the continuous 80 day K2 observations, the stars drift across the CCD by approximately 1 pixel every 6 hr due to the spacecraft's pointing jitter. As the stars drift through pixel-phase, intra-pixel sensitivity variations and errors in the flatfield cause the apparent brightness of the target star to change. We detrend the apparent brightness variations of our target using the telescope roll angle between the target frame and an arbitrary reference observation. For each of the two transits, we then extract transit light curves ranging from three hours before the transit ingress to three hours after ingress, providing sufficient baseline for the transit light curve analysis.

As shown in this work, cosmic ray hits can substantially affect the astrophysical results derived from K2 photometry. This is particular troubling because the effects of cosmic ray hits are generally too small to result in obvious outliers in the long-cadence, 30 minute K2 photometry. As a result, cosmic ray hits have stayed unnoticed to date, and at least for K2-18b have resulted in inaccurate estimates of the transit parameters. To address this issue and detect cosmic ray hits in the K2 photometry, we introduce a generally applicable algorithm to efficiently identify cosmic ray hits in the presence of substantial telescope jitter as we see for K2.

We apply the algorithm to the K2-18 data as follows. First, we remove the sky background from each individual frame by subtracting the median pixel value outside the point-spread function (PSF) of the target star. For each background-subtracted frame in the K2 image series, we then find a "most similar frame" by identifying the frame in the K2 image series that minimizes the weighted sum of the square differences in x and y centroid location and telescope roll angle. This most similar frame will generally appear virtually identical to the original frame because the PSF falls virtually identical on the detector pixels. The most similar frame can, therefore, be used to subtract the target star from the image and isolate potential cosmic ray hits. Finally. we automatically identify cosmic ray hits by searching for $\gt 10\sigma$ outliers in the time series of each pixel in the difference images.

As an example, Figure 4 shows the detection of the cosmic ray hit near the ingress of the second K2 transit observation of K2-18b. Panel (a) shows the background subtracted image of K2-18b for the cosmic-ray affected observation near the ingress of the second transit. For comparison, panel (b) shows the frame in the K2-18 data with the most similar centroid position (BJD = 2456875.66). Panels (a) and (b) appear virtually identical despite the strong color stretching. The cosmic ray hit near x = 4.5 and y = 5.5 becomes apparent, however, in the difference image (panel (c)). The five pixels near the cosmic ray hit are identified as outliers with a significance of $\gt 20\sigma$ in the difference image light curves for these pixels. For comparison, panel (d) shows the difference image for a regular frame without a cosmic ray hit. Small residuals (∼1%) in the difference image remain due to slight differences in the centroid position and shape of the PSF. However, the residuals are not mistaken as cosmic ray hits because the pixel values are within the variances of the difference image light curves for those pixels.

Zoom In Zoom Out Reset image size

Our Spitzer/IRAC instrument model accounts for intrapixel sensitivity variations and temporal sensitivity changes using a modified version of the systematics model proposed by Deming et al. (2015). Our instrument model is

$$\begin{eqnarray}&&S({t}_{i})=\displaystyle \frac{{\displaystyle \sum }_{k=1}^{9}{w}_{k}{D}_{k}({t}_{i})}{{\sum }_{k=1}^{9}{D}_{k}({t}_{i})}+m\cdot {t}_{i},\end{eqnarray} \tag{ 1 }$$

where the sensitivity function $S({t}_{i})$ is composed of the pixel-level decorrelation (PLD) term introduced by Deming et al. (2015) and a linear sensitivity gradient in time. The ${D}_{k}({t}_{i})$'s in the PLD term are the raw counts in the 3 × 3 pixels, $k=1\ldots 9$, covering the central region of the PSF. In the numerator, these raw data values are multiplied by the nine time-independent PLD weights, $\{{w}_{1}\ldots {w}_{9}\}$, fitted as free parameters in the light curve analysis. Together with the linear slope m, the instrument model therefore includes 10 free instrument fitting parameters to capture the intrapixel sensitivity variations and temporal sensitivity changes. The differences between Equation (1) in this work and Equation (4) in Deming et al. (2015) are that we do not include the offset constant (h) and we apply $S({t}_{i})$ as a multiplicative correction factor corresponding to a variation in sensitivity rather than an additive term. The log-likelihood function for fitting the Spitzer raw photometry is then

$$\begin{eqnarray}&&\mathrm{log}L=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{D({t}_{i})-S({t}_{i})\cdot f({t}_{i})}{\sigma }\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where $f({t}_{i})$ is the median-normalized flux summed over all pixels of the target's PSF, $S({t}_{i})$ is aforementioned instrument sensitivity, $f({t}_{i})$ is the model transit light curve, and σ is the photometric scatter parameter simultaneously fit with the instrument model and transit light curve parameters.

We do not include the constant term h from Deming et al. (2015) because only the relative sizes of the PLD weights carry information about the intrapixel sensitivity variation. The sum of the PLD weights ${\sum }_{k}{w}_{k}$ uniformly scales the entire light curve up or down, which is perfectly equivalent to adding a constant h. As a result, if an extra term h was included, it would be perfectly degenerate with ${\sum }_{k}{w}_{k}$ in the fitting, resulting in 100% degeneracies between the nine PLD weight and h.

We choose to include $S({t}_{i})$ as a multiplicative correction factor rather than an additive term as done by Deming et al. (2015) because the multiplicative factor matches more closely the underlying detector behavior, which is a variation in sensitivity (or quantum efficiency) across the pixel area. The difference between a multiplicative factor and an additive term is generally small because $S({t}_{i})$ is near unity and the multiplication can be approximated by $1+\epsilon$. However, inaccuracy due to the negligence of the cross-term $\delta f\cdot \delta S$ can be as high as $0.01\cdot 0.005=50\,{\rm{parts \mbox{-} per \mbox{-} million}}$ for 1% transit depth and typical sensitivity variations. We choose to be on the safe side by introducing a multiplicative correction factor that correctly captures the cross-term $\delta f\cdot \delta S$.

We compute the transit light curve $f({t}_{i})$ using the Batman implementation (Kreidberg 2015) of the equations derived in Mandel & Agol (2002). The transit parameters fitted in our Spitzer analysis are the planet-to-star radius ratio ${R}_{P}/{R}_{* }$, the mid-transit time T_C, and the impact parameter b. We fix the stellar radius and mass to the values derived based on stellar spectroscopy (Table 2) and assume a circular orbit. We use a quadratic limb darkening profile with coefficients interpolated from the tables provided by Claret & Bloemen (2011) for K2-18's stellar effective temperature and surface gravity (${u}_{1}=0.007\pm 0.007$, ${u}_{2}=-0.191\pm 0.005$). These limb darkening coefficients are computed specifically for the $4.5\,\mu {\rm{m}}$ Spitzer/IRAC Channel 2 bandpass from spherically symmetric Phoenix models (e.g., Hauschildt et al. 1999) using updated opacities. We account for the ∼30 minute cadence of the K2 observations by numerically integrating in time.

We compute the joint posterior distribution of the instrument and transit parameters using the emcee package (Foreman-Mackey et al. 2013), a Python implementation of the Affine Invariant MCMC Ensemble sampler (Goodman & Weare 2010). We seed 60 MCMC walkers with initial values widely spread in the prior parameter space. For convergence, we ensure that the chains for all parameters are well-mixed as indicated by Gelman–Rubin metrics smaller than 1.02 (Gelman & Rubin 1992). After an initial burn-in phase, we generally find good convergence after 3000–4000 iterations for each of the 60 walkers. Since the computational time is not a limiting factor in this work, we quadruple the number of iterations to obtain smooth posterior distributions (Figure 6). The final confidence intervals reported in this work are the 15.87% and 84.13% percentiles of each parameters' posterior distribution.

Zoom In Zoom Out Reset image size

Our Spitzer observations robustly reveal a transit event consistent in transit depth and duration with the super-Earth candidate K2-18b. The transit event, however, occurs approximately 1.85 hr ($7-\sigma$) before the transit time predicted for K2-18b by Montet et al. (2015), which will be discussed further in the following section.

We also find that, using our new modified PLD systematics model, the final systematics-corrected photometry is near the photon-noise limit and virtually free of red noise and systematics (Figures 2 and 3). All posteriors of Spitzer/IRAC systematics parameters are Gaussian-shaped and uncorrelated with the astrophysical parameters in the transit model, indicating that there is no dependency of the astrophysical parameters on the instrument parameters (Figure 6).

Since we fit the Spitzer photometry and two K2 transits simultaneously, our fit now includes 14 free parameters for instrument systematics. As in Section 3.1, we include nine PLD weights and one linear slope to account for the intrapixel sensitivity variations and temporal sensitivity changes in the Spitzer photometry (Section 3.1.1). In addition, two free parameters are included for each of the two K2 transits to account for the linear slope and offset in the K2 baseline. The log-likelihood function for simultaneously fitting the Spitzer transit and the two K2 transits is

$$\begin{eqnarray}\mathrm{log}L & = & -\,\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{{N}_{{\rm{S}}}}{\left(\displaystyle \frac{{D}_{\mathrm{Spitzer}}({t}_{i})-S({t}_{i})\cdot f({t}_{i})}{{\sigma }_{\mathrm{Spitzer}}}\right)}^{2}\\ & & -\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{{N}_{1}}{\left(\displaystyle \frac{{D}_{{\rm{K}}\mathrm{2,1}}({t}_{i})-({m}_{1}{t}_{i}+{b}_{1})\cdot f({t}_{i})}{{\sigma }_{{\rm{K}}2}}\right)}^{2}\\ & & -\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{{N}_{2}}{\left(\displaystyle \frac{{D}_{{\rm{K}}\mathrm{2,2}}({t}_{i})-({m}_{2}{t}_{i}+{b}_{2})\cdot f({t}_{i})}{{\sigma }_{{\rm{K}}2}}\right)}^{2},\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{\mathrm{Spitzer}}$ and ${\sigma }_{{\rm{K}}2}$ are the fitted photometric scatter values for the Spitzer and K2 photometry, N_s, N₁, and N₂ are the number of data points for each of the transit light curves, and ${{mt}}_{i}+b$ is the linear baseline for the individual K2 transit light curves.

Our global light curve fit includes a wavelength-independent planet-to-star radius ratio $p={R}_{P}/{R}_{* }$, the impact parameter b, the mid-transit time T_C, and the orbital period P. As in Section 3.1.2, we fix the stellar radius and stellar mass to the values derived based on stellar spectroscopy (Table 2) and assume a circular orbit. For the $4.5\,\mu {\rm{m}}$ Spitzer/IRAC Channel 2 observations, we use the same quadratic limb darkening profile (${u}_{1}\ =0.007\pm 0.007$, ${u}_{2}=-0.191\pm 0.005$) as in Section 3.1.2. For the visible Kepler bandpass, we derive the quadratic limb darkening coefficients ourselves using Phoenix models and obtain ${u}_{1}=0.153\pm 0.004$ and ${u}_{2}=0.261\pm 0.007$.

Our joint analysis of the Spitzer and K2 data reveals that a previously undetected outlier point in the K2 photometry near the ingress of second K2 transit is the reason for the 1.85 hr ($7\sigma$) discrepancy between the transit time predicted from the K2 data and the transit time observed by Spitzer (Figure 7). Using our newly developed cosmic ray detection algorithm we find that the outlier is caused by a cosmic ray hit near the edge of the target star's PSF on the detector (Figure 4). After removal of the outlier point, we find that a Keplerian orbit with periodic transit events provides a good joint fit to the K2 and Spitzer transit light curves (Figure 7). Our new best estimate for the transit ephemeris of K2-18b is ${T}_{0}={2457264.39131}_{-0.00067}^{+0.00060}$ and $P={32.939614}_{-0.000084}^{+0.000101}$.

Zoom In Zoom Out Reset image size

As for the Spitzer-only fit in Section 3.1.4, we find that the posteriors of all transit light parameters are Gaussian-shaped and uncorrelated with the parameters in the systematics model. This indicates that the derived astrophysical parameters are independent of the fitted instrument parameters in our joint Spitzer/K2 fit. We further find that including the K2 in the fit does not compromise the excellent noise characteristic of the fit to the Spitzer data (compare Figures 3 and 7(right panel)).