Two Bright M Dwarfs Hosting Ultra-Short-Period Super-Earths with Earth-like Compositions^*

TESS observed TOI-1634 and TOI-1685 at a 2 minute cadence in Sectors 18 and 19, respectively. The observations were conducted from UT 2019 November 3 to 2019 December 23, resulting in photometry spanning approximately 27 days for each target, with gaps of about 1 day for data downlink in the middle of each observing sequence. Near the beginning of Sector 18, there is an additional 6.2 hr data gap due to the instrument being shut down for Earth eclipse. Light curves were produced by the Science Processing Operations Center (SPOC) photometry pipeline (Jenkins et al. 2016) using the apertures shown in Figure 1. For our transit analyses, we used the PDCSAP light curves produced by the SPOC pipeline (Smith et al. 2012; Stumpe et al. 2012, 2014). An error in the SPOC pipeline resulted in oversubtraction of the sky background, causing fractional changes (e.g., transits) in the light curves of TOI-1634 and TOI-1685 to be artificially deeper by 2.2% and 2.9%, respectively (Jon Jenkins, private communication). To account for this, we applied a correction to the R_p /R_⋆ values from our fits to the TESS data before combining them with our ground-based photometric measurements (see Section 3.2); we note the effect is smaller than the uncertainty of the R_p /R_⋆ values derived from the TESS light curves and has negligible impact on the final values. The SPOC pipeline applies a photometric dilution correction based on the CROWDSAP metric, which we independently confirmed by computing dilution values based on Gaia DR2 magnitudes (approximating G_RP as the TESS bandpass, and assuming an FWHM of 25''). For TOI-1634 there are two significantly contaminating sources in the aperture (Gaia DR2 IDs 223158499176634112 and 223158808416782208), which are 2.9 and 4.8 mag fainter in the G_RP band, respectively. For TOI-1685 there are three significantly contaminating sources (Gaia DR2 IDs 252366613254979328, 252366578895244672, and 252366578895245696), which are 3.7, 6.0, and 6.5 mag fainter in G_RP, respectively; an additional source (Gaia DR2 ID 252363589598010240) located just outside and to the south of the aperture is 0.25 mag brighter than TOI-1685 and thus also significantly contaminating despite contributing less than 10% of its flux.

Zoom In Zoom Out Reset image size

TOI-1634 has a resolved companion star separated by 25 from the primary star (Gaia DR2 ID 223158499176634112), with a TESS magnitude of 14.368 ± 0.010 mag (i.e., about 3.3 mag fainter than TOI-1634). The Gaia astrometry indicates the companion star has a parallax of 28.62 ± 0.11 mas and proper motions of ${\mu }_{\alpha }\cos \delta =80.64\pm 0.13$ mas yr⁻¹ and μ_δ = 14.539 ± 0.091 mas yr⁻¹, respectively (Gaia Collaboration et al. 2021), suggesting that TOI-1634 and the companion star share almost the same parallax and common proper motions. Thus, they are likely bound to each other, which was also reported in the visual-binary catalog for TOI's (Mugrauer & Michel 2020) as well as the more recent catalog by El-Badry et al. (2021) based on Gaia EDR3. Light-curve dilutions due to this companion star are taken into account when we perform the light-curve analyses. The impact of the companion on the estimation of the stellar properties as well as the long-term RV drift for TOI-1634 will be discussed in Sections 3.1 and 3.4. Other than this companion star, no stars were identified within ${1}^{{\prime} }$ in the Gaia EDR3 catalog having proper motions in common with TOI-1634 and TOI-1685.

The signature of TOI-1634.01 was initially detected by the TESS SPOC in a transiting-planet search of sector 18 that occurred UT on 2019 December 12, yielding a 1.8R_⊕ planet in a 0.98933 day orbit about its host star. The signal was detected at 10.6σ with an adaptive, noise-compensating matched filter (Jenkins 2002; Jenkins et al. 2010, 2020), passed all the diagnostic tests performed and published in the resulting Data Validation reports, and was fitted with a limb-darkened transit model (Twicken et al. 2018; Li et al. 2019). These included tests for eclipsing binaries, such as an odd/even depth test, a weak secondary test, and a ghost diagnostic test. The difference imaging centroid test showed that the source of the transit signature was consistent with the target star, TIC 201186294, with a measured offset from the target star of 81 ± 29 (we take 3σ as the confusion radius). The SPOC pipeline search removed the signature of TOI-1634.01 from the light curve and performed a search for additional transit signatures, which were not found. An alert for TOI-1634.01 was issued by the TESS Science Office (TSO) on UT 2020 January 14.

The signature of TOI-1685.01 was detected by the SPOC pipeline in a transiting-planet search of Sector 19 that occurred on UT 2020 January 17, resulting in a 1.47R_⊕ planet in a 0.6669 day orbit. This transit signature passed all the diagnostic tests performed and reported in the Data Validation reports archived to MAST and the TSO alerted the community to this planet candidate on UT 2020 January 30. The difference imaging centroid test showed that the source of the transit signature was consistent with the target star, TIC 28900646, with a measured offset from the target star of 279 ± 266. As was done for TOI-1634, the SPOC pipeline removed the signature of TOI-1685.01 from the light curve and performed a search for additional transit signatures, which were not found. We note that these difference imaging centroid measurements are complementary to the high-resolution imaging reported in Section 2.2, which is limited to separations of 12 and 30 from each target.

We independently confirmed the transit signals of each planet candidate using a second-order polynomial Savitzy–Golay filter to remove stellar variability and instrumental systematics from each light curve, then used the transit least-squares algorithm (TLS; Hippke & Heller 2019) ⁴⁵ to search them for transit signals, resulting in a signal detection efficiency (SDE) of 17.9, orbital period of 0.989 ± 0.003 days, and transit depth of 1.6 parts per thousand (ppt) for TOI-1634.01, and SDE of 18.6, orbital period of 0.669 ± 0.001 days, and transit depth of 1.0 ppt for TOI-1685.01. We subtracted each signal and repeated the transit search, but no additional transit signals with SDE above 10 were found in either light curve. TLS also reports the approximate depths of each individual transit; we note that these transit depths and uncertainties are useful for diagnostic purposes only, as they are simplistically determined from the mean and standard deviation of the in-transit flux. The depths of the odd transits are within 0.5σ of the even transits for both signals, suggesting a low probability of either signal being caused by an eclipsing binary at twice the detected period. These signals are consistent with those reported by the TESS team on ExoFOP-TESS. ⁴⁶ The TLS detections are shown in Figure 2.

Zoom In Zoom Out Reset image size

We observed four transits of TOI-1685.01 on UT 2020 November 24, UT 2021 January 10, UT 2021 January 12, and UT 2021 January 14 using the multiband imager MuSCAT (Narita et al. 2015) mounted on the 188 cm telescope at Okayama Astro-Complex in Japan. MuSCAT has three channels for the g, r, and z_s bands, enabling three-band simultaneous imaging observations. Each channel is equipped with a 1024 × 1024 pixel CCD camera with a pixel scale of 036 pixel⁻¹, which provides a field of view (FOV) of 61 square. We observed the target field with exposure times of 6–30 s depending on the band and sky condition. The obtained images were corrected for dark and flat in a standard manner, and aperture photometry was performed by a custom-built photometry pipeline (Fukui et al. 2011) to produce normalized light curves, in which the combinations of comparison stars and aperture radius were optimized such that the light-curve dispersion was minimized. The adopted aperture radius ranges from 8 to 14 pixels (from 29 to 51) depending on the band and night.

We observed five transits of TOI-1634.01 on UT 2020 February 7, UT 2020 February 10, UT 2020 February 11, UT 2021 February 14, and UT 2021 February 16 using the multiband imager MuSCAT2 (Narita et al. 2019) mounted on the 1.52 m TCS telescope at Teide Observatory in Spain. MuSCAT2 is a sibling of MuSCAT but has four channels for the g, r, i, and z_s bands. The CCD cameras of MuSCAT2 are identical to those of MuSCAT, but the pixel scale is 044 pixel⁻¹, which provides a $7\buildrel{\,\prime}\over{.} 4\times 7\buildrel{\,\prime}\over{.} 4$ FOV. The observations were conducted with exposure times of 3–60 s depending on the band and sky condition. The obtained data were reduced in the same way as for the MuSCAT data. We adopted aperture radii of 8–12 pixels (35–52) depending on the band and night, which means that the companion star at 25 away is contaminated into the photometric apertures in all bands.

We observed one transit of TOI-1685.01 on UT 2021 February 1 using the brand-new multiband imager MuSCAT3 (Narita et al. 2020), which was installed on the 2 m Faulkes Telescope North (FTN) at Haleakala Observatory in Hawaii in late 2020. The telescope and instrument are operated by Las Cumbres Observatory. As with MuSCAT2, MuSCAT3 has four channels for the g, r, i, and z_s bands, but has wider format CCD cameras with a size of 2k × 2k. The pixel scale of each camera is 0266 pixel⁻¹, which provides an FOV of 91 × 91. The observation was done slightly out of focus and with the exposure times of 25, 9, 8, and 20 s for the g, r, i, and z_s bands, respectively. The obtained raw images were processed by the BANZAI pipeline (McCully et al. 2018b) for dark and flat corrections, and then aperture photometry was performed in the same way as for the MuSCAT and MuSCAT2 data. The adopted radii of photometric aperture were 14, 18, 14, and 16 pixels (36, 47, 36, and 42) for the g, r, i, and z_s bands, respectively.

We observed a full transit of TOI-1634.01 on UT 2020 September 30 in Pan-STARRS z-short band and a full transit of TOI-1685.01 on UT 2020 November 11 in the Sloan $i^{\prime}$ band from the Las Cumbres Observatory Global Telescope (LCOGT) (Brown et al. 2013) 1.0 m network node at McDonald Observatory. We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen 2013), to schedule our transit observations. The 4096 × 4096 LCOGT SINISTRO cameras have an image scale of 0389 per pixel, resulting in a $26^{\prime} \times 26^{\prime}$ field of view. The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018a), and photometric data were extracted with AstroImageJ (Collins et al. 2017). The TOI-1634.01 observation was slightly defocused and used 40 s exposures and a photometric aperture radius of 58 to extract the differential photometry, resulting in a photometric precision of ∼500 ppm model residuals in 5 minute bins. The TOI-1685.01 observation was mildly defocused and used 50 s exposures and a photometric aperture radius of 78 to extract the differential photometry, resulting in a photometric precision of ∼410 ppm model residuals in 5 minute bins.

We observed a full transit of TOI-1685.01 at Observatoire du Mont-Mégantic, Canada, on UT 2020 March 8. The observations were made in the ${i}^{{\prime} }$ filter with a 15 s exposure time using the 1.6 m telescope of the observatory equipped with the 1024 × 1024 PESTO camera. PESTO has an image scale of 0466 per pixel, which provides an on-sky 795 × 795 FOV. The light-curve extraction via differential photometry was accomplished using an aperture radius of 70 and AstroImageJ. This software was also used for image calibration (bias subtraction and flat field division).

On UT 2020 December 2 TOI-1634 was observed with the 'Alopeke speckle imager (Scott 2019), mounted on the 8 m Gemini North telescope on Maunakea. 'Alopeke simultaneously acquires data in two bands centered at 562 nm 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 image achieved a contrast of Δmag = 8 at a separation of 1'' in the 832 nm band (see Figure 3). No secondary source was identified within 12 from TOI-1634.

Zoom In Zoom Out Reset image size

We observed TOI-1685 with near-infrared (IR) high-resolution adaptive optics (AO) imaging at the Keck Observatory. We carried out the AO imaging using the NIRC2 instrument on Keck II behind the natural guide star AO system. The observations were made on UT 2020 September 9 in the standard three-point dither pattern that is used with NIRC2 to avoid the left lower quadrant of the detector, which is typically noisier than the other three quadrants. The dither pattern step size was set to 3'' and was repeated twice, with each dither offset from the previous dither by 05.

The observations were made in the narrowband Br − γ filter (λ_o = 2.1686 μm; Δλ = 0.0326 μm) with an integration time of 1.5 s with one coadd per frame for a total of 13.5 s on target. The camera was in the narrow-angle mode with a full FOV of ≈ 10'' and a pixel scale of ≈ 000994 per pixel. The FWHM of the target in the combined image was ≈ 0052, and no additional stellar companions were detected in the 6'' × 6'' FOV (Figure 4).

Zoom In Zoom Out Reset image size

The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 20° at separations of integer multiples of the central source's FWHM. Following, e.g., David et al. (2019), we computed the 5σ sensitivity limit as a function of the radial distance from the target. The near-IR AO sensitivity curve for TOI-1685 is shown in Figure 4 along with an inset image zoomed to the primary target showing no other companion stars.

We obtained reconnaissance spectra of TOI-1634 on UT 2020 February 1 and UT 2020 September 4 and of TOI-1685 on UT 2020 February 2 and UT 2020 February 3 using the Tillinghast Reflector Echelle Spectrograph (TRES; Furesz 2008) located at the Fred Lawrence Whipple Observatory in Arizona, USA. TRES has a resolving power of ≈44,000 and a wavelength coverage of 385–910 nm, and the spectra were extracted as described in Buchhave et al. (2010).

RVs were determined from the TRES spectra using methods outlined in Winters et al. (2018). Briefly, molecular bands due to TiO in the wavelength range 7065–7165 Å found in aperture 41 of the TRES spectra were cross-correlated with an observed template spectrum of Barnard's Star (Gl 699). We conducted a search for maximum cross-correlation over a range of values of the rotational broadening $v\sin i$ applied to the template spectrum prior to correlation. As a result, we concluded there was no rotational broadening detectable in either target and therefore fixed the rotational broadening to zero for the final analysis. There is a systematic uncertainty in the velocity zero point of approximately 0.5 km s⁻¹, which may be important when considering the absolute barycentric RV, rather than the relative velocity differences between the epochs. We obtained RV = −17.066 km s⁻¹ (2020 February 1) and −17.105 km s⁻¹ (2020 September 4) for TOI-1634, and RV = −43.306 km s⁻¹ (2020 February 2) and −43.219 km s⁻¹ (2020 February 3) for TOI-1685. For each target, the two spectra were secured at near opposite quadratures in the orbital phase based on the TESS ephemerides. Therefore, the absence of large RV variations ( ≳ 0.5 km s⁻¹) ruled out stellar and brown-dwarf companions as the source of the transits for both targets.

For precise RV measurements of TOI-1634 and TOI-1685, we carried out near-IR observations of those two M dwarfs using Subaru/IRD between 2020 September and 2021 February under the Subaru IRD TESS intensive follow-up program (ID: S20B-088I). Every month during the period, we observed the two targets on two to three different nights when the program was assigned. On some of those nights, we visited the target stars twice within a night (two visits separated by a few hours) in order to mitigate the impact of the 1 day observing window, which happens to be close to the period of TOI-1634.01. IRD is a fiber-fed spectrograph placed in a temperature-stabilized chamber, which can simultaneously cover broadband near-IR wavelengths from 930 to 1740 nm with a spectral resolution of ≈ 70,000 (Tamura et al. 2012; Kotani et al. 2018). Stellar light collected by the telescope is first squeezed by the AO system on Subaru (Hayano et al. 2008), which is then injected into the spectrograph through a multimode fiber. For TOI-1634, the companion star at 25 was resolved in IRD's fiber injection module camera, and we ensured that only the primary (brighter) star was injected into the fiber. To trace the temporal instrumental stability, a secondary fiber is inserted into the spectrograph for simultaneous wavelength calibration, to which the laser-frequency comb (LFC) is usually injected. The integration times for both targets were set to 720–1200 s for each exposure, depending on the observing condition. We also observed at least one telluric standard star (A0 or A1 star) on each night to correct for the telluric lines in extracting the template spectrum for the RV analysis.

Raw IRD data were reduced by the standard procedure using IRAF (Tody 1993) as well as our custom codes to process the detector's bias and wavelength calibrations by LFC spectra (Kuzuhara et al. 2018; Hirano et al. 2020). The reduced one-dimensional spectra have a typical signal-to-noise ratio (S/N) of 60–95 per pixel at 1000 nm for both targets. Analyzing these reduced spectra, we extracted the RV for each frame. The RV analysis pipeline for IRD is described in Hirano et al. (2020); in short, individual observed spectra are first processed to create the stellar template spectrum, which is free from the telluric features and instrumental broadening. Using this stellar template as well as the instantaneous instrumental profile (IP) of the spectrograph (based on each LFC spectrum), each spectrum is fitted with the forward modeling technique. The typical RV internal errors are 3–4 m s⁻¹ for both targets.

To estimate the basic stellar parameters, we independently analyzed the optical high-resolution spectra taken by TRES and near-IR spectra by IRD. For the TRES spectra, we made use of SpecMatch-Emp (Yee et al. 2017) to determine the effective temperature T_eff, radius R_⋆, and iron abundance [Fe/H] of the stars. The code attempts to fit an observed (input) high-resolution spectrum to a number of library spectra, whose stellar parameters were well determined, and find the best-matched stars in the library, by which the stellar parameters for the input spectrum are determined by interpolations. SpecMatch-Emp returned T_eff = 3474 ± 70 K and 3468 ± 70 K, R_⋆ = 0.435 ± 0.044 R_⊙ and 0.417 ± 0.042 R_⊙, and [Fe/H] = 0.13 ± 0.12 dex and 0.03 ± 0.12 dex, for TOI-1634 and TOI-1685, respectively.

To estimate the atmospheric parameters for the two targets, we also analyzed the IRD spectra. Because many parts of the original IRD spectra suffer from significant telluric features (both absorptions and airglow emissions), we used the template spectra extracted for the RV analyses (Section 2.3), in which telluric features were removed and multiple frames were combined. The template spectra were then subjected to the analysis tool developed by Ishikawa et al. (2020). The analysis is based on a line-by-line comparison between the equivalent widths (EWs) from observed spectra and those from synthetic spectra. The synthetic spectra were calculated with a one-dimensional LTE spectral synthesis code that is based on the same assumptions as of the model atmosphere program of Tsuji (1978). For the atmospheric layer structure, we interpolated the grid of MARCS models (Gustafsson et al. 2008). The surface gravity $\mathrm{log}g$ and microturbulent velocity were needed to be assumed for the analysis. We referred to TIC for the $\mathrm{log}g$ values calculated from masses and radii (Stassun et al. 2019), which were estimated from the mass–M_K relation in Mann et al. (2019) and the radius–M_K relation in Mann et al. (2015), respectively. The microturbulent velocity was fixed at 0.5 ± 0.5 km s⁻¹ for both objects for simplicity.

First, we used the FeH molecular lines in the Wing–Ford band at 990–1020 nm for the T_eff estimation. The band consists of more than 1000 FeH lines, of which 57 lines with relatively clear line profiles were selected for the analysis. The adopted spectral line data are available from the MARCS web page. ⁴⁷ We measured the EW of each FeH line by fitting the Gaussian profile and found the T_eff at which the synthetic spectra best reproduce the EW by an iterative search. Throughout this first step, we assumed the solar value for the metallicity. The average of the T_eff estimates for each of the 57 lines was taken as the best estimate here. Its uncertainty was given as the line-to-line scatter calculated by the standard deviation over the estimates from all the lines. Those procedures will be provided in more detail in Ishikawa et al. (2021, in preparation).

As a second step, adopting the T_eff value estimated above, we determined the elemental abundances of Na, Mg, Si, Ca, Ti, Cr, Mn, and Fe from the corresponding atomic lines. The details of the abundance analysis are given in Ishikawa et al. (2020), although they adopted literature values for T_eff. The spectral line data were taken from the Vienna Atomic Line Database (Kupka et al. 1999; Ryabchikova et al. 2015). We selected the lines based on three criteria: (1) not suffering from blending of other absorption lines, (2) sensitive to elemental abundances, and (3) continuum level can be reasonably determined. The EWs were measured by fitting synthetic spectra on a line-by-line basis. We searched for an elemental abundance until the synthetic EW matches the observed one for each line and took the average for all the lines to estimate [X/H] for an element X.

Subsequently, we adopted the iron abundance [Fe/H] determined in the second step as the metallicity of the atmospheric model grid to redetermine the T_eff by the same procedure as in the first step. Then, we adopted the resulting T_eff to finally determine the elemental abundances including [Fe/H] again in the same way as in the second step. The procedure up to this point allows the results of T_eff and abundances to converge well within the measurement errors. Based on these analyses of IRD spectra, we obtained T_eff = 3432 ± 99 K and 3428 ± 97 K and [Fe/H] = 0.27 ± 0.12 dex and 0.27 ± 0.12 dex for TOI-1634 and TOI-1685, respectively. The abundances for the other elements are listed in Table 1.

We also performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia EDR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021) in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016), Stassun et al. (2017), Stassun et al. (2018). We pulled the JHK_S magnitudes from 2MASS (Skrutskie et al. 2006), the W1 – W4 magnitudes from WISE (Wright et al. 2010), the G, G_BP, G_RP magnitudes from Gaia (Gaia Collaboration et al. 2021), and the y-band magnitudes from Pan-STARRS (Flewelling et al. 2020). Together, the available photometry spans the full stellar SED over the wavelength range 0.4–20 μm (see Figure 5). We performed a fit using NextGen stellar atmosphere models, with T_eff and [Fe/H] as the free parameters; the extinction A_V was fixed at zero due to the proximity of the stars. Integrating the (unreddened) model SEDs gives the bolometric flux at Earth, F_bol. Finally, taking the F_bol and T_eff together with the Gaia parallax gives the stellar radius, R_⋆. The SED analysis provided T_eff = 3500 ± 85 K and 3475 ± 75 K, [Fe/H] = 0.0 ± 0.5 dex and 0.0 ± 0.5 dex, F_bol = (7.64 ± 0.27) × 10⁻¹⁰ erg s⁻¹ cm⁻² and (6.65 ± 0.15) × 10⁻¹⁰ erg s⁻¹ cm⁻², and R_⋆ = 0.466 ± 0.024 R_⊙ and 0.473 ± 0.021 R_⊙ for TOI-1634 and TOI-1685, respectively.

Zoom In Zoom Out Reset image size

The three measurements (optical spectroscopy, near-IR spectroscopy, and SED fitting) of T_eff and [Fe/H] yielded consistent results within their errors, and thus we computed the weighted means of those parameters to gain the final values (Table 1) used in the subsequent analyses. Because these measurements ultimately rely on similar stellar atmosphere models or the same calibration sources, we conservatively adopted the representative errors for the mean values of the two parameters (i.e., 70 K for T_eff and 0.12 dex for [Fe/H]). Based on the basic parameters derived above, we further estimated the other stellar parameters (i.e., the stellar mass M_⋆, radius R_⋆, surface gravity $\mathrm{log}g$, mean density ρ_⋆, and luminosity L_⋆), as well as refined the basic parameters (i.e., the stellar metallicity [Fe/H] and distance d) by combining all observed quantities in a consistent manner. In doing so, we took an approach described in Hirano et al. (2018), but with the inclusion of Gaia parallaxes; because the observed quantities are redundant (e.g., there are two sets of estimates for the stellar radius) and can be correlated with each other through the empirical relations, we performed Markov Chain Monte Carlo (MCMC) simulations in which the χ² statistic of the likelihood function ($\propto \exp (-{\chi }^{2}/2)$) is defined as

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \displaystyle \frac{{({R}_{\star ,\mathrm{TRES}}-{R}_{\star })}^{2}}{{\sigma }_{{R}_{\star },\mathrm{TRES}}^{2}}+\displaystyle \frac{{({R}_{\star ,\mathrm{SED}}-{R}_{\star })}^{2}}{{\sigma }_{{R}_{\star },\mathrm{SED}}^{2}}\\ & & +\displaystyle \frac{{({m}_{{K}_{s},2\mathrm{MASS}}-{m}_{{K}_{s}})}^{2}}{{\sigma }_{{m}_{{K}_{s}},2\mathrm{MASS}}^{2}},\end{array}\end{eqnarray} \tag{ 1 }$$

where ${R}_{\star ,\mathrm{TRES}}$ and R_⋆,SED, are the stellar radii estimated by the optical spectroscopy and SED integration, and ${\sigma }_{{R}_{\star },\mathrm{TRES}}$ and ${\sigma }_{{R}_{\star },\mathrm{SED}}$ are their errors, respectively. The apparent K_s -band magnitude by 2MASS and its error are denoted by ${m}_{{K}_{s},2\mathrm{MASS}}$ and ${\sigma }_{{m}_{{K}_{s}},2\mathrm{MASS}}$, respectively. The fitting parameters in the MCMC analysis are the absolute K_s magnitude ${M}_{{K}_{s}}$, stellar metallicity [Fe/H], and the distance d to the system. The modeled quantities R_⋆ and ${m}_{{K}_{s}}$ in the right-hand side of Equation (1) are calculated from ${M}_{{K}_{s}}$, [Fe/H], and d through the empirical relation by Mann et al. (2015) and ${m}_{{K}_{s}}-{M}_{{K}_{s}}=5.0\mathrm{log}d\,-5.0$. We assume A_V = 0, given the proximity of the two stars to Earth. We imposed Gaussian priors on [Fe/H] and d based on the weighted mean value and its error for [Fe/H] derived above, and the Gaia parallax (Gaia Collaboration et al. 2021). In implementing the MCMC analysis, we computed M_⋆ via the empirical relation of Mann et al. (2019) from ${M}_{{K}_{s}}$ and [Fe/H], as well as the surface gravity $\mathrm{log}g$, the mean density ρ_⋆, and the luminosity L_⋆ for each step of the chain. For L_⋆, we sampled the T_eff values with the Gaussian distribution based on the values in Table 1.

TOI-1634 has a companion star 25 away from the primary star, but we were unable to identify the companion star in the 2MASS catalog. We inspected the 2MASS image for TOI-1634 and found that the companion star was buried in the point-spread function of the primary star, whose FWHM was found to be 27–28). This suggests that the K_s magnitude listed in Table 1 may be contaminated by the companion star, and the true magnitude of the primary star could be slightly fainter. To roughly estimate its impact, we used the Dartmouth isochrone model (Dotter et al. 2008) and inferred the mass of the companion. Because the Dartmouth isochrones list the Gaia magnitudes as a function of stellar mass for a given set of stellar age and metallicity, we employed the Gaia G_RP magnitude to constrain the companion's mass. The magnitude difference of ΔG_RP = 2.959 between TOI-1634 and its companion translates to the companion's mass of ≈ 0.12 M_⊙ on the assumption that TOI 1634's mass is roughly ≈ 0.46 M_⊙. When those masses are adopted, the isochrones predict that the magnitude difference in the K_s band should be ${\rm{\Delta }}{m}_{{K}_{s}}\approx 2.8\mbox{--}3.0$ mag, implying that the true ${m}_{{K}_{s}}$ of the primary star is ≈ 0.07 mag fainter than the reported one. With this in mind, we adopted ${m}_{{K}_{s}}=8.67\pm 0.07$ instead of ${m}_{{K}_{s}}=8.600\pm 0.014$ for TOI-1634 (in addition to shifting the center value of the magnitude, we conservatively added the systematic error of 0.07 in ${m}_{{K}_{s}}$ in quadrature) and ran the MCMC analysis. For TOI-1685, we directly input the 2MASS K_s magnitude in the code. MCMC simulations were implemented using our custom code (e.g., Hirano et al. 2015) with a chain length of 10⁶ after the burn-in chains. The final derived parameters based on this MCMC analysis (d, [Fe/H], $\mathrm{log}g$, M_⋆, R_⋆, ρ_⋆, and L_⋆) are summarized in Table 1.

Using the Gaia EDR3 information as well as the absolute RVs from the TRES spectra, we also computed the Galactic space velocities (U, V, W) for the two stars with respect to the Sun (Table 1). The low space velocities for both targets indicate those stars belong to the thin disk. Velocity dispersions in the Galactic coordinate system are generally correlated with stellar age. Following the methodology described in Burgasser & Mamajek (2017), we computed the posterior distributions for the ages of the two stars. In doing so, we adopted the prescription given by Sanders & Binney (2015) for the velocity-dispersion evolution of the thin-disk stars with the Sun's peculiar velocity from Bland-Hawthorn & Gerhard (2016), and we used two different age priors: a uniform prior (0 < age ≤ 14 Gyr) and the age probability distribution in the Geneva–Copenhagen Survey (GCS) catalog (Casagrande et al. 2011). Based on the age posterior distributions, we found TOI-1634 has the age of ${3.2}_{-2.8}^{+6.2}$ Gyr (uniform prior) and ${5.2}_{-2.8}^{+4.0}$ Gyr (GCS prior) and that of TOI-1685 is ${5.0}_{-3.7}^{+5.6}$ Gyr (uniform prior) and ${5.7}_{-3.0}^{+3.8}$ Gyr (GCS prior), respectively. These results suggest the UVW velocities are not useful for constraining the ages of the two targets. We also confirmed that neither of the targets belong to nearby young associations based on the BANYAN Σ tool (Gagné et al. 2018).

The planetary transit was securely detected in the light curves by the ground-based photometry (Figure 8), in which the observed transit depths were consistent with the TESS photometry. However, the companion star at 25 was inside the photometric aperture ⁴⁹ , meaning that the ground-based photometry alone was not capable of ruling out the possibility that the transits are originating from the companion star (companion's flux contamination is larger than the transit depth). In order to check if our RV data alone indicate the presence of the USP planet around TOI-1634, we performed a period analysis using the GLS tool (Zechmeister & Kürster 2009) applied to the observed IRD-RV data. The upper panel of Figure 13 shows the GLS periodogram for TOI-1634's raw RV data. There are multiple peaks with very low FAPs (<0.1%), but the highest peak shows up at the correct period of the transiting planet (P = 0.989 days), which does not fall on the peaks of the window function (blue-shaded area). Therefore, our RV data indicate additional, independent evidence of the USP planet orbiting TOI-1634 and not orbiting its companion star.

Zoom In Zoom Out Reset image size

Next, we attempted the orbital fit to the observed RVs. In doing so, we first estimated the impact of the companion star around TOI-1634; given the proximity to the star, the companion star at 25 away might have a nonnegligible impact on the long-term RV baseline. With the distance of d = 35 pc for TOI-1634, the angular separation of 25 translates to the projected separation of 88 au, which approximately sets the lower limit to the semimajor axis of the binary orbit except for a highly eccentric orbit (i.e., a_binary ≳ 88 au). The RV acceleration of the primary star ($\dot{\gamma }$) around the center of mass of the system is expressed as

$$\begin{eqnarray}&&\dot{\gamma }=\displaystyle \frac{{{GM}}_{\mathrm{comp}}}{{a}_{\mathrm{binary}}^{2}}\sin {i}_{o}\sin (f+\omega ){\left(\displaystyle \frac{1+e\cos f}{1-{e}^{2}}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant, M_comp is the companion star's mass, i_o is the orbital inclination, f is the true anomaly, e is the orbital eccentricity, and ω is the argument of periastron. When we assume the companion's mass of ≈0.1 M_⊙ (see Section 3.1) and e = 0 for the binary orbit, the lower limit on a_binary gives the maximum RV acceleration as

$$\begin{eqnarray}&&| \dot{\gamma }| \leqslant \displaystyle \frac{{{GM}}_{\mathrm{comp}}}{{a}_{\mathrm{binary}}^{2}}\lesssim 7.3\times {10}^{-3}\,{\rm{m}}\,{{\rm{s}}}^{-1}\,{\mathrm{day}}^{-1}.\end{eqnarray} \tag{ 3 }$$

In the presence of moderate eccentricity, the upper limit of $\dot{\gamma }$ could be a few times larger than the above value, depending on the orbital phase. Hence, this order-of-magnitude estimation suggests that the stellar companion may lead to an RV drift of up to a few m s⁻¹ over the course of ≈5 months.

We constrained the visual-binary orbital parameters using the LOFTI_gaiaDR2 software package (Pearce et al. 2020). LOFTI_gaiaDR2 uses the instantaneous positions, proper motions, and masses of the components of visual-binary stars to estimate their orbital parameters. We used the astrometric parameters from Gaia EDR3 for this calculation, along with the stellar masses for the primary and secondary stars estimated in Section 3.1. The LOFTI_gaiaDR2 posterior probability distribution has a slight preference for highly eccentric solutions (68% confidence interval between e = 0.61 and 0.98) but remains consistent with circular orbits. We note that these parameters should be taken with some skepticism because the astrometric solution for the secondary star shows excess scatter (with a Renormalized Unit Weight Error, or RUWE, of 1.7), which can indicate that it is itself an unresolved binary companion that can significantly affect its proper motion. Regardless, we conclude that the Gaia positions and proper motions are not inconsistent with an eccentric visual binary orbit.

Based on these speculations, we modeled the observed RVs of TOI-1634 using the following equation, in which we allow for the presence of a possible RV trend:

$$\begin{eqnarray}&&\mathrm{RV}(t)=K[\cos (f+\omega )+e\cos \omega ]+\gamma +\dot{\gamma }(t-{t}_{0}),\end{eqnarray} \tag{ 4 }$$

where K is the RV semiamplitude and γ is the RV offset of our data set. The time t₀ is an arbitrary origin of time for which we adopt the time of the first RV point in the whole data set. We optimized the orbital parameters (K, $e\cos \omega$, $e\sin \omega$, γ, $\dot{\gamma }$) using MCMC (Hirano et al. 2015) with uniform priors for all parameters. In the fit, we fixed P and T_c based on the transit ephemeris (Table 2).

We attempted the orbital fits assuming both circular and eccentric orbits. The results of those fits are shown in Table 3 ("with $\dot{\gamma }$" columns). To discuss the significance of the nonzero eccentricity, we compared the Bayesian Information Criterion (BIC), which is computed by $\mathrm{BIC}={\chi }_{\mathrm{best}}^{2}+k\mathrm{ln}(N)$, where k is the number of fitting parameters and N is the number of data points. Comparing the two BIC values for the above solutions, we found ΔBIC = BIC_e=0 − BIC_e≠0 = 1.5, implying that the circular and eccentric orbital solutions are almost equally favored. In other words, no evidence for nonzero eccentricity is found in our data set. A near-zero orbital eccentricity is also expected from the tidal circularization timescale for USP planets; using Equation (17) of Patra et al. (2017) with the planetary tidal quality factor of Q_p ≈100 (for a terrestrial planet) (e.g., Ment et al. 2021), we obtain the tidal damping timescale of ≈5.5 × 10⁴ yr for TOI-1634b, implying that a nonzero eccentricity should have been damped in the past. Therefore, we concluded that the TOI-1634b has an almost circular orbit and adopt the fitting result for e = 0 in the subsequent analysis. The RV data and the best-fit orbital solution to the data are plotted in panels (a) and (b) of Figure 14.

Zoom In Zoom Out Reset image size

Table 3. Results of the Orbital Fits

	TOI-1634b			TOI-1685b
Parameter	with $\dot{\gamma }$ (e = 0)⋆	with $\dot{\gamma }$ (e ≠ 0)	no $\dot{\gamma }$ (e = 0)	with Feb-02 (e = 0)	no Feb-02 (e = 0)⋆	no Feb-02 (e ≠ 0)
K ( m s−1)	11.1 ± 1.0	11.31 ± 0.99	11.80 ± 0.91	3.3 ± 1.1	4.2 ± 1.1	${7.0}_{-1.6}^{+1.5}$
$e\cos \omega$	0 (fixed)	$-{0.118}_{-0.038}^{+0.040}$	0 (fixed)	0 (fixed)	0 (fixed)	0.278 ± 0.076
$e\sin \omega$	0 (fixed)	${0.010}_{-0.077}^{+0.075}$	0 (fixed)	0 (fixed)	0 (fixed)	${0.03}_{-0.21}^{+0.20}$
$\dot{\gamma }$ (m s−1 day−1)	−0.023 ± 0.012	−0.023 ± 0.012	0 (fixed)	0 (fixed)	0 (fixed)	0 (fixed)
BIC	70.6	69.1	70.0	54.4	43.1	39.0

Note. For each planet, the fitting result adopted to compute the planet mass is indicated by ⋆.

Download table as:  ASCIITypeset image

The best-fit RV acceleration $\dot{\gamma }$ is ≈3 times larger than the value in the right-hand side of Equation (3), but it is consistent with zero within 2σ. While this possibly large RV drift might be attributed to a moderate eccentricity of the binary orbit as discussed above, it could be an artifact caused by a small number of RV points around the beginning and/or end of our observing campaign spanning ∼5 months. Given the frequency of planet multiplicity for USP planets (Winn et al. 2018), it is also possible that there exists an outer planet in the system that gave systematic offsets at specific orbital phases for the inner USP planet. To discuss the significance of this RV trend, we next fitted the observed RV data in the absence of the RV trend $\dot{\gamma }$ assuming a circular orbit. Our MCMC analysis suggested K = 11.80 ± 0.91 m s⁻¹, which is compatible with the result in the presence of $\dot{\gamma }$. Comparing the BICs for the two fitting results, we found that the result without the trend is equally likely (${\rm{\Delta }}\mathrm{BIC}={\mathrm{BIC}}_{\mathrm{with}\,\dot{\gamma }}-{\mathrm{BIC}}_{\mathrm{no}\,\dot{\gamma }}=0.6$). We thus list both fitting results (with and without $\dot{\gamma }$) in Table 3 to take into account the uncertainty of the systematic RV offset. We employ the result with $\dot{\gamma }$ and e = 0, which is physically motivated from the dynamics of the system, in deriving the planet mass M_p as well as the mean density ρ_p from K (Table 2).

After removing the best-fit orbital model (e = 0, $\dot{\gamma }\ne 0$) for the observed RV data, we performed an extra periodogram analysis to search for additional planets in the system. The bottom panel of Figure 13 illustrates the GLS periodogram (red solid line) for the residual RV data. No significant peak was found in the residual RVs, suggesting either that no additional massive planet is present in the system with a period shorter than our observation span or that the signal of such unidentified planets was removed/minimized by the orbital fit of TOI-1634b and long-term RV trend. At this point, our RV data imply no evidence for additional planets in the system.

We ran a period analysis for the observed RV of TOI-1685 in a similar manner to TOI-1634. The upper panel of Figure 15 plots the GLS periodogram for the raw RV data. There are several significant peaks exceeding the FAP = 0.1% line, but the one at the period of TOI-1685b (P = 0.669 days) is not high enough to claim the detection of the orbital signal. After a preliminary orbital fit to the observed RV data using the transit ephemeris, we found that the RV points taken on UT 2021 February 2 (hereafter "Feb-02") are the primary outliers, deteriorating the fitting result for the planet. Although this could be indicative of the presence of an additional planet in the system, we also suspected that this sudden RV shift is caused by an instrumental systematic. The IRD spectrograph is known to exhibit a relatively large temporal RV drift, which is well correlated with the temperature instability at the camera lens inside the chamber (Kotani et al. 2018; Hirano et al. 2020). This instrumental RV drift is usually corrected by modeling the instantaneous IP of the spectrograph derived from the simultaneously taken wavelength-reference spectrum (i.e., LFC). However, if the variation in the IP is too fast compared to each integration time, it is theoretically expected that the LFC is unable to accurately trace the "effective" instantaneous IP of the spectrograph.

Zoom In Zoom Out Reset image size

To further investigate this possibility, we inspected the absolute instrumental drift of the spectrograph on February 2 and found that the IRD spectrograph indeed exhibits a large instrumental instability that night as shown in Figure 16. In particular, TOI-1685 was observed at the very beginning of the night (blue squares), when the instrumental RV variation was most significant; the spectrograph exhibits an RV drift of ≈8 m s⁻¹ for every integration. ⁵⁰ In addition, the observing condition during the twilight usually changes dramatically, and thus the combination of instrumental instability and variations in the twilight observing conditions may have affected the extraction and application of the effective IPs from the LFC spectra.

Zoom In Zoom Out Reset image size

The impact of IRD's instrumental RV drift, especially for the case of relatively long integrations, is under investigation, and therefore we decided to perform the orbital fits with and without including the Feb-02 data. We first computed the periodogram for the data set excluding the Feb-02 data. The lower panel of Figure 15 plots the resulting GLS periodogram. While the same peaks (FAP < 0.1%) identified for the original data set (upper panel) have similar GLS powers, the peak at the correct period of TOI-1685b (P = 0.669 days) now appears with a low FAP (<0.1%); whether instrumental or astrophysical, the absence of a significant peak at TOI-1685b's orbital period in the original periodogram is ascribed to the inclusion of the Feb-02 data. The two peaks around 0.70 and 0.72 days in Figure 15, which are higher than the 0.67 day peak, are likely aliases associated with the peak at 2.59 days. The window function has peaks at 1.0 day and 0.96 days (the highest and second-highest ones for P<10 days). When those window frequencies are coupled with the period at 2.59 days, the periodogram would exhibit alias peaks around 0.72 and 0.70 days, respectively. The 2.59 day periodicity will be discussed later.

For the RV data with and without the Feb-02 data, we next fitted the observed RVs with a single-planet model. Assuming either a circular or eccentric orbit, we performed the MCMC analysis as in the case of TOI-1634 for each data set. When the Feb-02 data were included, we obtained K = 3.3 ± 1.1 m s⁻¹ and $K={4.1}_{-1.4}^{+1.5}$ m s⁻¹ for the circular and eccentric orbits, respectively. The two fitting results yielded ΔBIC = BIC_e=0 − BIC_e≠0 = −2.3, implying that the circular orbit is slightly favored for this data set. We obtained larger K values in the absence of the Feb-02 data: K = 4.2 ± 1.1 m s⁻¹ and $K={7.0}_{-1.6}^{+1.5}$ m s⁻¹ for e = 0 and e ≠ 0, respectively. In this case, the two fits resulted in ΔBIC = BIC_e=0 − BIC_e≠0 = +4.0; unlike the case with the Feb-02 data, an eccentric orbit is a slightly favorable solution. Note that as in the case of TOI-1634, the tidal circularization timescale for TOI-1685.01 is estimated as ≈1.0 × 10⁴ yr for Q_p ≈100 (Earth-like rocky planet), indicating that e should be vanishingly low in the absence of an additional planet in the system. Those fitting results are shown in Table 3 and the phase-folded RVs are plotted in Figure 17. For the final planet mass M_p (Table 2), we adopt the K value for the case of e = 0 without the Feb-02 data.

Zoom In Zoom Out Reset image size

In order to search for an additional signal in the observed RV data, we computed the periodogram for TOI-1685's RVs after removing the best-fit single-planet model for each data set. Considering the short tidal circularization timescale for the USP planet, we removed the circular-orbit solutions derived above. Figure 18 plots the GLS periodograms for the whole RV data and the data subset without the Feb-02 data. For both panels, there are a few significant peaks (FAP<0.1%) that do not fall in the window function. The peak at 2.6 days is common to both periodograms, which was also seen in the original RVs without the Feb-02 data (lower panel of Figure 15). The high peaks at P<1.0 day in both panels are likely alias peaks associated with the 2.6 day peak and window functions.

Zoom In Zoom Out Reset image size

Given the limited phase coverage and unknown instrumental systematics, at this point we are not able to claim that the 2.6 day periodicity in the RV data represents an additional planet in the system; more RV measurements are essentially required to gain a robust conclusion on the presence of an additional body in the system. Nonetheless, we were tempted to fit the observed RV data with a two-planet model. In doing so, we ran the MCMC code and fitted the RV data (with and without the Feb-02 data) assuming two circular Keplerian orbits. We fixed the period of the USP planet at the one from the transit ephemeris and allowed the period of the outer planet P₂ and time of the inferior conjunction T_c,2 to float with uniform priors. The results of these fits are listed in Table 4. In the table, K₁ and K₂ represent the RV semiamplitudes for the inner (USP) and outer planets, respectively. The phase-folded RV curves (with the inclusion of Feb-02 data) after removing the Keplerian orbit for the other planet are shown in Figure 19. The RV semiamplitudes for the USP planet are consistent within ≈1σ with the values derived for the one-planet model (Table 3) in both cases, whereas the RV scatters around the best-fit models significantly improved with ΔBIC = BIC_one−planet − BIC_two−planet being greater than 10 for both fits.

Zoom In Zoom Out Reset image size

Table 4. Results of the Orbital Fits for TOI-1685 with a Two-planet Model

Parameter	with Feb-02 Data	no Feb-02 Data
K1 ( m s−1)	${4.8}_{-1.2}^{+1.1}$	${4.9}_{-1.5}^{+1.3}$
K2 ( m s−1)	6.2 ± 1.0	${5.6}_{-1.0}^{+1.0}$
P2 (days)	${2.5909}_{-0.0048}^{+0.0045}$	${2.5891}_{-0.0069}^{+0.0054}$
Tc,2 (BJDTDB)	${2458827.24}_{-0.64}^{+0.67}$	${2458827.46}_{-0.75}^{+0.90}$

Download table as:  ASCIITypeset image

We note that the 2.6 day signal is unlikely to be explained by stellar rotation. If the rotation period of the star is P_rot = 2.6 days, the equatorial rotation velocity must be ≈8.8 km s⁻¹, which also gives the projected rotation velocity $v\sin {i}_{s}$ for the case of spin–orbit alignment in the system. Both TRES optical spectra and IRD near-IR spectra, however, imply that the star is slowly rotating with $v\sin {i}_{s}\lt 5$ km s⁻¹. The slow rotation of TOI-1685 is also supported by the low-frequency light-curve modulation as discussed in Section 3.2. Therefore, we conclude that the 2.6 day periodicity does not indicate the rotational signal in the RV data, but represents any one of (1) an additional planet, (2) an instrumental/analysis artifact, or (3) an artifact caused by the mixture of (1) and (2) as well as the window function of our IRD observations. Again, further observations are required to test those possibilities.

If the 2.6 day signal indeed represents the period of the outer planet, K₂ ≈6 m s⁻¹ corresponds to the planetary mass of $M\sin {i}_{o}=7\mbox{--}8\,{M}_{\oplus }$. Although the two planets in the system have relatively small masses, the small orbital separation between the two planets prompted us to check for the orbital stability of the two planets. Because the outer one is not transiting and its orbital inclination (thus the true mass) is not known, currently there is little point in running detailed numerical simulations for the system. Instead, we simply compared the minimum separation between the two in terms of the mutual Hill sphere R_H, following Pu & Wu (2015). Inputting the semimajor axes of the two planets (a₁ = 0.0116 au and a₂ = 0.0285 au for the inner and outer planets, respectively) on the assumption that the planets are coplanar, we found R_H ≈0.00058 au. Thus, the minimum separation between the two planets (a₂ − a₁ = 0.0169 au) is about 29 times larger than the mutual Hill radius. Pu & Wu (2015) showed that if the minimum separation is larger than ≈12 R_H, the system should be stable on a billion-year timescale. Also considering that the periods of the two planets are not near a first-order mean-motion resonance, the addition of a super-Earth-mass planet at P = 2.6 days does not critically deteriorate the stability of the system.