Spitzer Transit Follow-up of Planet Candidates from the K2 Mission

The basis for this work is the initial identification of planet candidates in K2 light curves. This process is described in Cr16 and Petigura et al. (2018), but we briefly summarize it here. We use k2phot²⁰ to correct the instrumental systematics induced by the roll of the Kepler spacecraft. The resulting corrected light curves are publicly available on the community portal ExoFOP.²¹ We use TERRA²² to search these light curves for transit signals, and the resulting candidates are then vetted by eye to eliminate instrumental or astrophysical false positives (FPs). During this process, we also assess the utility of conducting follow-up transit observations with Spitzer. The planets we analyze here were all deemed interesting targets for Spitzer because they are relatively small, are temperate, and/or orbit late-type host stars. All of the planets and candidates in this work were previously published by Cr16, with the exception of K2-124 b, which was observed in K2 Campaign 5 and subsequently discovered by our team (Dressing et al. 2017; Livingston et al. 2018a; Petigura et al. 2018).

Spitzer presents several advantages over ground-based transit follow-up observations: its position in space enables precise photometry unaffected by Earth's atmosphere, its Earth-trailing orbit frees it from the scheduling constraints imposed by the day/night cycle on Earth, the diminished effects of limb darkening in the infrared enable precise estimation of transit model parameters, and the 4.5 $\mu {\rm{m}}$ Infrared Array Camera (IRAC; Fazio et al. 1998) bandpass (in conjunction with the Kepler bandpass) provides a relatively broad wavelength baseline that facilitates planet validation. In addition, our high-cadence Spitzer observations provide better sampling of the transit than the 30-minute cadence of K2.

We conducted our transit observations using the IRAC 4.5 $\mu {\rm{m}}$ channel as part of Spitzer cycle 11 GO program 11026 (PI: Werner). We chose integration times between 2 and 30 s to keep the detector in the linear regime and minimize downlink bandwidth. Target acquisition places the stars on the "sweet spot" of the detector, which has been well characterized for the purpose of precise time-series photometry (Ingalls et al. 2012) and falls within the region of the detector accessible in subarray mode (used for observing bright stars). Following the guidelines for high-precision Spitzer photometry (Grillmair et al. 2012), we performed ∼30 minutes of integrations on an empty field before each transit observation, which can help mitigate systematics induced by thermal settling of the spacecraft. See Table 1 for details of the observations.

High-resolution imaging is important for detecting stellar companions and constraining the probability of chance alignments with background sources within the K2 and Spitzer photometric apertures and thus plays a critical role in assessing the false-positive probability (FPP) of a planet candidate. In this work, we utilize imaging previously published by Cr16 along with additional AO imaging from Keck/NIRC2 (Wizinowich et al. 2014). Cr16 did not obtain imaging of EPIC 205084841, so we used NIRC2 in natural guide star mode to observe the star in K band on UT 2017 July 9. Using the image reduction and analysis methods described in Cr16, we find the star to be single and rule out companions above the contrast curve shown in Figure 1 at the 5σ confidence level. We also utilize a J-band image of K2-289, which was obtained with NIRC2 on UT 2015 April 1 and made available on ExoFOP but was not published in Cr16. As in the K -band image reported by Cr16, the companion is clearly detected in J band, thus providing useful color information (see Section 4.1).

Zoom In Zoom Out Reset image size

Our team has primarily used Keck/HIRES (Vogt et al. 1994) to obtain high-resolution spectra of candidate host stars, enabling the measurement of more robust planet properties, as well as detecting (or ruling out) double-lined spectroscopic binaries (see Cr16 and Petigura et al. 2018 for more details). We use SpecMatch-syn to measure precise stellar parameters from the spectra of stars hotter than ∼4200 K, which matches spectra to an interpolated grid of models from Coelho et al. (2005). For cooler stars, we use SpecMatch-emp, which matches spectra to a spectral library of 404 standard stars (Yee et al. 2017a).

We also conducted spectroscopic observations of K2-289 using the High Dispersion Spectrograph (HDS) mounted on the Subaru 8.2 m telescope between UT 2016 April 29 and May 2. We employed the standard I2a setup, which covers 4940–7590 Å, and image slicer #2, achieving a spectral resolution of R ∼ 80,000. For the RV measurements, we used the iodine (I₂) cell on three consecutive nights. We also obtained the stellar spectrum without the I₂ cell for the template of RV measurements. The HDS data were reduced using standard IRAF routines, by which we extracted one-dimensional, wavelength-calibrated spectra of K2-289. We then put those spectra through an RV analysis pipeline (Sato et al. 2002, 2012), which does a forward modeling of each observed spectrum to measure the RV relative to the template. Table 2 lists the extracted relative RVs and their internal errors.

As part of our team's large-scale transit search of the K2 data (Section 2.1), we correct systematics induced by the coupling of spacecraft motion and intrapixel gain variations with k2phot, which uses a Gaussian Process model (GP; Rasmussen & Williams 2005). The resulting light curves are essentially free of the large-amplitude position-dependent flux variations characteristic of raw K2 photometry, but stellar variability and potentially residual systematic trends must be accounted for to measure precise transit properties. A common approach is to first model and remove out-of-transit variability using various methods, such as a median filter, polynomial, or spline model. These approaches are simple and fast and usually do not significantly impact the results for stars with low levels of variability and residual systematics. In order to minimize the potential for biased parameter estimates while employing a uniform framework for a range of light-curve behaviors, we use the celerite GP framework with a Matern-3/2 kernel and take an approach similar to the "type II" maximum likelihood estimation (MLE) described in Gibson et al. (2012).

The GP is first trained on the out-of-transit light curve using scipy.optimize, and then the full light curve is analyzed with the same GP in conjunction with the transit model (Section 3.4). We first use the "L-BFGS-B" method (Byrd et al. 1995; Zhu et al. 1997) to find the MLE value of the kernel hyperparameters given the out-of-transit light curve, using the GP likelihood and gradient in celerite. During this stage, we perform iterative outlier rejection to minimize the possibility of biased kernel hyperparameters. Next, we use the Nelder–Mead simplex algorithm (Nelder & Mead 1965) to fit the joint GP and transit model to the full light curve, initialized with the MLE kernel hyperparameters and an initial set of transit parameters from previous analyses. Fitting the GP hyperparameters simultaneously with the transit model helps ensure that we find an optimal noise model that is valid during and out of transit. We then restrict the data to 1-day windows centered on each transit and sample the posterior of the joint GP-transit model, by running emcee for 500 steps initialized with the optimum found in the previous step. This brief sampling stage is especially important if the kernel hyperparameters are initially stuck in a local optimum, and it can be used to ensure that the posteriors are unimodal and sharply peaked around these optimal values.

Finally, we fix the hyperparameters to their optima and proceed to run the sampler for 10,000 steps (see also Section 3.4). This approach retains the benefits of a flexible noise model while minimizing computational complexity. Figures 2 and 3 show examples of these fits, with data points shown in gray if they were excluded during the iterative outlier rejection performed during the initial GP training stage.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We extract the Spitzer light curves following the approach taken by Knutson et al. (2012) and Beichman et al. (2016). In brief, we compute aperture photometry using circular apertures centered on the host star, for a range of radii between 2.0 and 5.0 pixels, corresponding to 24–60 owing to Spitzer's 12 pixel scale. We used a step size of 0.1 pixels from 2.0 to 3.0 and a step size of 0.5 from 3.0 and 5.0. An important component of precision photometry with Spitzer is the selection of an optimal aperture, due to the fact that significant levels of "red" (correlated) noise are present to a varying extent in each time series. Smaller radii tend to have less photon noise owing to the decreased sky background in the aperture, while slightly larger radii can sometimes mitigate the inter- and intrapixel gain variations that are responsible for the correlated noise. Ideally, an optimal radius minimizes both of these effects, although in practice it is common to attempt only to minimize correlated noise or the photon noise (e.g., Knutson et al. 2012; Lewis et al. 2013; Lanotte et al. 2014).

A typical transit data set contains significant time both in and out of transit, so we compute relevant noise metrics as a function of radius for different subsets, most of which are fully out of transit or between second and third contact (i.e., do not contain ingress or egress). Thus, for a given radius, the ensemble of these values is largely unaffected by the transit signal and thus reflects only photon noise and systematic noise. To quantify the level of red noise, we compute β, the factor by which the standard deviation of the observed binned residuals deviates from the theoretical value (Pont et al. 2006; Winn et al. 2008):

$\begin{eqnarray}&&\beta =\displaystyle \frac{{\sigma }_{M}}{{\sigma }_{0}}\sqrt{\displaystyle \frac{N(M-1)}{M}},\end{eqnarray} \tag{ 1 }$

where σ_M is the standard deviation of the binned residuals (in M bins), σ₀ is the standard deviation of the unbinned residuals, and N is the number of data points per bin. To ensure a robust estimate and focus on the timescales of red noise that could significantly impact transit parameter estimates, we compute the median β value for bin widths between 5 and 40 minutes. We divide each flux time series into 10 equal-sized segments and compute both the standard deviation (i.e., overall noise level) and the β value (i.e., red-noise level) for each segment. The optimal aperture is then the one that minimizes each metric. We then compute the median of the optimal aperture radii for each metric over all segments. Finally, the aperture radius adopted for subsequent analysis is the mean of these two "optimal" radii; the selected aperture thus only "approximately" minimizes both metrics in cases where these two radii are not equal. We chose 10 segments as a trade-off between having more robust statistics and having enough light curve in each segment to compute red noise on a range of different timescales. We find that the minimum red-noise aperture is frequently consistent with the minimum standard deviation aperture to within a few tenths of a pixel, and the optimal radius is typically 2.2–2.4 pixels, which is consistent with the optimal apertures found in previous analyses of Spitzer transit data (e.g., Knutson et al. 2012), in which the residuals computed from the best-fit transit and systematics model are analyzed instead of the raw light curve.

In principle, aperture selection could be handled via Bayesian model selection (i.e., computing the Bayesian evidence), although any improvements might not be significant enough to justify the computational cost. Additionally, it may be fruitful to simultaneously estimate the white- and red-noise levels using a GP; a sufficient choice of kernel could more fully disentangle these two noise signals, as compared to the standard deviation and β factor. We leave an investigation of these possibilities for a future work.

We model the systematics inherent to the Spitzer light curves using the pixel-level decorrelation (PLD) method (Deming et al. 2015). In comparisons between various methods used to correct Spitzer systematics (Ingalls et al. 2016), PLD was among the top performers, displaying both high precision and repeatability. PLD uses a linear combination of (normalized) pixel light curves to model the effect of point-spread function (PSF) motion on the detector coupled with intrapixel gain variations; thus, it does not require the calculation of centroids. The parameterization of the full model for the transit light curve, including PLD, is

$\begin{eqnarray}&&{\rm{\Delta }}{S}^{t}=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{9}{c}_{i}{P}_{i}^{t}}{{\displaystyle \sum }_{i=1}^{9}{P}_{i}^{t}}+{M}_{\mathrm{tr}}({\boldsymbol{\theta }},t)+\varepsilon (\sigma ),\end{eqnarray} \tag{ 2 }$

where M_tr is the transit model, c_i are the PLD coefficients, and ε(σ) are zero-mean Gaussian errors with width σ. To form a valid set of basis vectors for the instrumental systematics component of the light curve, the astrophysical signal in each individual pixel light curve is removed by normalization (the sum in the denominator of Equation (2)). We show an illustrative example of these normalized pixel light curves in Figure 4. In testing, we found that using a 3 × 3 pixel grid sufficiently captures the information content corresponding to the motion of the PSF on the detector (which is typically ≲ a few tenths of a pixel). However, Spitzer target acquisition occasionally misses the "sweet spot" pixel, which may yield data sets with more pronounced systematics that could benefit from using a larger pixel grid (such as 5 × 5). In Figure 5 we plot the full model fit for two data sets (top panels), as well as the data corrected by subtracting the best-fit PLD noise model (bottom panels). In Figure 6 we plot the corrected data and transit models for the remaining Spitzer data sets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To model the transits in both the K2 and Spitzer light curves, we use the analytic model of Mandel & Agol (2002), assuming a circular orbit and a quadratic limb-darkening law, as implemented in batman. The free parameters ${\boldsymbol{\theta }}$ of the transit model are the planet-to-star radius ratio R_p/R_⋆, scaled semimajor axis a/R_⋆, midtransit time T₀, orbital period P, impact parameter b ≡ a cos(i)/R_⋆, and modified quadratic limb-darkening coefficients q₁ and q₂, which efficiently sample the space of physically allowed limb-darkening coefficients using the triangular sampling method of Kipping (2013). These transformed coefficients are computed directly from the quadratic limb-darkening coefficients u₁ and u₂ using Equations (17) and (18) of Kipping (2013), reproduced here for convenience:

$\begin{eqnarray}&&{q}_{1}\equiv {\left({u}_{1}+{u}_{2}\right)}^{2},\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{q}_{2}\equiv \displaystyle \frac{{u}_{1}}{2({u}_{1}+{u}_{2})}.\end{eqnarray} \tag{ 4 }$

Following Cr16, we use Gaussian limb-darkening priors in our light-curve analysis, which we derive from the coefficients for the Kepler and IRAC bandpasses tabulated by Claret et al. (2012). We use a Monte Carlo approach in which we sample the stellar parameters (${T}_{\mathrm{eff}}$, $\mathrm{log}g$, [Fe/H]) of each star and interpolate the tabulated coefficients at the sampled stellar values. Because we sample in q-space, we convert the tabulated values of u₁/u₂ to q₁/q₂ to enable using a Gaussian prior in q-space. We then use the mean and standard deviation of the sampled coefficients to define the Gaussian priors (see the Appendix for more details). For the final set of parameter estimates listed in Table 3 we also impose a prior on the mean stellar density determined from our isochrones analysis (see Section 3.7), which yields more precise parameter estimates by leveraging more information about the host star. However, we also perform a parallel set of identical analyses without this prior, which provides the opportunity to compare the density from our stellar characterization to the independent measurement of the mean stellar density from the light curve ρ_⋆,LC (see Equation (7)). We discuss this in detail in Section 4.1.

Table 3.  Transit Parameter Estimates

Name	a	b	i	Rp,K	Rp,S	Rp,C	log(σK)	log(σS)	T14	T23	η	Rp,max
	(${R}_{\star }$)		(deg)	(${R}_{\star }$)	(${R}_{\star }$)	(${R}_{\star }$)			(days)	(days)		(${R}_{\star }$)
K2-52 b	${5.2}_{-0.3}^{+0.2}$	${0.32}_{-0.18}^{+0.14}$	${86.46}_{-1.91}^{+2.04}$	${0.072}_{-0.001}^{+0.001}$	${0.078}_{-0.002}^{+0.002}$	${0.075}_{-0.001}^{+0.001}$	−7.01	−5.70	${0.222}_{-0.003}^{+0.003}$	${0.187}_{-0.003}^{+0.002}$	${0.84}_{-0.02}^{+0.01}$	${0.085}_{-0.008}^{+0.013}$
K2-53 b	${26.3}_{-0.4}^{+0.4}$	${0.31}_{-0.10}^{+0.07}$	${89.32}_{-0.16}^{+0.22}$	${0.028}_{-0.001}^{+0.001}$	${0.031}_{-0.002}^{+0.002}$	${0.030}_{-0.001}^{+0.001}$	−8.83	−6.61	${0.145}_{-0.003}^{+0.003}$	${0.136}_{-0.003}^{+0.004}$	${0.94}_{-0.00}^{+0.00}$	${0.033}_{-0.002}^{+0.002}$
205084841.01	${17.8}_{-0.6}^{+0.5}$	${0.62}_{-0.03}^{+0.03}$	${87.99}_{-0.17}^{+0.16}$	${0.127}_{-0.002}^{+0.002}$	${0.117}_{-0.002}^{+0.002}$	${0.122}_{-0.002}^{+0.002}$	−6.36	−5.25	${0.189}_{-0.002}^{+0.003}$	${0.125}_{-0.003}^{+0.003}$	${0.66}_{-0.02}^{+0.02}$	${0.203}_{-0.014}^{+0.016}$
K2-289 b	${22.5}_{-0.3}^{+0.3}$	${0.60}_{-0.02}^{+0.02}$	${88.48}_{-0.07}^{+0.07}$	${0.081}_{-0.001}^{+0.001}$	${0.082}_{-0.002}^{+0.001}$	${0.081}_{-0.001}^{+0.001}$	−7.42	−6.24	${0.168}_{-0.002}^{+0.002}$	${0.130}_{-0.002}^{+0.002}$	${0.77}_{-0.01}^{+0.01}$	${0.128}_{-0.005}^{+0.005}$
K2-174 b	${26.6}_{-0.2}^{+0.2}$	${0.02}_{-0.02}^{+0.03}$	${89.95}_{-0.05}^{+0.04}$	${0.033}_{-0.000}^{+0.000}$	${0.038}_{-0.001}^{+0.001}$	${0.036}_{-0.001}^{+0.001}$	−9.25	−6.80	${0.242}_{-0.002}^{+0.002}$	${0.225}_{-0.002}^{+0.002}$	${0.93}_{-0.00}^{+0.00}$	${0.036}_{-0.001}^{+0.001}$
K2-87 b	${14.4}_{-0.3}^{+0.3}$	${0.57}_{-0.03}^{+0.03}$	${87.74}_{-0.17}^{+0.17}$	${0.045}_{-0.001}^{+0.001}$	${0.055}_{-0.002}^{+0.002}$	${0.050}_{-0.001}^{+0.001}$	−7.84	−5.94	${0.190}_{-0.002}^{+0.002}$	${0.164}_{-0.003}^{+0.003}$	${0.86}_{-0.01}^{+0.01}$	${0.074}_{-0.004}^{+0.005}$
K2-90 b	${33.4}_{-0.3}^{+0.3}$	${0.44}_{-0.05}^{+0.05}$	${89.24}_{-0.09}^{+0.09}$	${0.036}_{-0.001}^{+0.001}$	${0.041}_{-0.002}^{+0.002}$	${0.038}_{-0.001}^{+0.001}$	−8.33	−6.37	${0.123}_{-0.003}^{+0.003}$	${0.112}_{-0.003}^{+0.003}$	${0.91}_{-0.01}^{+0.01}$	${0.048}_{-0.003}^{+0.004}$
K2-124 b	${28.0}_{-1.7}^{+1.5}$	${0.47}_{-0.14}^{+0.10}$	${89.03}_{-0.27}^{+0.32}$	${0.069}_{-0.002}^{+0.002}$	${0.068}_{-0.003}^{+0.003}$	${0.069}_{-0.002}^{+0.002}$	−6.89	−5.73	${0.070}_{-0.002}^{+0.002}$	${0.058}_{-0.003}^{+0.003}$	${0.84}_{-0.03}^{+0.02}$	${0.089}_{-0.013}^{+0.016}$

Note. Transit parameters derived from the simultaneous analysis of the K2 and Spitzer data. The parameters R_p,K, R_p,S, and R_p,C are the K2, Spitzer, and combined values of the planet radius in units of the stellar radius. The best-fit logarithms of the Gaussian errors for the K2 and Spitzer time series are denoted by log(σ_K) and log(σ_S), respectively. The parameter η is the transit shape defined in Equation (5). The maximum planet radius allowed by the shape of the transit is denoted R_p,max and defined in Equation (6).

Download table as:  ASCIITypeset image

For Markov chain Monte Carlo (MCMC) parameter estimation, we use emcee, a Python implementation of the affine-invariant ensemble sampler (Goodman & Weare 2010). We performed an initial optimization with lmfit and positioned 128 "walkers" in a Gaussian ball centered on the optimum. We then ran the sampler for 10,000 steps, allowing it to evolve according to the MCMC. Finally, we checked for convergence by visual inspection of the trace, discarding the first 5000 samples as "burn-in," and computed the autocorrelation time of each parameter using the Python package acor,²³ to ensure we had collected a sufficient number of independent samples after burn-in.

In this work we simultaneously model the K2 and Spitzer light curves of each target. This is motivated by several factors. First, the posterior distributions of some transit model parameters are often distinctly non-Gaussian (e.g., a/R_⋆), so simply imposing Gaussian priors derived from previous analysis of the K2 light curves could bias the resulting fits to the Spitzer data. Second, we leverage the K2 data to model the Spitzer data, because we simultaneously fit the systematics and transit models. Some transit model parameters (a/R_⋆ and b) are shared between Spitzer and K2, while others are distinct. We fit R_p/R_⋆ separately for the Kepler and Spitzer bandpasses to enable the detection of FP scenarios (see, e.g., Figure 7).

Our procedure is as follows. We first fit the K2 data alone and form Gaussian priors on T₀ and P based on the mean and standard deviation of the posterior distributions from MCMC. These priors are then used in a fit to the Spitzer transit data alone, yielding an initial set of parameter estimates derived from the Spitzer data. We then fit the K2 and Spitzer data simultaneously, without priors on T₀ and P. This simultaneous fit is responsible for the improvement in the ephemeris estimates, for reasons discussed above, as well as more robust transit shape parameter estimates. We list the transit parameters in Table 3 and the ephemeris estimates in Table 4.

Table 4.  Ephemerides

	K2-only		K2 + Spitzer		Mid-2021 Timing Uncertainty
Name	P	T0	P	T0	K2-only	K2 + Spitzer	Improvementa
	(days)	(BKJD)	(days)	(BKJD)	(hr)	(minutes)
K2-52 b	${3.534890}_{-0.000181}^{+0.000187}$	${2063.02742}_{-0.00184}^{+0.00190}$	${3.535055}_{-0.000016}^{+0.000017}$	${2063.02618}_{-0.00111}^{+0.00108}$	3.1	17	11×
K2-53 b	${12.207253}_{-0.000971}^{+0.000940}$	${2063.38751}_{-0.00373}^{+0.00400}$	${12.207720}_{-0.000113}^{+0.000123}$	${2063.38582}_{-0.00186}^{+0.00189}$	4.7	35	8×
205084841.01	${11.310159}_{-0.000557}^{+0.000564}$	${2060.85877}_{-0.00221}^{+0.00224}$	${11.310099}_{-0.000045}^{+0.000045}$	${2060.85913}_{-0.00109}^{+0.00110}$	3.0	14	12×
K2-289 b	${13.157344}_{-0.000634}^{+0.000622}$	${2064.14222}_{-0.00198}^{+0.00205}$	${13.156969}_{-0.000051}^{+0.000047}$	${2064.14325}_{-0.00094}^{+0.00097}$	2.9	13	13×
K2-174 b	${19.564172}_{-0.000966}^{+0.000996}$	${2250.77803}_{-0.00099}^{+0.00100}$	${19.562307}_{-0.000076}^{+0.000078}$	${2250.77917}_{-0.00094}^{+0.00095}$	2.8	13	13×
K2-87 b	${9.726793}_{-0.000646}^{+0.000663}$	${2239.30189}_{-0.00242}^{+0.00229}$	${9.726618}_{-0.000055}^{+0.000055}$	${2239.30232}_{-0.00171}^{+0.00171}$	3.8	19	12×
K2-90 b	${13.733225}_{-0.001187}^{+0.001206}$	${2245.66111}_{-0.00190}^{+0.00185}$	${13.733314}_{-0.000099}^{+0.000083}$	${2245.66112}_{-0.00172}^{+0.00170}$	4.9	22	13×
K2-124 b	${6.413721}_{-0.000272}^{+0.000271}$	${2309.18063}_{-0.00182}^{+0.00187}$	${6.413651}_{-0.000031}^{+0.000036}$	${2309.18106}_{-0.00107}^{+0.00105}$	2.3	17	8×

Notes. BKJD is the Barycentric Julian Date offset by the beginning of the Kepler mission, i.e., BJD –2,454,833.

^aRelative improvement in period measurement precision from our joint analysis of the K2 and Spitzer data, as compared to K2-only.

Download table as:  ASCIITypeset image

An additional benefit to this approach is that it allows the high cadence and diminished limb darkening of our Spitzer light curves to yield improved constraints on transit parameters that are sensitive to ingress and egress (e.g., the impact parameter b). For example, even though our Spitzer observation (just barely) missed the ingress of K2-289 b (see Figure 5), the joint analysis of the K2 and Spitzer data yields better constraints on the transit geometry than either the Spitzer-only or K2-only analyses. Figure 8 shows several key posterior distributions for this system from K2-only, Spitzer-only, and our joint analysis of both data sets, illustrating the improved constraints in transit geometry. The leftmost panels of Figure 8 show the stellar density from our isochrones analysis as gray bands for comparison to the mean stellar density derived from Equation (7) and the transit fit posteriors (without a density prior). In this case the density estimated from the K2 data alone can be seen to be in mild disagreement with the isochrones density, but the Spitzer data yield an improved density estimate in good agreement. The improvement in parameter estimates afforded by this simultaneous modeling approach can be thought of as a type of "Bayesian shrinkage," in which the high cadence of Spitzer and high photometric precision of K2 work together to extract higher measurement precision from the data.

The addition of even a single follow-up transit measurement can result in significant refinement to the ephemeris of K2 planet candidates, due to the relatively small number of transits detected in each ∼80-day K2 observing campaign. The ephemerides from our simultaneous fit to the K2 and Spitzer data are often a highly significant improvement over the K2-only ephemerides. To quantify this improvement, we compute the factor by which the precision on the estimate of the orbital period increases between the K2-only and joint fits to the K2 and Spitzer data, which ranges from 8× to 13× (see Table 4).

The benefits of updating the ephemeris can be especially important in the context of planning future transit observations. Large uncertainties in the ephemeris estimates from the K2 light curves are common, due to the relatively short time span of each K2 observing campaign. The problem is particularly pronounced for planets with longer periods, which may transit only a small number of times in a given campaign. For those interested in studying the atmospheres of these planets, this complicates scheduling of future transit observations, due to the need to lengthen the observation window to ensure that the transit is fully observed.

However, the addition of a single transit observation with Spitzer dramatically reduces the length of the necessary observation windows in the JWST era. First, the longer time baseline yields an order-of-magnitude improvement in the precision of our orbital period estimates. Second, the new, combined transit midtime results in smaller propagated timing uncertainties as a result of being closer in time to the planned transit observation than the K2 epoch. Our Spitzer observations reduce the typical uncertainty on the predicted transit time in mid-2021 from hours to minutes (see Table 4), enabling follow-up observations to be much more efficiently scheduled; even in 2025, the uncertainties will still be less than an hour.

We used the Python package isochrones to infer a uniform set of stellar parameters using priors from spectroscopy (when available), Two Micron All Sky Survey (2MASS) JHK photometry (Skrutskie et al. 2006), and Gaia DR2 parallaxes (Gaia Collaboration et al. 2016, 2018). We list the inputs to isochrones in Table 5. We sampled the posteriors using MultiNest (Feroz et al. 2013) in conjunction with the MIST stellar evolution models (Choi et al. 2016). The resulting posteriors are listed in Table 6.

Table 5.  Gaia DR2 Parallax, 2MASS Photometry, and Spectroscopic Priors Used to Estimate Stellar Parameters

EPIC	π	J	H	K	${T}_{\mathrm{eff}}$	$\mathrm{log}g$	$[\mathrm{Fe}/{\rm{H}}]$	Note
	(mas)	(mag)	(mag)	(mag)	(K)	(cgs)	(dex)
203776696	0.9394 ± 0.1169	12.729 ± 0.026	12.126 ± 0.021	11.853 ± 0.019	⋯	⋯	⋯	⋯
204890128	7.2173 ± 0.1161	10.306 ± 0.024	9.826 ± 0.024	9.664 ± 0.021	5278 ± 100	4.546 ± 0.100	−0.03 ± 0.06	1
205084841	1.1594 ± 0.1117	13.443 ± 0.026	12.797 ± 0.026	12.612 ± 0.027	⋯	⋯	⋯	⋯
205686202	3.6170 ± 0.1037	11.435 ± 0.021	10.847 ± 0.022	10.641 ± 0.021	5365 ± 110	⋯	0.16 ± 0.08	3
210558622	9.9783 ± 0.1062	10.231 ± 0.021	9.649 ± 0.022	9.496 ± 0.017	4310 ± 70	⋯	0.11 ± 0.09	2
210731500	2.0071 ± 0.1040	11.811 ± 0.020	11.359 ± 0.022	11.197 ± 0.019	5694 ± 100	4.067 ± 0.100	0.36 ± 0.06	1
210968143	7.4234 ± 0.1026	11.168 ± 0.023	10.495 ± 0.027	10.360 ± 0.020	4465 ± 100	4.553 ± 0.100	−0.25 ± 0.06	1
212154564	7.0953 ± 0.1198	12.838 ± 0.023	12.227 ± 0.022	11.975 ± 0.018	3443 ± 70	⋯	−0.13 ± 0.09	2

Note. (1) Keck/HIRES SpecMatch-syn; (2) Keck/HIRES SpecMatch-emp; (3) Subaru/HDS SpecMatch-emp. Following Luri et al. (2018), we added 0.1 mas in quadrature to the Gaia DR2 parallax uncertainties to account for systematic errors.

Download table as:  ASCIITypeset image

Table 6.  Stellar Parameters

EPIC	${T}_{\mathrm{eff}}$	$\mathrm{log}g$	$[\mathrm{Fe}/{\rm{H}}]$	${M}_{\star }$	${R}_{\star }$	Distance	${\rho }_{\star }$
	(K)	(cgs)	(dex)	(M⊙)	(R⊙)	(pc)	(cgs)
203776696	${7147}_{-506}^{+483}$	${3.979}_{-0.085}^{+0.087}$	${0.055}_{-0.139}^{+0.154}$	${1.691}_{-0.143}^{+0.157}$	${2.192}_{-0.232}^{+0.280}$	${1041.7}_{-104.1}^{+132.7}$	${0.22}_{-0.06}^{+0.08}$
204890128	${5263}_{-86}^{+88}$	${4.555}_{-0.028}^{+0.022}$	$-{0.037}_{-0.053}^{+0.052}$	${0.851}_{-0.038}^{+0.032}$	${0.807}_{-0.015}^{+0.016}$	${139.1}_{-2.0}^{+2.2}$	${2.29}_{-0.18}^{+0.15}$
205084841	${6245}_{-284}^{+369}$	${4.271}_{-0.066}^{+0.062}$	${0.019}_{-0.175}^{+0.150}$	${1.168}_{-0.094}^{+0.120}$	${1.313}_{-0.113}^{+0.126}$	${860.4}_{-76.4}^{+84.1}$	${0.73}_{-0.16}^{+0.18}$
205686202	${5529}_{-74}^{+77}$	${4.393}_{-0.022}^{+0.024}$	${0.162}_{-0.074}^{+0.087}$	${0.953}_{-0.051}^{+0.037}$	${1.025}_{-0.028}^{+0.030}$	${270.7}_{-7.1}^{+7.6}$	${1.24}_{-0.09}^{+0.10}$
210558622	${4455}_{-29}^{+31}$	${4.625}_{-0.017}^{+0.017}$	${0.113}_{-0.071}^{+0.065}$	${0.700}_{-0.023}^{+0.025}$	${0.676}_{-0.008}^{+0.008}$	${100.1}_{-1.0}^{+1.0}$	${3.21}_{-0.14}^{+0.15}$
210731500	${5747}_{-49}^{+58}$	${4.208}_{-0.035}^{+0.032}$	${0.328}_{-0.045}^{+0.051}$	${1.164}_{-0.044}^{+0.042}$	${1.401}_{-0.062}^{+0.083}$	${503.6}_{-21.9}^{+29.9}$	${0.59}_{-0.08}^{+0.07}$
210968143	${4484}_{-50}^{+58}$	${4.653}_{-0.014}^{+0.017}$	$-{0.225}_{-0.053}^{+0.053}$	${0.629}_{-0.018}^{+0.022}$	${0.619}_{-0.010}^{+0.009}$	${134.9}_{-1.7}^{+1.7}$	${3.74}_{-0.14}^{+0.19}$
212154564	${3570}_{-70}^{+63}$	${4.858}_{-0.088}^{+0.012}$	${0.023}_{-0.106}^{+0.061}$	${0.390}_{-0.042}^{+0.013}$	${0.388}_{-0.009}^{+0.016}$	${140.4}_{-2.3}^{+2.3}$	${9.58}_{-2.12}^{+0.47}$

Download table as:  ASCIITypeset image

The open-source Python package vespa (Morton 2015b) has been used to validate planets via statistical FPPs in numerous recent works (e.g., Montet et al. 2015; Morton et al. 2016) and is similar to previous methods developed for Kepler and CoRoT, such as BLENDER (Torres et al. 2011) and PASTIS (Díaz et al. 2014). Cr16 used vespa to compute FPPs for planet candidates from K2's first five observing campaigns, including those we analyze here (with the exception of K2-124 b). vespa uses the TRILEGAL Galaxy model (Girardi et al. 2005) to simulate populations of FP scenarios (i.e., eclipsing binaries [EBs], background EBs [BEBs], and hierarchical triple EB [HEB] systems) and then compares these to the observed phase-folded light curve to compute statistical likelihoods for each type of false FP, as well as the planetary scenario. For a more detailed description of how we use vespa, see Cr16.

In common practice among planet hunters, the "by-eye" transit shape is used as a first line of defense in the initial stages of vetting newly detected planet candidates in transit surveys: EBs are usually more obviously "V-shaped" than a planetary transit, as well as being deeper (in the absence of significant dilution). However, the 30-minute cadence of K2 can make the transits of real planets appear more V-shaped, so a more quantitative assessment based on precise parameter estimates is more reliable. The ratio of the full transit duration (T₂₃) to the total transit duration (T₁₄) is directly determined by the transit geometry, regardless of dilution from any additional sources within the photometric aperture (Seager & Mallén-Ornelas 2003). We denote this ratio η, the transit "shape"—a number that scales between 0 and 1, where 0 corresponds to a "V-shaped" transit and 1 corresponds to a "box-shaped" transit:

$\begin{eqnarray}&&\eta =\displaystyle \frac{{T}_{23}}{{T}_{14}}.\end{eqnarray} \tag{ 5 }$

The transit shape thus sets an upper limit on R_p/R_⋆:

$\begin{eqnarray}&&{R}_{p,\max }/{R}_{\star }=\displaystyle \frac{1-\eta }{1+\eta }.\end{eqnarray} \tag{ 6 }$

As $\eta \to 1$, ${R}_{p}/{R}_{\star }\to {R}_{p,\max }/{R}_{\star }$, meaning that the constraint is stronger for more "box-shaped" transits. This constraint can be particularly useful in the case of dilution from a known stellar companion (detected either photometrically or spectroscopically). vespa implicitly uses the shape information content of the phase-folded transit light curve to compute the likelihoods of FP and planet models. Thus, an independent consideration of FP scenarios is enabled by quantifying the shape η of each transit. Furthermore, a constraint on the maximum allowed dilution results directly from comparison of the observed values of R_p/R_⋆ and R_p,max/R_⋆. The transit shape η and the maximum radius ratio R_p,max/R_⋆ are listed in Table 3, along with other parameters of interest.

Another useful constraint derived from the transit fit is the estimate of the mean stellar density. We compute this estimate of the stellar density directly from the observed transit light-curve fit parameters using Equation (4) of Kipping (2014):

$\begin{eqnarray}&&{\rho }_{\star ,\mathrm{LC}}=\displaystyle \frac{3\pi {\left(a/{R}_{\star }\right)}^{3}}{{\mathrm{GP}}^{2}}.\end{eqnarray} \tag{ 7 }$

This estimate can then be compared directly to independent estimates of the mean stellar density, i.e., from spectroscopy. Significant disagreement could arise from the violation of any of the assumptions inherent to this estimate (i.e., non-negligible blending, a noncircular orbit, M_p ∼ M_⋆).

The planets in our sample with established "validated" dispositions (K2-52 b, K2-53 b, K2-87 b, K2-90 b, Cr16; K2-124 b, Dressing et al. 2017; Livingston et al. 2018a; K2-174 b, Mayo et al. 2018) do not exhibit suspiciously chromatic transit depths that would indicate that they are actually FPs. With the exception of K2-174, our transit analyses (without density priors) yield mean stellar densities (ρ_⋆,LC) within 2σ of the values from our independent assessment (Section 3.7). This agreement provides an additional layer of confidence in the disposition of these planets, as well as the quality of our transit analysis and stellar characterization.

On the other hand, K2-174 b yields ${\rho }_{\star ,\mathrm{LC}}={0.53}_{-0.21}^{+0.23}$ g cm⁻³, which is ∼10σ discrepant with the value we derive for the host star, ${\rho }_{\star }={3.21}_{-0.14}^{+0.15}$ g cm⁻³. Furthermore, the addition of a density prior does not significantly change this result, which indicates that the data strongly constrain the density to this value. We interpret this result as an indication of eccentricity in the system, which could potentially be measured via RVs (see Section 4.2).

In the next two subsections we examine the dispositions of the two previously unvalidated planet candidates in our sample, EPIC 205084841.01 and K2-289 b.

4.1.1. EPIC 205084841.01

EPIC 205084841.01 does not warrant validation at present, but neither should it be considered a "confirmed" FP; the current disposition of planet candidate is appropriate. Using the contrast curve derived from our Keck/NIRC2 AO image (see Figure 1), along with the isochrones inputs listed in Table 5, we compute an FPP of 15.4% for this target with vespa, somewhat higher than the FPP of 2.7% reported by Cr16, and significantly above the commonly used validation threshold of 1%. We discuss below several considerations pertinent to the disposition of this system.

Caution is warranted by the candidate's large radius of ∼17 ${R}_{\oplus }$, which is comparable to the radii of low-mass stars (Shporer et al. 2017). Additionally, there is moderate tension between R_p/${R}_{\star }$ in the Kepler and Spitzer bandpasses (see Figure 7). The posterior of R_p/${R}_{\star }$ is ∼3.3σ larger in the Kepler bandpass, which is consistent with an occultation by a lower-mass star. Although the FPP implies that the candidate is more likely a planet than an FP, vespa does not take into account chromaticity of transit depth. The second-most likely scenario reported by vespa is an EB, and vespa essentially rules out a hierarchical or background EB scenario.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Assuming that the candidate is actually an EB with the same orbital period, a secondary eclipse could be measurable in the K2 data; the absence of any apparent secondary eclipses (deeper than ∼0.1% at any phase; see Figure 2) is therefore suggestive of an eccentric orbit. Our transit analysis (without a density prior) yields ${\rho }_{\star ,\mathrm{LC}}={1.43}_{-0.34}^{+0.29}$ g cm⁻³, which is 1.8σ higher than the value we measure for the host star (see Table 6) and thus modestly suggestive of eccentricity. There is no apparent variation in the transit depth, and the transit geometry implies a radius ratio much smaller than unity (≲20%; see Table 3), so an EB at twice the estimated orbital period is unlikely. We conclude that EPIC 205084841.01 is potentially an FP caused by the eclipses of a low-mass star in an eccentric orbit, but further observations are required to confirm this scenario (e.g., RV monitoring); based on its FPP alone, the candidate is more likely to be an inflated gas giant planet.

4.1.2. K2-289 b

Cr16 reported an FPP of 1.3 × 10⁻¹¹ for this candidate but did not validate the candidate owing to the presence of a nearby stellar companion revealed by AO imaging. Using the contrast curve for this star in Figure 1 and the isochrones inputs listed in Table 5, we find a higher (but still rather low) FPP of 3.7 × 10⁻⁴. The companion is at a separation of 08 and is 3.8 mag fainter (in K band) than the primary star. The 12 pixel scale of IRAC is thus insufficient to resolve the companion in our follow-up Spitzer transit photometry. Intriguingly, the companion was not detected by Gaia DR2, but K2-289 is listed with zero excess astrometric noise, which is often associated with binarity (see, e.g., Evans 2018). This underscores the need for high-resolution imaging, even in the Gaia era.

From our AO imaging (Section 2.3), we compute deblended J-band magnitudes of 11.461 ± 0.024 and 15.554 ± 0.055 for the primary and secondary stars, respectively; in K band the deblended magnitudes are 10.674 ± 0.021 and 14.446 ± 0.024. The J − K color of the secondary star is thus 1.108 ± 0.060, which is consistent with a late M dwarf.

From the observed transit depth, an EB with 100% eclipse depth that is less than 5.4 mag fainter than the primary could reproduce the observed transit depth. However, the J − K color of the companion suggests that it is a late M dwarf, so the contrast in the Kepler bandpass should be ≳4 mag greater than in K band. In this likely scenario, we can rule out the secondary as the source of the signal, as the dilution from the primary star in the Kepler bandpass would require a depth larger than 100%. Furthermore, we can leverage our measurement of the transit geometry to eliminate considerations of the secondary star's spectral type. From our joint analysis of the K2 and Spitzer data, we have the constraint R_p,max/R_⋆ < 0.153 (5σ), which corresponds to a maximum undiluted transit depth of ∼2.3%. Given the observed transit depth, this translates to an upper limit on the amount of dilution of about 1 mag. The secondary star would thus need to have a nonphysical color of V − K ≲ −2.8 to be the source of the signal. We therefore conclude that the observed signal comes from the primary star.

From the above considerations we can confidently rule out scenarios in which the signal comes from the secondary source, i.e., the signal is due to a BEB or eclipsing HEB. We now consider the possibility that the primary is itself an EB. While low-mass stars can potentially be roughly Jovian in size, such massive bodies would induce a large-amplitude Doppler signal in time-series RV measurements of the primary star, as well as potentially significant secondary eclipses.

The RVs obtained with Subaru/HDS (see Section 2.4) show no significant variation on a timescale of ∼3 days, ruling out the possibility that K2-289 is an EB. Assuming a circular orbit, we used RadVel to estimate the RV semi-amplitude exerted by the transiting companion as $K={0.03}_{-0.03}^{+5.47}$ m s⁻¹, which is consistent with a null detection. The semi-amplitude best-fit value and 3σ upper limit are 11.8 and 46.1 m s⁻¹, respectively. This result is consistent with the vespa result, as well as our assessment of the K2 and Spitzer transit data, which both suggest a planetary origin for the observed transit signal, so we conclude that this is a valid planet.

Since we do not account for the dilution from the companion in our transit fits, the planet radius we measure may be underestimated. However, the delta-magnitude of 3.82 (in K band) implies a planet radius only ∼1.5% larger than we list in Table 7, which is a factor of ∼6 smaller than the precision of our planet radius measurement. Furthermore, the dilution is likely to be lower in Kepler band than K band, given that the secondary is probably a bound late-type companion. Future studies may yield precise enough measurements of R_p/R_⋆ that dilution from the companion will need to be accounted for.

Table 7.  Derived Planet Properties

Name	Rp	Mp	Kpred	a	${T}_{\mathrm{eq}}$	g	H	δTS
	(${R}_{\oplus }$)	(${M}_{\oplus }$)	(m s−1)	(au)	(K)	(g⊕ )	(km)	(ppm)
K2-52 b	18.0 ± 2.2	128.6 ± 109.9	38.1 ± 32.6	0.054 ± 0.002	2004 ± 184	0.4	2129	1049
K2-53 b	2.6 ± 0.1	9.2 ± 2.5	2.8 ± 0.8	0.098 ± 0.001	665 ± 13	1.3	211	112
205084841.01	17.5 ± 1.7	123.4 ± 96.3	31.8 ± 25.0	0.104 ± 0.003	979 ± 73	0.4	1029	1374
K2-289 b	9.1 ± 0.3	48.6 ± 18.6	13.6 ± 5.2	0.107 ± 0.001	753 ± 16	0.6	545	623
K2-174 b	2.6 ± 0.1	9.2 ± 2.5	2.8 ± 0.8	0.126 ± 0.001	455 ± 5	1.3	145	110
K2-87 b	7.6 ± 0.4	37.9 ± 12.7	10.3 ± 3.4	0.094 ± 0.001	979 ± 28	0.7	633	323
K2-90 b	2.6 ± 0.1	9.0 ± 2.5	3.3 ± 0.9	0.096 ± 0.001	502 ± 8	1.4	157	139
K2-124 b	2.9 ± 0.1	10.5 ± 2.6	6.8 ± 1.8	0.049 ± 0.002	442 ± 17	1.2	151	383

Note. M_p comes from the probabilistic mass–radius relation of Wolfgang et al. (2016); we apply this relation uniformly for the sake of homogeneity, but note that M_p (and its uncertainty) may be underestimated for the larger planets in our sample. K_pred is the predicted RV semi-amplitude, a is semimajor axis, and ${T}_{\mathrm{eq}}$ is equilibrium temperature assuming a Bond albedo of 0.3; the uncertainties in these parameters propagate from the formal transit and stellar parameter estimates but are likely underestimated owing to uncertainties in the underlying models and assumptions. g is surface gravity, H is atmospheric scale height, and δ_TS is the predicted amplitude of atmospheric features accessible via transmission spectroscopy; these estimates have large uncertainties owing to a number of factors, such as uncertainties in the stellar parameters, unknown planet masses, and assumed atmospheric compositions.

Download table as:  ASCIITypeset image

In order to assess the potential for future characterization studies of these planets, we computed their physical properties using the parameter estimates from our transit analyses in Tables 3 and 4, as well as the stellar parameters in Table 6. We first computed updated planet radii and used the probabilistic mass–radius relation of Wolfgang et al. (2016) to predict the masses of these planets. We then computed the orbital semimajor axis (in physical units), predicted RV semi-amplitude (assuming circular orbits), equilibrium temperature (assuming a Bond albedo of 0.3), surface gravity, and atmospheric scale height (assuming a hydrogen-dominated atmosphere, i.e., a mean molecular weight of 2). Finally, we used the formalism of Miller-Ricci et al. (2009) to estimate the amplitude of features in the planets' transmission spectra. We list these derived planet properties in Table 7.

K2-53 and K2-174 are moderately bright in the optical (Kp ∼ 12), and each hosts an R_p ∼ 2.7 ${R}_{\oplus }$ planet with a predicted RV semi-amplitude of K_pred ∼ 3 m s⁻¹. These are thus potentially good targets for precision RV mass measurements with current optical spectrographs such as HIRES and HARPs, which would enable studies of the densities and bulk composition of temperate sub-Neptunes. K2-289 b and K2-87 b are both sub-Saturns orbiting moderately IR-bright stars (J < 12) with K_pred = 13.6 ± 5.2 and K_pred = 10.3 ± 3.4 m s⁻¹, respectively, making them good targets for one of the high-precision "red" (IR) spectrographs currently under development (e.g., Subaru/IRD, Tamura et al. 2012; Spirou, Artigau et al. 2014; HPF, Mahadevan et al. 2015). Mass measurements for these planets would add crucial data points to the relatively small population of similarly sized planets with well-measured densities, which exhibit a diverse range of core mass fractions (Petigura et al. 2017). K2-90 and K2-124 are similarly IR-bright but host temperate sub-Neptunes with K_pred = 3–6 m s⁻¹, making them also potentially interesting targets for red spectrographs.

Due to the small size and moderate infrared brightness of the host star, K2-124 b could be an interesting target for atmospheric studies via transmission spectroscopy with JWST. The predicted amplitude of features in the planet's transmission spectrum would be roughly 350 ppm (with large uncertainties owing to the unknown planet mass and atmospheric mean molecular weight). By combining multiple transit observations, these features should be detectable with JWST. This 2.5 ${R}_{\oplus }$ planet receives only ∼8 times the incident radiation as Earth, so such studies would probe the atmospheric properties of temperate sub-Neptunes.

The combination of K2 and Spitzer data in this work demonstrates the synergy between transit detection and follow-up with space-based instruments. Because TESS has a pixel scale five times larger than Kepler (21'' vs. 4''), the frequency of stellar blends will likely be significantly larger, which will increase the utility of transit follow-up. The upcoming CHEOPS mission has a similar pixel scale to Spitzer and will thus prove useful for TESS follow-up in much the same way that Spitzer has for K2.