An Earth-sized Planet around an M5 Dwarf Star at 22 pc

K2-415 was observed by K2 in Campaign 16 from UT 2017 December 13 to UT 2018 February 25 in the long-cadence (30 minutes) mode. We reduced the target pixel files downloaded from the Mikulski Archive for Space Telescopes (MAST) website. ³³ Our pipeline for reducing K2 data, searching for transiting planets, and vetting the planet candidates was detailed in Hirano et al. (2018). The major challenge in reducing K2 data was to mitigate the systematic variation due to the rolling motion of the telescope along the boresight (Howell et al. 2014). We decorrelated the flux variation with the flux centroid motion using a method similar to that described by Vanderburg & Johnson (2014). We then detrended the light curve with a cubic spline as a function of time with a width of 0.75 day to remove long-term systematics and stellar activity. We note that this spline-detrending was only used to empirically normalize the light curve for planet detection. A more physically motivated modeling of the light curve using Gaussian Process regression was also conducted as described in Section 4.4. A (BLS; Kovács et al. 2002) search returned a strong detection with a signal detection efficiency (SDE) of 12.9 at an orbital period of 4.02 day. We checked this transit signal for odd–even variation and deep secondary eclipse, which are typical signs of a false positive due to eclipsing binaries. We only detected a 2.4 σ odd–even variation and a 0.9σ secondary eclipse (with the derived depth of 0.00014 ± 0.00015). K2-415 hence passed our initial vetting and was promoted for further follow-up observations.

K2-415 was also observed by TESS at a 2 minutes cadence with the TIC ID of 323687123 in Sectors 44, 45, and 46, between UT 2021 October 21 and December 30. In order to see if the transit signals by the same planet candidate (P = 4.02 days) are identifiable in the TESS data, we downloaded the PDCSAP FLUX light-curve data generated by the SPOC pipeline (Smith et al. 2012; Stumpe et al. 2012, 2014). After applying a 3.5 σ clipping to the flux data as well as correcting for the flux offsets between different sectors (panel (b) of Figure 1), we implemented the BLS analysis (Kovács et al. 2002) as in the case of the K2 light curve. We identified the same transiting-planet candidate (P = 4.02 days) as the highest power in the BLS periodogram, with an SDE of 10.0 (Figure 2). This P = 4.02 days planet candidate was independently detected by the TESS project, by which K2-415 was named "TOI-5557."

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to search for an additional transiting-planet candidate in the system, we performed a joint BLS analysis using both K2 and TESS light curves. However, no significant peak other than the P = 4.02 day candidate was identified in the BLS periodogram.

On UT 2018 June 18, we performed high-resolution imaging for K2-415 using the InfraRed Camera and Spectrograph (IRCS; Kobayashi et al. 2000) and the adaptive-optics system AO188 (Hayano et al. 2008), both mounted on the Subaru 8.2 m telescope. Adopting the fine-sampling mode (1 pix = 20 mas) and the ${K}^{{\prime} }-$band filter (≈2.1 μm), we imaged K2-415 with a five-point dithering. The exposure time for each dithering position was set to 3 s × 3 coadd = 9 s so that the peak count of K2-415's image stays within the linearity count regime. We obtained two sequences of the dithering pattern, giving a total on-sky integration time of 90 s.

Raw IRCS frames were reduced as in Hirano et al. (2016), and we aligned and median-combined the reduced images in order to suppress the background noise and thus achieve a higher flux contrast. The background flux scatter in the combined image was found to be 5.9 ADU, while the flux peak of K2-415 was ≈11,400 ADU. The full width at the half maximum (FWHM) of K2-415's image was measured to be 0133, and no secondary source was detected in the field of view (FoV) of 20'' × 20''. To gain a constraint on the magnitude of any secondary source, we computed a 5 σ contrast as a function of angular separation from the target, following the procedure in Hirano et al. (2016). The 5σ contrast curve is plotted in Figure 3, in which the inset displays the 4'' × 4'' image of K2-415. IRCS AO imaging achieved ${\rm{\Delta }}{m}_{{K}^{{\prime} }}$ of 6–7 mag at the angular separation of 1''.

Zoom In Zoom Out Reset image size

On the night of UT 2019 January 20, K2-415 was observed with the NESSI speckle imager (Scott 2019), mounted on the 3.5 m WIYN telescope at Kitt Peak. NESSI simultaneously acquires data in two bands centered at 562 and 832 nm using high-speed electron-multiplying CCDs (EMCCDs). We collected and reduced the data following the procedures described in Howell et al. (2011). The resulting reconstructed 832 nm band image achieved respective contrasts of Δmag ∼ 4 and Δmag ∼ 6 at separations of 02 and 1'' (see Figure 4).

Zoom In Zoom Out Reset image size

As part of the follow-up campaign for transiting-planet candidates identified by K2 and TESS, we conducted high-resolution spectroscopy using the InfraRed Doppler spectrograph (IRD; Tamura et al. 2012; Kotani et al. 2018) on the Subaru telescope. IRD covers the near-IR wavelengths from 970 to 1730 nm with a spectral resolution of R ≈ 70, 000. The IRD observations were carried out between UT 2019 January 15 and 2022 May 25, and a total of 42 IRD frames were obtained, including a few frames taken in the Subaru Strategic Program (SSP), which is a blind Doppler survey to find exoplanets around nearby mid- to late-M dwarfs (see, e.g., Harakawa et al. 2022). For the nine frames secured in 2022 January, only H-band spectra (1450 nm < λ) were obtained, due to technical issues in the YJ-band detector. We set each integration time to 300–1800 s with a median exposure of 1200 s. For each integration, we simultaneously took the reference spectrum of the laser-frequency comb (LFC), which is used for the estimation of the instantaneous instrumental profile (IP) of the spectrograph.

Raw IRD frames were reduced as in Hirano et al. (2020), and we extracted the wavelength-calibrated one-dimensional (1D) spectra for both stellar and reference (LFC) fibers. Based on those reduced spectra, we measured precise RVs using IRD's RV-analysis pipeline (Hirano et al. 2020); briefly speaking, the instantaneous IP for each spectral segment was first estimated from the LFC spectrum, with which individual stellar spectra were deconvolved. The stellar template spectrum for K2-415 was then generated by combining those IP-deconvolved spectra, in which telluric lines were also removed via theoretical telluric model fits or using a rapid rotator's spectrum. The combined template was found to have good signal-to-noise (S/N) ratios (typically >100 per pixel) for the whole spectral range covered by IRD. Finally, relative RV (v_⋆) with respect to this template (denoted by S(λ)) was measured by forward-modeling each observed spectral segment (f_obs(λ)) as

$$\begin{eqnarray}&&{f}_{\mathrm{obs}}(\lambda )=k(\lambda )\left[S\left(\lambda \sqrt{\displaystyle \frac{1+{v}_{\star }/c}{1-{v}_{\star }/c}}\right)T\left(\lambda \right)\right]\ast \mathrm{IP},\end{eqnarray} \tag{ 1 }$$

where T(λ) is the theoretical telluric model and k(λ) is a second-order polynomial describing the spectrum continuum for each segment (see Hirano et al. 2020, for more details). The resulting relative RVs are listed in Table 2. The typical RV uncertainty (internal error) was 4–7 m s⁻¹. We also estimated the absolute RV of K2-415 by fitting the individual molecular lines in the IRD spectra via Gaussian functions and comparing their centers with the vacuum line positions in the literature. Based on the analysis of 37 relatively deep OH lines in the H band, we obtained a mean absolute RV of 22.5 ± 0.1 km s⁻¹.

Using K2-415's template spectrum produced in Section 3.3, we performed a line-by-line analysis to estimate the stellar atmospheric parameters. We measured equivalent widths of FeH molecular lines and atomic lines of Na, Mg, Ca, Ti, Cr, Mn, Fe, and Sr to derive an effective temperature T_eff and abundances of individual elements [X/H]. For the T_eff estimation, 47 well-isolated FeH lines in the Wing–Ford band at 990–1020 nm were employed in the procedure described by Ishikawa et al. (2022). For the abundance analysis, we analyzed 33 atomic lines by following the procedure of Ishikawa et al. (2020).

We iterated the T_eff estimation and the abundance analysis alternately until both results became consistent with each other. First, we derived a provisional T_eff by adopting the solar metallicity as the initial value, and then by adopting this provisional T_eff, we derived the individual abundances of the eight elements. Subsequently, we adopted the iron abundance [Fe/H] derived in the previous step as the input metallicity of the T_eff estimation. Finally, adopting the new T_eff value, we redetermined the elemental abundances. The procedure up to this point allows the inputs and outputs of the analyses to be consistent within the measurement errors.

The final results of the T_eff and the elemental abundances are listed in Table 1. Note that the uncertainty of the T_eff here is calculated by the standard deviation of the estimates from the individual FeH lines, while it may also have a systematic error of less than 100 K as discussed in Ishikawa et al. (2022).

Based on those atmospheric parameters, we also derived the other physical parameters of K2-415. Inputting the above-derived metallicity [Fe/H] as well as 2MASS's K_s -band magnitude (Skrutskie et al. 2006) and Gaia parallax (Gaia Collaboration et al. 2021), we derived the stellar mass M_⋆ and radius R_⋆ using the empirical relations by Mann et al. (2019) and Mann et al. (2015), respectively. The uncertainties of those parameters were calculated with Monte Carlo simulations, assuming Gaussian distributions for the above input parameters and systematic errors in the empirical formula (Mann et al. 2019). In the Monte Carlo calculations, we derived the surface gravity $\mathrm{log}g$, mean stellar density ρ_⋆, and luminosity L_⋆. We also computed the Galactic (U, V, and W) velocities with respect to the Sun, using the Gaia DR3 information as well as the systemic RV (Table 1). All of those derived parameters are also summarized in Table 1. The low kinematic velocities indicate that K2-415 is a thin disk star.

Both K2 and TESS light curves in Figure 1 exhibit flux variations that are representative of spot-induced rotational modulations. To estimate the rotation period of K2-415, we computed the generalized Lomb–Scargle (LS) periodogram (Zechmeister & Kürster 2009) for the two light curves. As shown in Figure 5, both periodograms show strong peaks at similar periods, whose powers correspond to false-alarm probabilities (FAP) much lower than 10⁻³. Inspecting the peaks of the periodograms as well as the shapes of original light curves, we estimated the rotation period of the star to be P_rot = 4.36 ± 0.15 days for K2 and P_rot = 4.26 ± 0.12 days for TESS light curves, respectively. The uncertainties of these periods were derived based on the FWHMs of the peaks.

Zoom In Zoom Out Reset image size

The clear photometric variability (with the amplitude of 0.2%–0.4%) as well as the relatively short period of rotation suggests that K2-415 is a moderately active star. Inspecting an archived optical spectrum of K2-415 observed by the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Luo et al. 2015), we checked for its chromospheric activity. As shown in Figure 6, K2-415's low-resolution spectrum exhibits a moderate emission at Hα, whose EW is measured to be EWHα = −1.4600 ± 0.0036 Å ³⁴ (a negative value of EW indicates that it is an emission line), in agreement with the moderate variability seen in the light curves. We note that the variability amplitude of light-curve modulations seems to evolve rapidly in time, as evidenced by the vanishing variation in the second half of the TESS data. This fact suggests that activity-induced spectroscopic variations such as RV jitters do not produce a coherent pattern even on a timescale of 1–2 months. We will discuss the spot-induced photometric and spectroscopic variabilities in more detail in Section 4.4.

Zoom In Zoom Out Reset image size

Because two different sets of light curves (K2 and TESS) are available for K2-415, we first analyzed each data set independently, to check for consistency in the derived parameters (e.g., transit depth). For the transit analysis of K2 data, we ended up using the public light curve, K2SFF (Vanderburg & Johnson 2014), as K2SFF generally delivers better-quality light curves than our own ones (Section 2.1) in terms of flux scatters and behavior of outliers. Following the procedure in Hirano et al. (2015), we extracted light-curve segments around the transits (assuming a constant period, with a period based on a preliminary fit), and simultaneously fitted all the extracted segments while allowing for possible transit timing variations (TTVs). Each extracted segment covers approximately 5× the transit duration, thus typically involving 8–10 flux points. We implemented a Markov Chain Monte Carlo (MCMC) analysis to fit the light-curve segments after a global optimization of the fitting parameters using Powell's conjugate direction method as in Hirano et al. (2015). In the analysis, we employed the transit model by Ohta et al. (2009) and a linear function of time for the out-of-transit baseline. To take into account the ≈29 minute cadence of K2 photometry, for each observed flux we computed the transit model with one-minute sampling and binned the light curve to the K2 sampling. The fitting parameters are the scaled semimajor axis a/R_⋆, transit impact parameter b, star-to-planet radius ratio R_p /R_⋆, limb-darkening parameters in the quadratic law (u₁ and u₂), orbital period P, and mid-transit time for each transit ${T}_{c}^{(i)}$. We fixed the orbital eccentricity at e = 0 at this point. Only for the limb-darkening parameters, we imposed Gaussian priors as u₁ + u₂ = 0.81 ± 0.20 and u₁ − u₂ = 0.00 ± 0.20, based on the theoretical calculations by Claret et al. (2013). The result of this MCMC analysis is given in the leftmost column of Table 3. From the mid-transit times for individual transits, we also computed the zero-point mid-transit time T_c,0 for the K2 data set, assuming a linear ephemeris (i.e., a constant period).

Table 3. Results of Transit Analyses for K2-415b

Data Set	K2 (e = 0)	TESS (e = 0)	K2 + TESS (e = 0)	K2 + TESS (e prior)
(Fitting Parameter)
a/R⋆	${34}_{-16}^{+10}$	${27.4}_{-9.3}^{+4.2}$	${29.74}_{-0.91}^{+0.83}$	${29.63}_{-0.91}^{+0.86}$
b	${0.52}_{-0.36}^{+0.39}$	${0.43}_{-0.30}^{+0.37}$	${0.39}_{-0.18}^{+0.11}$	${0.37}_{-0.23}^{+0.19}$
Rp /R⋆	${0.0513}_{-0.0045}^{+0.0138}$	${0.0514}_{-0.0033}^{+0.0053}$	0.0470 ± 0.0018	${0.0472}_{-0.0018}^{+0.0020}$
${u}_{1}^{({\rm{K}}2)}+{u}_{2}^{({\rm{K}}2)}$	0.81 ± 0.20	⋯	0.79 ± 0.18	0.79 ± 0.19
${u}_{1}^{({\rm{K}}2)}-{u}_{2}^{({\rm{K}}2)}$	−0.01 ± 0.20	⋯	0.00 ± 0.20	−0.01 ± 0.20
${u}_{1}^{(\mathrm{TESS})}+{u}_{2}^{(\mathrm{TESS})}$	⋯	0.72 ± 0.20	0.74 ± 0.20	0.72 ± 0.19
${u}_{1}^{(\mathrm{TESS})}-{u}_{2}^{(\mathrm{TESS})}$	⋯	−0.17 ± 0.20	−0.16 ± 0.20	−0.17 ± 0.20
$\sqrt{e}\cos \omega$	0 (fixed)	0 (fixed)	0 (fixed)	${0.00}_{-0.30}^{+0.29}$
$\sqrt{e}\sin \omega$	0 (fixed)	0 (fixed)	0 (fixed)	${0.00}_{-0.25}^{+0.20}$
P (days)	4.01780 ± 0.00019	4.01803 ± 0.00024	4.0179682 ± 0.0000021	4.0179694 ± 0.0000027
Tc,0 − 2454833 (BJD)	3267.0727 ± 0.0018	4665.3230 ± 0.0033	3267.07137 ± 0.00054	3267.07116 ± 0.00069

Download table as:  ASCIITypeset image

We also analyzed the TESS light curve (PDCSAP FLUX), following the same procedure as above. The only differences in the modeling are that we did not integrate the theoretical transit model for TESS data (2 minutes cadence) and we adopted different priors for the limb-darkening parameters (u₁ + u₂ =0.73 ± 0.20 and u₁ − u₂ = − 0.17 ± 0.20 for the TESS band). The result of the MCMC fit to the TESS light curve is presented in the second column of Table 3.

The transit parameters (R_p /R_⋆ and P in particular) derived from K2 and TESS light curves are consistent with one another within 1σ, suggesting that no significant dilution is blended in those light curves, given that they were obtained with different apertures and different pass bands. The respective stellar densities estimated from the transit modeling are ${47}_{-39}^{+51}$ g cm⁻³ for K2 and ${24}_{-17}^{+13}$ g cm⁻³ for TESS light curves. Although these are both compatible with the mean stellar density derived in Section 4.1 (${\rho }_{\star }={30.3}_{-2.6}^{+2.9}$ g cm⁻³), the constraint is rather weak due to the degeneracy in a/R_⋆, b, and R_p /R_⋆ in the transit modeling. Figure 7 plots the observed minus calculated (O − C) residuals from the linear ephemeris for the transit times derived from the K2 (top) and TESS (bottom) data sets. No significant TTVs are seen for either data set; the χ² statistics for the transit times are 15.2 (K2) and 16.3 (TESS), with the degrees of freedom being 15 and 16, respectively. We also ensured that the two T_c,0 values derived from the K2 and TESS data sets are compatible with each other once P's error propagation is taken into account.

Zoom In Zoom Out Reset image size

Satisfied with the results of these consistency checks, we next performed joint analyses of K2 and TESS light curves. Based on the lack of TTVs, in the joint analyses we did not allow for a variable T_c for each transit, but instead allowed only the period P and global zero-point T_c,0 to float freely. In addition, to help break the degeneracy between a/R_⋆ and b, we imposed a Gaussian prior on the stellar density as ρ_⋆ =30.5 ± 2.7 g cm⁻¹, based on Table 1 in the MCMC analysis. We jointly modeled and fitted all the light-curve segments from the K2 and TESS data, and we derived the global posterior distribution for the fitting parameters.

In this joint analysis, we tried two different fits: one with e = 0 and the other with a floating e. The result of the e = 0 model fit is presented in the third column of Table 3. The derived parameters are in good agreement with the ones derived from fitting the K2 or TESS light curve alone, but with much smaller uncertainties. For the floating e model, we first attempted an MCMC analysis allowing $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$ to vary freely, but we found the fit did not converge, likely owing to the degeneracy in the fitting parameters together with the shallow transit and very short transit duration (see, e.g., Serrano et al. 2022). Orbital eccentricities of close-in planets have been studied in the past, and low or moderate eccentricities were suggested for small, close-in planets (e.g., Mayor et al. 2011; Van Eylen et al. 2019). Thus, we imposed weak Gaussian priors on $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$, with a mean and standard deviation of 0 and 0.335, respectively, so that the 1σ upper limit of e becomes 0.45 (Mayor et al. 2011). These priors were used only to account for the uncertainty in e and derive realistic errors for the other fitting parameters. The result of this MCMC fit for the e ≠ 0 model is shown in the rightmost column of Table 3. The phase-folded transit light curves after corrections for the out-of-transit baseline are plotted in Figure 8, along with the best-fit transit models.

Zoom In Zoom Out Reset image size

We computed the false-positive probability (FPP) of K2-415b using the Python package vespa (Morton 2015), which was developed for the statistical validation of planets from the Kepler mission (Morton et al. 2016). vespa employs a robust statistical framework to compare the likelihood of the planetary scenario to the likelihoods of several astrophysical false-positive scenarios involving eclipsing binaries, relying on simulated eclipsing populations based on the TRILEGAL galaxy model (Girardi et al. 2005). As inputs to vespa, we used an exclusion radius of 4'', WIYN/NESSI contrast curves (Section 3.2), phase-folded TESS photometry, and a 3σ upper limit on the secondary eclipse depth (Section 2.1), as well as broadband optical and NIR photometry, the Gaia parallax, and our spectroscopically derived estimates of the effective temperature and metallicity of the host star (Table 1). The FPP from vespa for K2-415b is 2 × 10⁻⁴, well below the commonly used validation threshold of 1%, which is in good agreement with the nondetection of any massive companions in our high-resolution imaging and RV measurements.

To confirm the low FPP, we also implemented another Python package, TRICERATOPS (Giacalone et al. 2021), which was developed to vet and validate TESS planet candidates. TRICERATOPS also returned a low FPP of 0.15% for K2-415b, confirming the result by vespa. K2-415b is therefore incompatible with any known false-positive scenarios, so we conclude that it has to be a real planet.

RV data presented in Table 2 show a relatively large scatter with a root-mean-square (rms) of 11.9 m s⁻¹, which is significantly larger than the mean RV internal error of 6.5 m s⁻¹. The expected RV variation by K2-415b is only a few m s⁻¹, based on the mass–radius relation (e.g., Otegi et al. 2020). In order to see if we can detect any planet signal associated with the period of K2-415b, we computed the generalized LS periodogram (Zechmeister & Kürster 2009) for the observed RV data. However, no significant peak showed up at P = 4.02 days, as drawn in Figure 9 (more generally, no peak exceeded the FAP = 1% threshold). While a fraction of the excess RV scatter might be ascribed to the instrumental systematics of IRD (e.g., Hirano et al. 2021; Harakawa et al. 2022), one should notice that K2-415 has a relatively short period of rotation (P_rot = 4.1–4.5 days, according to the LS periodograms) and both K2 and TESS light curves exhibit significant variations due to star spots. The stellar rotation period and the radius translate to the equatorial rotation velocity of ≈2.3 km s⁻¹, which places an upper limit on the $v\sin i$ of K2-415. Assuming that the system has a spin–orbit alignment (i.e., $v\sin i\approx 2.3$ km s⁻¹) and that the star has an effective surface-spot area of ≈0.2%–0.4% (from the light curves in Figure 1), one can roughly estimate the expected RV jitter amplitude as 8–16 m s⁻¹ (e.g., Desort et al. 2007; Boisse et al. 2012), which is comparable to the amplitude of the excess RV scatter. Therefore, in the following discussion, we model the observed RVs while taking into account the spot-induced RV jitters, to obtain an accurate constraint on the mass of K2-415b.

Zoom In Zoom Out Reset image size

To model the observed RVs, we employed the Gaussian Process (GP) based approach described in Rajpaul et al. (2015); in short, they applied the GP regression to the observed RVs together with the auxiliary parameters (the $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ index and the inverse slope of the cross-correlation bisector (BIS)) measured from the same spectra. To implement a similar GP regression, we measured some activity indices from the IRD spectra in a manner similar to that of Harakawa et al. (2022). Because the Ca HK line is not covered by IRD, instead of $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$, we adopted the FWHM of the mean line profile, which is known to display a behavior similar to that of $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ (Rajpaul et al. 2015). The mean line profile for each IRD spectrum is extracted using the least-squares deconvolution method (LSD; Donati et al. 1997; Asensio Ramos & Petit 2015), and in doing so, we used only OH lines, which are the most dominant opacity sources in the H-band spectra of early- to mid-M dwarfs. ³⁵ From the same LSD profile, we also measured the BIS. For both FWHM and BIS measurements, we subtracted the temporal mean values for individual lines and averaged over many different lines in the same spectrum in order to focus on the relative variations in FWHM and BIS, which we call ΔFWHM(t) and ΔBIS(t). More details of activity-index measurements for IRD spectra will be presented in our upcoming papers. The time series and LS periodograms of ΔFWHM(t) and ΔBIS(t) are also presented in Table 2 and Figure 9, respectively.

Following Rajpaul et al. (2015), we model the stochastic component of the observed RVs (ΔRV(t)) as well as ΔFWHM(t) and ΔBIS(t) by the following equations:

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{RV}(t)={V}_{c}G(t)+{V}_{r}\dot{G}(t),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{FWHM}(t)={WG}(t),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{BIS}(t)={B}_{c}G(t)+{B}_{r}\dot{G}(t),\end{eqnarray} \tag{ 4 }$$

where G(t) is a latent function associated with the fractional spot area and its derivative. The coefficients V_c , V_r , W, B_c , and B_r are determined by fitting the observed data. In a multidimensional GP regression, fitting parameters are optimized so that the logarithm of the likelihood function ${ \mathcal L }$ is maximized:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{ \mathcal L } & = & -\displaystyle \frac{1}{2}{\left({\boldsymbol{y}}-{\boldsymbol{m}}\right)}^{{\rm{T}}}{{\rm{\Sigma }}}^{-1}({\boldsymbol{y}}-{\boldsymbol{m}})\\ & & -\,\displaystyle \frac{1}{2}\mathrm{log}(| {\rm{\Sigma }}| )-\displaystyle \frac{n}{2}\mathrm{log}(2\pi ),\end{array}\end{eqnarray} \tag{ 5 }$$

where y and m are the vectors (with n components) representing the observed variables and mean functions, respectively, and Σ is the covariant matrix of the input variables. In our formulation, n = 3N_data, with N_data being the number of IRD spectra (= 42 in the present case). We adopt the mean functions m given by

$$\begin{eqnarray}&&{m}_{\mathrm{RV}}(t)=K\{\cos (f+\omega )+e\cos \omega \}+{\gamma }_{\mathrm{RV}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{m}_{{\rm{\Delta }}\mathrm{FWHM}}(t)={\gamma }_{{\rm{\Delta }}\mathrm{FWHM}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{m}_{{\rm{\Delta }}\mathrm{BIS}}(t)={\gamma }_{{\rm{\Delta }}\mathrm{BIS}},\end{eqnarray} \tag{ 8 }$$

where f is the true anomaly and K is the RV semi-amplitude due to K2-415b's Keplerian motion. The constant offset γ in each variable is also optimized in the fit. For the RV modeling, a fixed orbital period P is used, based on the transit light-curve analyses in Section 4.2.

The only remaining quantity to be determined a priori is the form of the covariant matrix Σ. Here, we adopt the quasi-periodic kernel for the covariance between the latent functions G(t_i ) and G(t_j ):

$$\begin{eqnarray}&&{k}_{\mathrm{qp}}({t}_{i},{t}_{j})\equiv \exp \left[-\displaystyle \frac{{\sin }^{2}\{\pi ({t}_{i}-{t}_{j})/{P}_{\mathrm{rot}}^{{\prime} }\}}{2{\lambda }_{{\rm{p}}}^{2}}-\displaystyle \frac{{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{\lambda }_{{\rm{e}}}^{2}}\right],\,\,\end{eqnarray} \tag{ 9 }$$

and compute the covariant matrix Σ following the expressions given by Rajpaul et al. (2015). The hyperparameters ${P}_{\mathrm{rot}}^{{\prime} }$, λ_p, and λ_e are also fitting parameters in the regression.

We attempted to implement a multidimensional GP regression to the observed quantities (RV, ΔFWHM(t), and ΔBIS(t)) using our custom MCMC code (Hirano et al. 2016), in which all the fitting parameters are allowed to float freely, but we immediately found that the fit did not converge, often resulting in walkers being stuck to different local minima. This is likely attributable to the small number of data points (only 42), spread over many years. In particular, the sparseness of the data seemed to prohibit reliable estimations of the periodicity (${P}_{\mathrm{rot}}^{{\prime} }$) and evolution timescale (λ_e) of RV jitters. We thus resorted to a two-step approach employed by Grunblatt et al. (2015), who used the light curves to constrain the hyperparameters in the covariance kernel (Equation (9)), and performed GP regressions to the K2 and TESS data before analyzing the spectral data. Adopting the covariance matrix for the flux values, whose elements are given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{ij}}={A}^{2}{k}_{\mathrm{qp}}({t}_{i},{t}_{j})+{\sigma }_{i}^{2}\delta ({t}_{i}-{t}_{j}),\end{eqnarray} \tag{ 10 }$$

where σ_i is the ith flux error, we ran MCMC analyses to estimate ${P}_{\mathrm{rot}}^{{\prime} }$, λ_p, and λ_e as well as the correlation amplitude A for each of K2 and TESS light curves (binned after masking the transits). The results of these analyses are presented in Table 4, and GP regressions to the light curves are plotted by the red solid lines in Figure 1. The derived parameters show a moderate disagreement (≈3σ except A) between K2 and TESS data, but these are likely due to different properties of the two data sets (e.g., different observing bands) as well as data processing. The GP-based period (${P}_{\mathrm{rot}}^{{\prime} }\approx 4.2\mbox{--}4.5$ days) roughly agrees with the rotation period derived from the LS periodograms (P_rot ≈ 4.3–4.4 days; see Section 4.1), which justifies the use of the quasi-periodic kernel for the GP regression.

Table 4. GP Regressions to the Light Curves

Data Set	K2	TESS
(Fitting Parameter)
A	${0.00288}_{-0.00031}^{+0.00039}$	${0.00132}_{-0.00013}^{+0.00017}$
${P}_{\mathrm{rot}}^{{\prime} }$ (days)	${4.164}_{-0.065}^{+0.084}$	${4.54}_{-0.11}^{+0.10}$
λp	${0.692}_{-0.058}^{+0.069}$	${0.411}_{-0.055}^{+0.069}$
λe (days)	4.71 ± 0.23	${3.65}_{-0.39}^{+0.41}$

Download table as:  ASCIITypeset image

We computed the weighted mean for each hyperparameter in Table 4, and used the results from the light-curve analysis as priors in GP regressions for the spectral data; employing Gaussian priors of ${P}_{\mathrm{rot}}^{{\prime} }=4.291\pm 0.061$ days, λ_p = 0.554 ± 0.044, and λ_e = 4.45 ± 0.20 days, we implemented the MCMC fit to the observed RVs along with ΔFWHM and ΔBIS. In the analysis, we introduced additional fitting parameters C_RV, C_ΔFWHM, and C_ΔBIS, to take into account underestimation/overestimation for internal errors of input variables, which are inserted into the diagonal components of the covariance matrix Σ as

$$\begin{eqnarray}&&{C}_{l}^{2}{\sigma }_{i}^{(l)2}\delta ({t}_{i}-{t}_{j}),\end{eqnarray} \tag{ 11 }$$

where l = {RV, ΔFWHM, ΔBIS} and ${\sigma }_{i}^{(l)}$ is the input statistical error at time t_i . We employed this formulation to model the white noise components in the fit, as opposed to adding extra terms of white noise in quadrature, because we found the input statistical errors of BIS we measured from the spectral analysis are slightly overestimated (i.e., C_ΔBIS < 1.0).

As in the case of transit photometry analyses, we tried two different fits with e = 0 and with priors on e. A fit with completely free $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$ did not result in a meaningful constraint on e, likely due to the small number of data points and significant RV jitters, and hence we imposed Gaussian priors on these parameters, as in Section 4.2, only to obtain realistic uncertainties for the other fitting parameters. With all the fitting parameters other than ${P}_{\mathrm{rot}}^{{\prime} }$, λ_p, λ_e, $\sqrt{e}\cos \omega$, and $\sqrt{e}\sin \omega$ being allowed to vary freely, we performed GP regressions to the spectral data by implementing MCMC analyses (Hirano et al. 2016). The results of these analyses are summarized in Table 5.

Table 5. GP Regressions to the IRD Data (RV + Spectral Indices)

Condition	e = 0	e prior
(Fitting Parameter)
K (m s−1)	${4.4}_{-3.3}^{+3.1}$	${4.1}_{-3.6}^{+3.5}$
$\sqrt{e}\cos \omega$ *	0 (fixed)	0.03 ± 0.29
$\sqrt{e}\sin \omega$ *	0 (fixed)	0.00 ± 0.32
γRV (m s−1)	1.6 ± 2.4	${1.6}_{-2.5}^{+2.4}$
γΔFWHM (km s−1)	−0.038 ± 0.018	−0.038±0.018
γΔBIS (km s−1)	${0.098}_{-0.042}^{+0.043}$	0.096 ± 0.042
${P}_{\mathrm{rot}}^{{\prime} }$ * (days)	${4.295}_{-0.060}^{+0.061}$	4.298 ± 0.061
λp *	0.564 ± 0.044	0.562 ± 0.044
λe * (days)	4.46 ± 0.20	4.46 ± 0.20
Vc (m s−1)	${5.0}_{-3.3}^{+3.4}$	5.2 ± 3.4
Vr (m s−1)	$-{2.3}_{-3.5}^{+4.2}$	$-{2.5}_{-3.5}^{+4.5}$
W (km s−1)	${0.008}_{-0.049}^{+0.041}$	${0.008}_{-0.048}^{+0.039}$
Bc (km s−1)	${0.00}_{-0.13}^{+0.12}$	$-{0.01}_{-0.13}^{+0.12}$
Br (km s−1)	${0.104}_{-0.052}^{+0.042}$	${0.102}_{-0.053}^{+0.043}$
CRV	${1.56}_{-0.24}^{+0.31}$	${1.55}_{-0.24}^{+0.30}$
CΔFWHM	${1.93}_{-0.25}^{+0.27}$	${1.94}_{-0.25}^{+0.27}$
CΔBIS	${0.61}_{-0.10}^{+0.15}$	${0.61}_{-0.10}^{+0.15}$

Note. Gaussian priors are imposed on the fitting parameters with asterisks.

Download table as:  ASCIITypeset image

Figure 10 plots the time series of the observed RVs along with the two spectral indices, and the red solid line for each panel indicates the best-fit GP regression (with the e prior) to the observed data. The phase-folded RV curve after removing the GP stochastic components is given in Figure 11. The derived RV semi-amplitude K of ${4.1}_{-3.6}^{+3.5}$ m s⁻¹ well agrees with that expected from the mass–radius relation (K = 1–2 m s⁻¹; Otegi et al. 2020), but it is also compatible with a nondetection of the planet within 1.1 σ, which is consistent with the lack of a significant peak at the orbital period in the LS periodogram of the RV data. Nonetheless, the 1 σ upper limit of K translates to the planet mass constraint of <5.7 M_⊕, which, along with the lack of a long-term RV trend, allows us to completely rule out FP with a stellar companion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size