GROUND-BASED TRANSIT OBSERVATION OF THE HABITABLE-ZONE SUPER-EARTH K2-3D

K2, the extended Kepler mission (Howell et al. 2014), is increasing the number of known small transiting planets orbiting nearby M dwarfs, including potentially habitable ones, providing promising targets for future atmospheric studies of small planets by, e.g., the James Webb Space Telescope (JWST). However, the survey duration of K2 per field (∼80 days) is much shorter than Kepler's prime mission (four years), decreasing the number of observable transits for a given planet by a factor of ∼18. The smaller time baseline per field inhibits precise measurement of the orbital period of transiting planets, as well as transit timing variations (TTVs). Transit followup observations by other telescopes are thus an important way to improve the ephemeris, which is essential in planning future observations, as well as to study the dynamics of multi-planet systems via TTVs.

K2-3d is one such planet worthy of additional followup. K2-3d was discovered around a nearby M0 dwarf in Campaign 1 of the K2 mission (Crossfield et al. 2015, hereafter C15), as the outermost planet among three small planets transiting the same host star. Independent detections of the planetary signal were also reported by Huang et al. (2015) and Vanderburg et al. (2016), and radial velocity measurements of this system were reported by Almenara et al. (2015) and Dai et al. (2016). K2-3d has an orbital period of 44.6 days, which corresponds to a semimajor axis of 0.208 AU, where the planet receives an incident flux of ${1.51}_{-0.43}^{+0.57}$ times that of the Earth (C15). This means that this planet is probably located at the inner edge of (according to the recent Venus model, Kasting et al. 1993) or well within (in the case of synchronous rotation, Yang et al. 2014; Kopparapu et al. 2016) the habitable zone around the M-dwarf host star. In addition, the composition of this planet could be dominated by rock, given its radius of ${1.52}_{-0.20}^{+0.21}$ R_⊕ (C15), which is on the boundary between rocky and volatile-rich planets for close-in planets (e.g., Weiss & Marcy 2014; Rogers 2015), although the boundary for cooler planets is less clear. Given the brightness of the host star (V = 12.17, J = 9.42), K2-3d is currently one of the best targets for spectroscopic characterizations of potentially habitable planets. However, the relatively long orbital period allowed K2 to observe only two transits within its observing campaign, leaving large uncertainties in the transit parameters, particularly in the orbital period.

Beichman et al. (2016, hereafter B16) observed two additional transits of K2-3d using Spitzer (IRAC 4.5 μm), aiming at improving the transit ephemeris and other parameters. Combining the transit timing data of K2 and Spitzer, they attempted to reduce the uncertainty of the transit timing prediction in the JWST era (around 2019 December) from ∼6 hr to ∼0.4 hr. However, because of the shallow transit depth of K2-3d (∼0.7 mmag), which is similar in amplitude to the systematic noise afflicting Spitzer light curves, a great deal of caution must be exercised when analyzing the data, and additional transit observations can yield significant further improvement of the ephemeris. Furthermore, there is the potential for dynamic interactions between the planets in the system to cause measurable TTVs, which could only be revealed by additional observations.

In this paper we report the first ground-based transit observation of K2-3d, which was achieved by a 2 m class telescope and a multi-band imager. Although observing the tiny transit of K2-3d is challenging for a ground-based 2 m class telescope, the multi-band simultaneous data help us not only to increase the statistical significance but also to reduce the systematics caused by the Earth's atmospheric extinction.

The rest of this paper is organized as follows. We describe the observation and reduction in Section 2. In Section 3 we reanalyze the K2 data in order to refine the transit parameters, with which we search for the transit signal in the MuSCAT data in Section 4. We then reanalyze the Spitzer data and refine the transit ephemeris in Section 5. Finally, our results are discussed and summarized in Section 6.

We conducted a photometric observation of K2-3 on 2016 March 2 UT, when a primary transit of K2-3d was predicted to occur, by using the Multi-color Simultaneous Camera for studying Atmospheres of Transiting exoplanets (MuSCAT, Narita et al. 2015a) mounted on the 188 cm telescope at Okayama Astrophysical Observatory (OAO) in Japan. MuSCAT is a three-channel imager for the Sloan ${g}_{2}^{\prime }$, ${r}_{2}^{\prime }$, and ${z}_{{\rm{s}},2}$ bands (hereafter g, r, and z bands, respectively), recently developed for validating (Hirano et al. 2016) and characterizing (Narita et al. 2015b; Fukui et al. 2016) transiting planets. MuSCAT provides a field of view of 61 × 61 for all channels, within which only two sufficiently bright comparison stars, namely TYC 4923-662-1 (V = 11.88, J = 10.67) and GSC2.3 S4RH000119 (V = 13.16, J = 12.07), were simultaneously imaged with the relatively bright target star K2-3. The exposure time was set to 60 s for all bands. The telescope was defocused such that the FWHM of the stellar point-spread function (PSF) was within the ranges of 15–21, 18–24, and 19–25 pixels for the g, r, and z bands, respectively. We started the observation at 12:29 UT and continued it through 19:06 UT, which resulted in the total time coverage of 6.6 hr. During the observation the sky was clear, with a waning moon (age of 23 days) for the last two hours.

All the observed images were dark-subtracted and flat-fielded in a standard manner. To create flat-field images, we median-combined 50 dome-flat images for each band that were obtained on the observing night. After the dark–flat correction, we applied a nonlinearlity correction for each CCD. Then aperture photometry was performed for the target and the two comparison stars by using the customized code (Fukui et al. 2011), which calculates the theoretical error for each flux by taking account of the photon noise, detector noise, and scintillation noise. The aperture radius was optimized for each band such that the dispersion of the differential light curve, which was created by taking the magnitude difference between the target star and the ensemble of the two comparison stars, was minimum. As a result, radii of 20, 22, and 21 pixels were selected for the g, r, and z bands, respectively. All the time stamps were then converted into barycentric JD (BJD) using the code by Eastman et al. (2010). The derived differential light curves are shown in Figure 1.

Zoom In Zoom Out Reset image size

In those light curves, we found a 20 minutes long, 2–4 mmag dip (depending on wavelengths) near the end of the observation (BJD ∼ 2457450.28). We attributed this feature to systematics of the Earth's atmospheric origin rather than an astrophysical phenomenon, because it coincided with a distortion of the raw light curves of both the target and comparison stars (see the middle panel of Figure 1 for the case of the comparison stars). This systematics was not correctable, even with the new approach to correct the second-order extinction effect described in Section 4.2, and we simply discarded this part (BJD > 2457450.255, shaded in Figure 1) from the rest of the analyses. The numbers of data points to be used for the analyses are 307, 307, and 308 for the g, r, and z bands, respectively.

Because the transit signal is not obvious in the MuSCAT light curves, it is difficult to identify it without any prior information about the shape of the transit light curve. We therefore utilize the transit parameters from K2 as the prior information. Since the original discovery of K2-3d by C15, several new K2 data pipelines have been developed that may achieve higher photometric precision than previous efforts. In our testing we found the light curve for K2-3 produced by the EVEREST pipeline (Luger et al. 2016) to be of a higher quality than others, as measured by the out-of-transit (OOT) scatter. We therefore update the transit parameters using this light curve. First, we reproduce the systematic-corrected light curve for K2-3 (EPIC 201367065) using the EVEREST pipeline by modeling the systematics with the transit parts masked to avoid underestimation of the transit depth, following the suggestion by Luger et al. (2016). Next, we extract individual transit light curves, each covering ±5 × T₁₄ from the transit center, where T₁₄ is the transit duration (∼4 hr). For each light curve we assign rms of the OOT residuals to each flux error bar. We note that the time-correlated-noise (red-noise) factor β, which is the ratio of the observed rms for a binned residual light curve to the rms expected from the dispersion of the unbinned residual light curve (Pont et al. 2006; Winn et al. 2008), becomes unity for both light curves, indicating no sign of red noise.

We then fit the two transit light curves simultaneously using a customized Markov chain Monte Carlo (MCMC) code (Narita et al. 2007, 2013). The code uses the transit model by Ohta et al. (2009) with the following parameters: the mid-transit time T_c for each transit, the orbital period P, the ratio of planetary radius to stellar radius R_p/R_s, the scaled semimajor axis a/R_s, the transit impact parameter b, and the quadratic limb-darkening coefficients u₁ and u₂. We fix P at the value in B16 (44.55765 days) and assume a circular orbit. For u₁ and u₂, we let them be free but impose Gaussian priors that we determine from the tabulated values of Claret et al. (2012), with location and scale parameters equal to the mean and twice the standard deviation of all entries that satisfy 3700 K ≤ T_eff ≤ 4100 K and log₁₀g ≥ 4.5, respectively. The adopted priors are listed in Table 1. Additionally, we impose spectroscopic constraints from C15 in the form of the Gaussian prior a/R_s ∼ 79.6 ± 10.4. Each light curve is also fit with a linear function to take into account the gradient of the baseline:

$$\begin{eqnarray}&&f(t)={{ \mathcal F }}_{\mathrm{tr}}(t)\times ({k}_{0}+{k}_{t}{\rm{\Delta }}t),\end{eqnarray} \tag{ 1 }$$

where f(t) is normalized flux at time t, ${\rm{\Delta }}t\equiv t-{T}_{c}^{\prime }$ is the differential time from a fixed initial guess for the mid-transit time ${T}_{c}^{\prime }$, ${{ \mathcal F }}_{\mathrm{tr}}(t)$ is the transit light-curve model on a flux scale, and k₀ and k_t are free parameters. Note that the model flux is computed by integrating time sampling models of 1 minute over the K2 integration time (30 minutes). We perform five independent MCMC runs with 10⁶ steps each, and create merged posterior probability distributions for the parameters. The resultant median and 1σ uncertainties of R_p/R_s, a/R_s, b, and T₁₄ are summarized in Table 2. We confirm that our refined values are consistent with those reported by C15 within ∼1σ.

Table 2.  Refined Transit Parameters for K2-3d

Parameter	This work	C15
a/Rs	${78.7}_{-5.0}^{+4.0}$	78.7 ${}_{-13.0}^{+6.7}$
b	${0.15}_{-0.14}^{+0.19}$	${0.45}_{-0.28}^{+0.23}$
Rp/Rs	${0.02484}_{-0.00054}^{+0.00056}$	${0.0248}_{-0.0010}^{+0.0014}$
T14 (hr)	${4.17}_{-0.07}^{+0.09}$	3.98 ${}_{-0.15}^{+0.17}$

Download table as:  ASCIITypeset image

4.1. 1D Search for ${T}_{{\rm{c}}}$

To extract the transit signal from the MuSCAT light curves, we conduct a one-dimensional (1D) search for T_c. In this search we continuously vary T_c with a step size of 0.002 days through the time window $[{t}_{1}-{T}_{14}/2$, ${t}_{2}+{T}_{14}/2$], where t₁ and t₂ are the start and end times of the observation, respectively. At each step we minimize the χ² value for a transit+systematic model for each light curve, fixing all the transit parameters except for T_c at either the values derived in Section 3 (a/R_s, b, R_p/R_s, and P) or the theoretical values given by Claret et al. (2012) (u₁ and u₂, listed in Table 1). For the systematic model, we first apply a simple linear function in the form of

$$\begin{eqnarray}&&{\rm{\Delta }}{m}_{{\rm{t}}}={{ \mathcal M }}_{\mathrm{tr}}(t)+{k}_{0}+{k}_{t}{\rm{\Delta }}t+{k}_{Z}Z,\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Delta }}{m}_{{\rm{t}}}$ is the differential magnitude between the target star and the ensemble of the comparison stars, ${{ \mathcal M }}_{\mathrm{tr}}$ is the transit light-curve model in magnitude, Z is the theoretical airmass calculated from zenith distance, and k₀, k_t, and k_Z are coefficients for normalization, time, and airmass, respectively. Note that we do not incorporate any instrumental variables for modeling systematics, such as the stellar displacements on the detectors (Δx and Δy) and FWHM of the stellar PSF, because we find no correlation between the light curves and these variables.

In the top panel of Figure 2, we show the χ² difference between the transit and null-transit models, Δχ², as a function of T_c. A negative value means that the transit model is superior to the null-transit one. We note that the absolute (reduced) χ² values for the null-transit model are 362.5 (1.19), 346.9 (1.14), and 344.9 (1.13) for the g, r, and z bands, respectively. We also calculate the red-noise factor β, for which we take the median value for the binning sizes ranging from 5 through 20 minutes, as a function of T_c. The results are shown in the bottom panel of the same figure.

Zoom In Zoom Out Reset image size

All three Δχ² plots show similar patterns, with roughly three local minima (valleys). We mark these valleys (for the r band as a representative) by dashed vertical lines, and label them as Valley1 through Valley3 from left to right. Of these, Valley1 gives the minimum Δχ² values of 14.6 and 14.3 for the r and z bands, respectively. The β values are also minimal around Valley1 with values close to unity for the r and z bands. Moreover, the location of Valley1 is marginally (<3σ) consistent with the expected T_c from B16 as indicated by an orange bar in Figure 2.

On the other hand, the minimum Δχ² and β values for the g band both appear around Valley3, which is inconsistent with the other bands. However, β for the g band is relatively large as a whole (>1.2), indicating that the result for this band is less reliable than the other two bands due to the likely higher red noise. Therefore, at this point Valley1 is the most likely transit signal, although its significance with respect to the other valleys is not high enough, given that the χ² difference between Valley1 and Valley2 (Valley3) is 1.4 (4.3) and 6.2 (9.7) for the r and z bands, respectively. In the top panels of Figure 3 we show the systematic-corrected light curves where we assume that the transit signal is present at Valley1.

Zoom In Zoom Out Reset image size

4.2. New Approaches to Correct Second-order Extinction

The relatively large β value for the g-band light curve indicates that it is significantly affected by systematic noise. We attribute the majority of the noise to second-order extinction by the Earth's atmosphere. In this section we introduce new approaches to reduce this effect, aiming at enhancing or refuting the evidence for the tentative transit signal. We begin with a brief description of the second-order extinction effect in Section 4.2.1, then introduce the new approaches based on single-band data (Section 4.2.2) and multi-band data (Section 4.2.3).

4.2.1. General Description of Second-order Extinction

The second-order extinction effect arises when the target and comparison stars have different spectral types. In such a case, because the extinction coefficient (or the attenuation by the Earth's atmosphere) is a function of wavelength (see Figure 4), the effective extinction coefficient, which is the weighted average of the extinction coefficient over the passband of a given filter for a given stellar spectrum, differs between the two stars. Consequently, the differential extinction between the two stars varies along with the change in atmospheric conditions that is dominated by airmass, causing systematic trends in the differential light curves (e.g., Young et al. 1991; Everett & Howell 2001; Bailer-Jones & Lamm 2003; Mann et al. 2011).

Zoom In Zoom Out Reset image size

This effect can be expressed in mathematical form as follows. The relation between the observed and intrinsic (unattenuated) brightness of the target and comparison stars for a given filter can be written as

$$\begin{eqnarray}&&{m}_{{\rm{t}}}={m}_{{\rm{t}},0}+{A}_{{\rm{t}}}Z+{I}_{{\rm{t}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{m}_{{\rm{c}}}={m}_{{\rm{c}},0}+{A}_{{\rm{c}}}Z+{I}_{{\rm{c}}},\end{eqnarray} \tag{ 4 }$$

where ${m}_{{\rm{t}}}$ (${m}_{{\rm{c}}}$) and ${m}_{{\rm{t}},0}$ (${m}_{{\rm{c}},0}$) are the observed and intrinsic magnitudes of the target (comparison) star integrated over the passband, respectively, ${A}_{{\rm{t}}}$ (${A}_{{\rm{c}}}$) is the effective extinction coefficient for the target (comparison) star at airmass = 1, Z is airmass, and ${I}_{{\rm{t}}}$ (${I}_{{\rm{c}}}$) is a constant value for the target (comparison) star that is accounted for by the instrumental throughput including the optics and quantum efficiency of the detector. Assuming that ${m}_{{\rm{c}},0}$ is constant over time, one can express the differential brightness of the target star as

$$\begin{eqnarray}{\rm{\Delta }}{m}_{{\rm{t}}} & \equiv & {m}_{{\rm{t}}}-{m}_{{\rm{c}}}\\ & = & {m}_{{\rm{t}},0}+({A}_{{\rm{t}}}-{A}_{{\rm{c}}})Z+C,\end{eqnarray} \tag{ 5 }$$

where C ($={I}_{{\rm{t}}}-{I}_{{\rm{c}}}-{m}_{{\rm{c}},0}$) is a constant. The second-order extinction effect arises when $({A}_{{\rm{t}}}-{A}_{{\rm{c}}})Z$ changes significantly during an observation compared with the measurement errors ${\rm{\Delta }}{m}_{{\rm{t}}}$.

The effect is expected to be large on our g-band observation for the following reasons. First, the differential extinction becomes larger at shorter wavelengths where the gradient of atmospheric transmittance is large due to scattering by atmospheric molecules and aerosols (see Figure 4). Thus, the effect on the g band is generally larger than on redder bands. Second, the color difference between the target and the ensemble comparison stars in our observation is moderately large, with ${\rm{\Delta }}(V-J)=1.33$. Such a situation is common for targets that are bright M dwarfs, which are only sparsely distributed on the sky. Last, our observation was conducted at a low-altitude (372 m) observatory, where the overhead atmosphere is thicker than at most of the great astronomical sites in the world.

Assuming that ${A}_{{\rm{t}}}$ and ${A}_{{\rm{c}}}$ are constant over the night, we estimate ${A}_{{\rm{t}}}$ and ${A}_{{\rm{c}}}$ for each band for our observation by using Equations (3) and (4). The results are summarized in Table 3. We estimate the differential extinction coefficient, ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$, to be −8, 4, and 2 mmag for the g, r, and z bands, respectively, indicating that the g band gives the largest difference between ${A}_{{\rm{t}}}$ and ${A}_{{\rm{c}}}$ as expected. Because the airmass value differs by more than 0.7 during the observation, the second-order extinction effect is expected to be considerable in all three bands, which in fact can be seen as long-term trends, prominently in the g band, in Figure 1. We note that although ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$ is expected to be negative for all bands given the color difference between the target and comparison stars, we can attribute the positive values for the r and z bands to the deep TiO absorptions expected to exist around the redder side in each band in the target star's spectrum.

Table 3.  Effective Extinction Coefficients for Our Observation

	g band		r band		z band
	Value	rms	Value	rms	Value	rms
	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)
${A}_{{\rm{t}}}$	0.445	0.024	0.290	0.016	0.166	0.010
${A}_{{\rm{c}}}$	0.453	0.025	0.286	0.016	0.164	0.010
${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$	−0.008	⋯	0.004	⋯	0.002	⋯

Download table as:  ASCIITypeset image

Such a long-term trend, however, can be corrected at the first order by estimating ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$ using the airmass function, under the assumption that ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$ is constant over the night. This approach is most commonly used and was also applied in Section 4.1.

However, in reality, ${A}_{{\rm{t}}}$ and ${A}_{{\rm{c}}}$ can temporarily fluctuate due to the temporal changes in atmospheric conditions such as barometric pressure, precipitable water, and the amount and nature of aerosols (e.g., Hill et al. 1994; Stubbs et al. 2007; Burke et al. 2010). In fact, the raw light curves deviate from a smooth function of airmass (Equations (3) and (4)) by 25, 16, and 10 mmag in rms for the g, r, and z bands, respectively (see Table 3), indicating that the extinction coefficients varied considerably with time at the level of several per cent. The fluctuations of ${A}_{{\rm{t}}}$ and ${A}_{{\rm{c}}}$ could also cause time variations of ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$, which could leave residual systematics that cannot be corrected by the conventional approach using the airmass function. For example, fluctuations of 5% in ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$ could cause systematic signals in the g-band light curve with amplitudes of 0.5–0.7 mmag in the airmass range Z = 1.2–1.8, which are comparable to the transit depth of K2-3d (∼0.7 mmag).

4.2.2. Single-band Approach

Here we introduce new approaches to correct the systematic noise originating from the second-order extinction. First we assume that fluctuations in the atmospheric transparency happen achromatically. In this case the ratio of the effective extinction coefficients for the target and comparison stars is constant over time, such that

$$\begin{eqnarray}&&{A}_{{\rm{t}}}={k}_{c}{A}_{{\rm{c}}},\end{eqnarray} \tag{ 6 }$$

where k_c is a constant factor. In this case, Equation (3) can be transformed into

$$\begin{eqnarray}&&{m}_{{\rm{t}}}={m}_{{\rm{t}},0}+{k}_{c}{m}_{{\rm{c}}}+C^{\prime} ,\end{eqnarray} \tag{ 7 }$$

where $C^{\prime}$ ($={I}_{{\rm{t}}}-{k}_{c}{I}_{{\rm{c}}}-{k}_{c}{m}_{{\rm{c}},0}$) is a constant. For a transit light curve, this can be rewritten as

$$\begin{eqnarray}&&{m}_{{\rm{t}}}(t)={{ \mathcal M }}_{\mathrm{tr}}(t)+{k}_{0}+{k}_{t}{\rm{\Delta }}t+{k}_{c}{m}_{{\rm{c}}}(t),\end{eqnarray} \tag{ 8 }$$

in which we add a linear function of ${k}_{0}+{k}_{t}{\rm{\Delta }}t$ for the baseline gradient. This equation differs from Equation (2) in two ways: (1) the raw magnitude of the target star is directly modeled as a function of the raw magnitude of the comparison stars, and (2) it does not rely on the theoretical airmass Z. Using Equation (8), we rerun the 1D search for T_c, allowing k_c to be free. For a computational reason, the uncertainty of ${m}_{{\rm{c}}}(t)$ in Equation (8) is propagated beforehand to that of ${m}_{{\rm{t}}}(t)$ so that we can treat the uncertainty of ${m}_{{\rm{c}}}(t)$ as zero during the fitting process. The results are plotted with orange lines in Figure 5. The effects are visible in the g band; both the depths of the suspicious valleys in the Δχ² plot and the overall β values are slightly reduced. Nevertheless, the strongest peak is still present around Valley3. On the other hand, for the r and z bands, we obtain almost identical results to those derived in Section 4.1. In the panels in the middle row of Figure 3, we show the differential light curves that are corrected using the new formula and the transit model that assumes that the transit signal is at Valley1.

Zoom In Zoom Out Reset image size

4.2.3. Multi-band Approach

Next, we consider a more realistic situation in which the extinction coefficients of ${A}_{{\rm{t}}}$ and ${A}_{{\rm{c}}}$ are variable in time. In general, the variations of the extinction coefficients depend on wavelength. In other words, the extinction coefficient for the target star for a given band at a given time can in principle be predicted from those for the comparison stars for multiple bands at the same moment. We assume here that A_t for a given band λ at time t can be expressed as a linear combination of A_c for multiple bands at the same time as follows:

$$\begin{eqnarray}&&{A}_{{\rm{t}},\lambda }(t)=\displaystyle \sum _{\lambda ^{\prime} =\{g,r,z\}}{k}_{c,\lambda ^{\prime} }{A}_{{\rm{c}},\lambda ^{\prime} }(t).\end{eqnarray} \tag{ 9 }$$

Substituting Equation (9) into Equation (3), one can rewrite Equation (8) as

$$\begin{eqnarray}{m}_{{\rm{t}},\lambda }(t) & = & {{ \mathcal M }}_{\mathrm{tr}}(t)+{k}_{0}+{k}_{t}{\rm{\Delta }}t\\ & & +\displaystyle \sum _{\lambda ^{\prime} =\{g,r,z\}}{k}_{c,\lambda ^{\prime} }{m}_{{\rm{c}},\lambda ^{\prime} }(t).\end{eqnarray} \tag{ 10 }$$

Note that this equation is equivalent to Equation (8) when ${k}_{c,\lambda ^{\prime} (\lambda ^{\prime} \ne \lambda )}=0$.

The point of this new formula is to correct the extinction for the target star by using the fluxes of the comparison stars in multiple bands. This new approach, however, has to be applied with caution. If the major sources of attenuation between two arbitrary bands are different, then the correlation between the extinction variations for these two bands could be weak. In such a case, incorporating fluxes of different bands in extinction correction could introduce additional, unwanted systematics. In fact, while scattering by aerosols is the dominant source at shorter wavelengths, water absorption become dominant in redder bands (see Figure 4), and the timescale of variability of these two sources could be different.

Keeping this risk in mind, we apply this new formula to the MuSCAT data. First, to determine which coefficients to include in ${k}_{c,\lambda }$, we fit each light curve with different combinations of ${k}_{c,\lambda }$, assuming that the transit signal is present at Valley1. We note that we always include ${k}_{c,\lambda ^{\prime} }$ when we fit the λ'-band light curve, because the extinction for the target star for a given band is primarily correlated with that for the comparison stars for the same band. In Table 4, we summarize the statistics of the fittings with different combinations of ${k}_{c,\lambda }$, in addition to those for the conventional approach (denoted as k_Z) for reference. The statistics include Bayesian information criteria BIC $\equiv \,{\chi }^{2}+k\mathrm{ln}N$, where k is the number of free parameters and N is the number of data points (Schwarz 1978), rms, β, and rms × β. Although BIC is often used to choose a baseline light-curve model from among several candidates, it does not adequately reflect red noise, which is now of particular concern. We therefore select the best baseline model such that (1) rms × β is minimum, and (2) BIC is minimum if there is another competitive model that has a similar value of rms × β to the minimum one. The selected combinations are {${k}_{c,g}$, ${k}_{c,r}$}, {${k}_{c,r}$, ${k}_{c,z}$}, and {${k}_{c,z}$} for the g, r, and z bands, respectively. We note that even when we conduct the same test assuming that no transit signal is present in the data, we obtain the same combinations of ${k}_{c,\lambda }$. The statistics for the null-transit model are appended in Table 4. This is also true when we assume that a transit signal is present at different valleys. Therefore, the choice of the baseline models does not produce any bias in the transit detection.

Table 4.  Comparison of Statistics between the Baseline Models^a

	Valley1 (Tc = 2457450.016)				Null Transit
Model	BIC	rms	β	rms × β	BIC	rms	β	rms × β
		(mmag)		(mmag)		(mmag)		(mmag)
g band
${k}_{{\rm{Z}}}$	377.7	1.31	1.50	1.97	379.6	1.30	1.52	1.98
${k}_{c,g}$	347.7	1.25	1.24	1.55	353.5	1.26	1.32	1.66
${{\boldsymbol{k}}}_{{\boldsymbol{c}},{\boldsymbol{g}}}+{{\boldsymbol{k}}}_{{\boldsymbol{c}},{\boldsymbol{r}}}$	319.1	1.19	1.00	1.19	320.9	1.20	1.00	1.20
${k}_{c,g}+{k}_{c,z}$	323.4	1.20	1.00	1.20	327.1	1.21	1.00	1.21
${k}_{c,g}+{k}_{c,r}+{k}_{c,z}$	323.5	1.19	1.00	1.19	325.7	1.20	1.00	1.20
r band
kZ	349.4	0.87	1.00	0.87	364.1	0.89	1.26	1.12
${k}_{c,r}$	348.5	0.87	1.00	0.87	364.5	0.89	1.24	1.10
${k}_{c,g}+{k}_{c,r}$	345.1	0.85	1.21	1.03	363.1	0.88	1.46	1.28
${{\boldsymbol{k}}}_{{\boldsymbol{c}},{\boldsymbol{r}}}+{{\boldsymbol{k}}}_{{\boldsymbol{c}},{\boldsymbol{z}}}$	348.2	0.86	1.00	0.86	363.7	0.88	1.22	1.08
${k}_{c,g}+{k}_{c,r}+{k}_{c,z}$	312.8	0.81	1.18	0.96	332.9	0.84	1.56	1.31
z band
kZ	352.6	1.21	1.12	1.36	362.1	1.23	1.28	1.58
${{\boldsymbol{k}}}_{{\boldsymbol{c}},{\boldsymbol{z}}}$	353.0	1.21	1.11	1.35	362.8	1.23	1.24	1.53
${k}_{c,g}+{k}_{c,z}$	346.1	1.19	1.35	1.60	355.0	1.21	1.46	1.76
${k}_{c,r}+{k}_{c,z}$	335.8	1.17	1.44	1.68	343.7	1.18	1.55	1.83
${k}_{c,g}+{k}_{c,r}+{k}_{c,z}$	316.1	1.12	1.35	1.51	320.8	1.13	1.42	1.60

Note.

^aThe final models are indicated by boldface.

Download table as:  ASCIITypeset image

The fact that the multi-coefficient models are selected for the g and r bands indicates that the multi-band approach is superior to the single-band approach for these two bands. On the other hand, for the z band, the single-coefficient model is selected, meaning that the multi-band approach does not improve but rather degrades the fit. In fact, while the BIC and rms values are reduced by the multi-coefficient models, the β value increases substantially so that rms × β increases, implying that in the case of the z band the information from different bands introduces additional systematic noise that outweighs its benefit.

We then rerun the 1D search for T_c with the multi-band approach. Just for comparison, we also try with the coefficients of {${k}_{c,g}$, ${k}_{c,r}$, ${k}_{c,z}$} for the z band, which is the best combination among the multi-coefficient models. The results are plotted with purple lines in Figure 5. The effect of the multi-band approach is obvious in the g band, where the suspicious peak around Valley3 is now well suppressed and β becomes unity for most T_c values. At the same time, however, no strong peak remains even at Valley1, which is probably a consequence of the fact that the true transit signal is overwhelmed by the photometric noise that is now dominated by statistical noise. Nevertheless, we can still see a possible transit signal in the g-band data (see Section 4.4).

In the case of the r band, the multi-band approach slightly enhances the significance of Valley1 over the other valleys. In the case of the z band, on the other hand, the multi-band approach does not enhance Δχ² at Valley1 but rather enhances the other valleys, probably due to the additional systematics. In the bottom panels of Figure 3, we show the systematic-corrected light curves from the multi-band approach for the transit model at Valley1.

4.3. p-value Calculation

4.4. Consistency Check for T_c and R_p/R_s

To see the consistency of the tentative transit signal at Valley1 within the three light curves and with what we expect from the past observations, we perform an MCMC analysis for individual light curves. Before the analysis we rescale the magnitude errors in each light curve such that the reduced ${\chi }^{2}$ value for the best-fit transit+systematic model at Valley1 becomes unity and we further rescale them by a factor of β (the value for Valley1 listed in Table 4). During the MCMC analysis we fix a/R_s and b to the values derived from the K2 data in Section 3, fix u₁ and u₂ to the theoretical values, and allow T_c, R_p/R_s, k₀, k_t, and the selected coefficients of ${k}_{c,\lambda }$ to be free. For the g band, to avoid local minima we loosely constrain T_c to the range between 2457450.000 and 2457450.032. We perform five independent MCMC runs with 10⁶ steps each for each light curve, and calculate the merged posterior probability distributions of the parameters. The resultant median and 1σ uncertainties for T_c and R_p/R_s as well as ${k}_{c,g}$, ${k}_{c,r}$, and ${k}_{c,z}$ are listed in Table 5. Figure 6 shows two-dimensional correlation plots between R_p/R_s and T_c and between R_p/R_s and k_t. We find that R_p/R_s and T_c for all bands are largely consistent with each other within 2σ. Remarkably, we detect a marginal signal in the g-band data at ∼1σ level, which has R_p/R_s and T_c values that are consistent with those in the other bands, further supporting the scenario in which the transit signal is present at Valley1. We also find that the R_p/R_s values for the g, r, and z bands are consistent with the value from K2 (indicated by blue horizontal lines in Figure 6) within 1σ, 2σ, and 3σ, respectively. The slightly large discrepancy between the z band and K2 might be caused by systematics in the z-band light curve; since only the partial transit was covered by MuSCAT, the transit parameters such as R_p/R_s can easily be affected by systematics. In fact, R_p/R_s is correlated with k_t (the coefficient of the linear function for the baseline) as shown in the right panels in Figure 6, indicating that R_p/R_s can be affected by systematics even in the out-of-transit part of the light curves.

Zoom In Zoom Out Reset image size

Table 5.  MCMC Results for the Individual Light Curves

Parameter	g band	r band	z band
Tc	${0.0153}_{-0.0097}^{+0.0105}$	${0.0157}_{-0.0024}^{+0.0029}$	${0.0096}_{-0.0038}^{+0.0020}$
(BJD – 2450000)
Rp/Rs	${0.0165}_{-0.0104}^{+0.0092}$	${0.0323}_{-0.0048}^{+0.0042}$	${0.0417}_{-0.0054}^{+0.0049}$
${k}_{c,g}$	${0.785}_{-0.032}^{+0.030}$	0 (fixed)	0 (fixed)
${k}_{c,r}$	${0.302}_{-0.048}^{+0.050}$	0.961 ± 0.015	0 (fixed)
${k}_{c,z}$	0 (fixed)	0.060 ± 0.026	0.9925 ± 0.0047

Download table as:  ASCIITypeset image

Next, to derive a final T_c value and its uncertainty, we fit all the three light curves simultaneously with MCMC. In order to properly reflect the uncertainties of the a priori information about the transit shape, we also simultaneously fit the K2 light curves in addition to the MuSCAT light curves. In this analysis we allow the following parameters to be free: a/R_s, b, and R_p/R_s as common parameters for all data, T_c for each transit, u₁ and u₂ for K2 with the same priors listed in Table 1, and the coefficients for the baseline function for each light curve. For the baseline function we use Equation (1) for the K2 data, and Equation (10) for the MuSCAT data with the selected coefficients in Section 4.2.3. After performing five independent MCMC runs with 5 × 10⁶ steps each, we merge all the MCMC steps and calculate the median and 1σ uncertainties, yielding

$$\begin{eqnarray}&&{T}_{c}({\mathrm{BJD}}_{\mathrm{TDB}}-2450000)\\ &&\quad =7450.0137\pm 0.0031(\mathrm{stat}.)\pm 0.0030(\mathrm{sys}.),\end{eqnarray} \tag{ 11 }$$

where we conservatively add 0.0030 in the uncertainty as a systematic error that corresponds to half of the offset between the best-fit T_c for the individual fit to the r- and z-band light curves. The estimated T_c is about 25 minutes earlier than predicted from the B16 ephemeris (T_c = 7450.0319 ± 0.0071), although the significance of the discrepancy is only 2.2σ. We will revisit this issue in Section 5.

4.5. Short Summary of the Signal Evidence

Motivated by the discrepancy of 25 minutes between T_c observed by MuSCAT and that expected from B16, we revisit the Spitzer photometry to investigate the possibility that the orbital period estimated by B16 could have been biased by uncorrected residual systematics in the Spitzer data. Systematic signals in the Spitzer photometry are produced by the motion of the PSF on the detector coupled with variations in intrapixel gain. Because the amplitude of the transit signal from K2-3d is small relative to the systematic signal in the Spitzer data, imperfect systematic correction can easily cause biases in the transit parameters. One clue that this might be the case can be found by comparing the transit duration (T₁₄) derived from the K2 light curve with that of Spitzer; based on the transit parameters reported in B16, T₁₄ for Spitzer is calculated to be 4.67 and 4.44 hr for the transit epochs of E = 6 and 10, respectively, which are slightly longer than the value for K2 (${4.17}_{-0.07}^{+0.09}$ hr). This might be the result of the conservative approach taken by B16, in which they intended to verify the original discovery by C15 by independently modeling the Spitzer data using uninformative priors in their MCMC analysis. Furthermore, they modeled the Spitzer systematics separately from the transit model, so the error bars they report may not accurately reflect the uncertainties in the parameters of the systematics model.

Because of this, we can hope to improve upon the quality of the extraction of the transit signal from the Spitzer data by simultaneously fitting the K2 and Spitzer data with transit+systematic models. By simultaneously modeling the K2 data (or alternatively using informative priors), the efficacy of modeling Spitzer's systematics may be enhanced due to the additional leverage provided by the constraints on the transit parameters from the K2 data.

We extract the Spitzer photometry and model its systematics using the pixel-level decorrelation (PLD) method, following Deming et al. (2015) and B16. PLD uses a linear combination of the individual pixel basis vectors (PBVs) to model the effect of PSF motion on the detector, thus it does not require the calculation of centroids. We use a similar parameterization to B16, with a multiplicative (instead of additive) linear baseline function:

$$\begin{eqnarray}&&{\rm{\Delta }}{S}^{t}=\left(\displaystyle \sum _{i=1}^{N}{c}_{i}{\hat{P}}_{i}^{t}+{F}_{\mathrm{tr}}(t)\right)\times ({k}_{0}+{k}_{t}{\rm{\Delta }}t),\end{eqnarray} \tag{ 12 }$$

where the c_i are the PBV coefficients and ${\hat{P}}_{i}^{t}$ is the ith pixel value of the normalized pixel grid at time t:

$$\begin{eqnarray}&&{\hat{P}}_{i}^{t}=\displaystyle \frac{{P}_{i}^{t}}{{\sum }_{i=1}^{N}{P}_{i}^{t}}\end{eqnarray} \tag{ 13 }$$

Instead of using residual rms as a merit function, we select the best circular aperture for the photometry in our analysis by minimizing residual correlated noise, as measured by β. We find that both the residual rms and residual correlated noise are approximately minimized (β ∼ 1) for a radius of 3 Spitzer/IRAC pixels in the case of the E = 6 data. However, for E = 10, which is more obviously affected by systematics (as seen in the raw light curves), we find that while a similar aperture (r = 2.9 pixels) minimizes residual rms, a smaller aperture of radius r = 2.2 pixels significantly reduces residual correlated noise (β drops from ∼1.6 to ∼1.3). We note that β for the E = 10 data can be further decreased by increasing the complexity of the systematics model, either by using a larger pixel grid (e.g., 5 × 5 versus 3 × 3) or a higher-order PLD model (e.g., including additional free parameters to allow for second-order terms in the Taylor expansion that PLD is based on). However, our testing of these alternative systematics models showed that the increase in model complexity was not warranted, as determined by BIC. We attempt to compensate for any biases resulting from the time-correlated noise by inflating the photometric uncertainties by their β factor. Note that rigorous comparisons between the various methods used to correct Spitzer systematics have been made (Ingalls et al. 2016) in which PLD was among the top performers, displaying both high precision and repeatability. However, for any individual data set, one generally high-performing method may be less well suited than another, e.g., if the amplitude of the systematics is unusually high during the observations due to above average thermal instability. A more comprehensive assessment of the Spitzer K2-3d light curves might entail implementing some of these other methods to check for agreement, but we leave this for a future work.

For the transit model and MCMC we use the open-source Python packages PyTransit (Parviainen 2015) and emcee (Foreman-Mackey et al. 2013). The free parameters of the transit model are R_p/R_s, a/R_s, T_c for each transit, and the orbital inclination i. We assume a circular orbit and a quadratic limb-darkening law. We impose Gaussian priors on the limb-darkening coefficients, which are determined in the same manner as in Section 3 and are listed in Table 1. We also impose the Gaussian prior a/R_s ∼ 79.5 ± 10.4 in the same way as in Section 3. Note that we ignore any wavelength dependence of R_p/R_s caused by a possible atmosphere, which is beyond the sensitivity of the data we model in this work.

In Figure 7 we present a corner plot showing the posterior probability distributions (PPDs) of T₁₄, ${T}_{c,6}$, and ${T}_{c,10}$, where ${T}_{c,6}$ and ${T}_{c,10}$ are T_c for E = 6 and 10, respectively, and the correlations between them. T₁₄ converged to ${4.166}_{-0.070}^{+0.096}$ hr, which is fully consistent with that from the K2-only fit but is significantly shorter than the values from B16. This indicates that the transit fits by B16 were most likely biased by systematics that we are now able to correct by simultaneously fitting the K2 data and systematics model.

Zoom In Zoom Out Reset image size

${T}_{c,6}$ also converges to a unimodal distribution with a median value and 1σ credible interval of ${T}_{c,6}={7093.5634}_{-0.0020}^{+0.0024}$, which is earlier than the value of B16 by ∼2σ. This discrepancy could arise from the same cause as the shrinkage of T₁₄, as indicated by the slight positive correlation between T₁₄ and ${T}_{c,6}$ (see Figure 7). On the other hand, the PPD for ${T}_{c,10}$ is significantly bimodal (${T}_{c}-2457271\,\sim$ 0.790 and 0.807), with more probability mass in the earlier mode. This can be interpreted as the result of the higher levels of "uncorrectable" systematics in the data set, as evident from the β factor. We find that the ${T}_{c,10}$ PPD can be well fit by a triple Gaussian model with two higher-amplitude modes and a low-amplitude mode in between, as shown in Figure 8. This analysis yields a mean and standard deviation of (μ₁, σ₁) = (0.7899, 0.0017) and (μ₂, σ₂) = (0.8073, 0.0026) for the dominant modes. In Table 6 we summarize the T_c values for the available five transits estimated in this work, along with the values reported in the literature.

Zoom In Zoom Out Reset image size

Table 6.  List of Observed T_c for K2-3d

Epoch	Instrument	Tc(BJD) – 2450000
		This Work	C15	B16
0	K2	6826.2274 ± 0.0020	${6826.2232}_{-0.0043}^{+0.0037}$	⋯
1	K2	${6870.7860}_{-0.0019}^{+0.0024}$	${6870.7863}_{-0.0070}^{+0.0073}$	⋯
6	Spitzer	${7093.5633}_{-0.0020}^{+0.0024}$	⋯	${7093.5752}_{-0.0050}^{+0.0036}$
10	Spitzer	7271.7899 ± 0.0017 (best solution)	⋯	${7271.8007}_{-0.0017}^{+0.0019}$
		7271.8073 ± 0.0026 (2nd solution)
14	MuSCAT	7450.0137 ± 0.0043a	⋯	⋯

Note.

^aBoth the statistical and systematic errors are added quadratically.

Download table as:  ASCIITypeset image

When we choose the mean value of the highest amplitude (earlier) mode of the ${T}_{c,10}$ PPD, all five T_c from this work are well aligned, giving a χ² value for the best-fit linear function of 1.41 over three degrees of freedom. We obtain the refined linear transit ephemeris of

$$\begin{eqnarray}&&{T}_{c}=6826.2282(14)+44.55612(21)\times E,\end{eqnarray} \tag{ 14 }$$

where the number in parentheses represents the last two digits of 1σ uncertainty. The residuals of T_c from the best-fit linear function along with the ephemerides of previous work are shown in Figure 9.

Zoom In Zoom Out Reset image size

On the other hand, when we choose the mean value of the second-highest amplitude (later) mode of the ${T}_{c,10}$ PPD, χ² becomes 30.1, with a maximum dispersion of ∼15 minutes from the linear ephemeris. Although these dispersions might be caused by the gravitational forces from other planets, it is unlikely, considering the orbital architecture of this system (at least to our current knowledge), for K2-3d to show TTVs with an amplitude of more than 5 minutes (B16). We therefore argue that the earlier ${T}_{c,10}$ with a linear ephemeris is the preferred scenario under Occam's razor, i.e., the simplest explanation should be selected among competing hypotheses. Furthermore, analysis of the E = 10 light curve resulting from an aperture of larger radius (r = 2.9 pixels), which minimizes residual rms instead of correlated noise, produced a bimodal PPD for ${T}_{c,10}$ with higher probability mass in the later mode; we interpret this as evidence that the later mode is positively correlated with the strength of uncorrected residual systematics, and is therefore less likely to be the result of the transit signal.

We have conducted the first ground-based transit observation of K2-3d, a planet of 1.5 Earth radii that is considered to be inside the habitable zone, using the OAO 188 cm telescope and the multi-band imager MuSCAT. Although the depth of the transit is only 0.7 mmag, which is smaller than the photometric precision per exposure and vulnerable to systematics, we have marginally but consistently detected the transit signal thanks to a priori knowledge of the transit shape from the K2 light curves and the multi-band observations with the new extinction-correction approaches.

We have also revisited the Spitzer data obtained by B16. Unlike B16, we have simultaneously fit both the Spitzer and K2 data with both the transit model and Spitzer's systematic models, which has allowed us to reduce the residual systematics that biased the results of B16. Although the posterior distribution of ${T}_{c,10}$ becomes bimodal, the earlier solution has a higher likelihood and can explain all the transit timing data including the MuSCAT one well by a linear ephemeris. Adopting the earlier solution for ${T}_{c,10}$, we have refined the orbital period to be 44.55612 ± 0.00021 days, which is 130 s shorter than the value in B16. This discrepancy would accumulate to ∼80 minutes in the predicted transit times in 2019, and is thus critical for planning future followup observations by, e.g., JWST. This work highlights the importance of a careful analysis for a tiny transit signal that is comparable to the systematic signals, as well as the importance of independent followup observations to verify the previous results.

The transit depth of K2-3d is the second shallowest among those observed by ground-based telescopes, following 55 Cnc e (0.4 mmag transit observed by the 2.5 m Nordic Optical Telescope, de Mooij et al. 2014). Our observation demonstrates that ground-based photometric observations can play an important role in improving the transit ephemeris of small-sized, long-period planets, including potentially habitable ones, discovered by on-going and future space-based transit surveys such as K2 and the Transiting Exoplanet Survey Satellite (Ricker et al. 2015), whose survey durations are limited.

We have also demonstrated that a multi-band imager like MuSCAT is useful not only for obtaining multi-band light curves simultaneously but also for reducing the systematic noise that originates from the second-order extinction effect, especially for shorter wavelengths and for observations conducted at low-altitude observatories. This means that our new technique will also be useful in studying planetary atmospheres via transmission spectrophotometry, because the shorter wavelengths in the optical region are important for seeing the effect of Rayleigh scattering (e.g., Pont et al. 2008). Note, however, that there can be multiple timescale and wavelength dependences in the variations of the differential extinction coefficient ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$, because there are various sources of variations in extinction, including the temporal changes of pressure, winds, the nature of aerosols, and precipitable water. In fact, an uncorrectable dip was seen near the end of the MuSCAT light curves (see Figure 1), which might be caused by variations of ${A}_{{\rm{t}}}-{A}_{{\rm{c}}}$ on a different timescale. These complex systematic features cannot be corrected by the simple form of Equation (10), and require higher-order models or non-parametric approaches such as Gaussian processes (e.g., Rasmussen & Williams 2006; Gibson et al. 2012), which is left as future work.

This work was supported by the Astrobiology Center Project of National Institutes of Natural Sciences (NINS) (Grant Numbers AB281012 and JY280092). This work was also supported by JSPS KAKENHI (Grant Numbers JP25247026 and JP16K17660).