TOI-1696: A Nearby M4 Dwarf with a 3 R_⊕ Planet in the Neptunian Desert

TESS observed TOI-1696 with a 2 minute cadence in Sector 19 from 2019 July 25 to August 22, resulting in photometry spanning approximately 27 days with a gap of about one day in the middle when the satellite reoriented itself for data downlink near perigee. Light curves were produced by the Science Processing Operations Center (SPOC) photometry pipeline (Jenkins 2002; Jenkins et al. 2010, 2016) using the aperture shown in Figure 1. We used the PDCSAP light curves produced by the SPOC pipeline (Stumpe et al. 2012; Smith et al. 2012; Stumpe et al. 2014) for our transit analyses. TOI-1696 is located in a fairly crowded field, owing to its low galactic latitude (b = −081). The SPOC pipeline applies a photometric dilution correction based on the CROWDSAP metric, which we independently confirmed by computing dilution values based on Gaia Data Release 2 (DR2) magnitudes. ³⁵

Zoom In Zoom Out Reset image size

TOI-1696.01 was detected by the SPOC pipeline in a transiting planet search, and the candidate was subsequently reported to the community by the TESS Science Office (TSO) on 2020 January 30 via the TESS Object of Interest (TOI; Guerrero et al. 2021) Releases portal. ³⁶ The candidate passed all data validation diagnostic tests (Twicken et al. 2018) performed by the SPOC. ³⁷ The SPOC pipeline removed the transit signals of TOI-1696.01 from the light curve and performed a search for additional planet candidates (Li et al. 2019), but none were reported.

We independently confirmed the transit signal found by the SPOC. After removing stellar variability and residual instrumental systematics from the PDCSAP light curve using a second-order polynomial Savitzky–Golay filter with a window =1001, we searched for periodic transit-like signals using the transit least-squares algorithm (TLS; Hippke & Heller 2019), ³⁸ resulting in the detection of TOI-1696.01 with a signal detection efficiency (SDE) of 11.6, a transit signal-to-noise ratio (S/N) of 7.4, orbital period of 2.50031 ± 0.00001 days, and transit depth of 10.6 parts per thousand (ppt), which is consistent with the values reported by the TESS team on the Exoplanet Follow-up Observing Program for TESS (ExoFOP-TES; ³⁹ ExoFOP 2019). We subtracted this signal and repeated the transit search, but no additional signals with SDE above 10 were found. 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 1.5σ of the even transits, suggesting a low probability of either signal being caused by an eclipsing binary at twice the detected period. The TLS detection is shown in Figure 2.

Zoom In Zoom Out Reset image size

We used KeplerCam, mounted on the 1.2 m telescope located at the Fred Lawrence Whipple Observatory (FLWO) atop Mt. Hopkins, Arizona, to observe a full transit on 2020 February 17. KeplerCam has a 231 × 231 field of view and operates in binned by 2 mode producing a pixel scale of 0672. Images were obtained in the i band with an exposure time of 300 s. A total of 29 images were collected over 144 minutes. The data were reduced using standard IDL routines, and photometry was performed using the AstroImageJsoftware package (Collins et al. 2017). These data are shown in Figure 9.

We observed a full transit on 2020 November 13, using Sinistro, an optical camera mounted on a 1 m telescope located at McDonald Observatory in Texas, operated by Las Cumbres Observatory (LCO; Brown et al. 2013). Sinistro has a 265 × 265 field of view with a pixel scale of 0389. We observed 62 images in total during 339 minutes, using a V-band filter, with an exposure time of 5 minutes. The data were reduced by the standard BANZAI pipeline by the Las Cumbres Observatory Global Telescope (LCOGT; McCully et al. 2018), and photometry was performed using AstroImageJsoftware. These data are shown in Figure 9.

Multicolor Simultaneous Camera for studying Atmospheres of Transiting exoplanets 3 (MuSCAT3) is a multiband simultaneous camera installed on the 2 m Faulkes Telescope North at LCO on Haleakala, Maui (Narita et al. 2020). It has four channels, enabling simultaneous photometry in the g (400–550 nm), r (550–700 nm), i (700–820 nm), and z_s (820–920 nm) bands. Each channel has a 2048 × 2048 pixel CCD camera with a pixel scale of 027, providing a 91 × 91 field of view. We observed a full transit of TOI-1696.01 on 2020 December 23, from BJD 2459206.703523 to 2459206.827246. We took 36, 41, 89, and 131 exposures with exposure times of 300, 265, 120, and 80 s in the g, r, i, and z_s bands, respectively.

The data reduction was conducted by the standard LCOGT BANZAI pipeline. Then differential photometry was conducted by a customized aperture-photometry pipeline for the MuSCAT series (Fukui et al. 2011). The optimized aperture radii are 8, 6, 10, and 8 pixels (216, 162, 27, and 216) for the g, r, i, and z_s bands, respectively. We optimized a set of comparison stars for each band to minimize the dispersion of the light curves. For computational efficiency, and to achieve a more uniform S/N, we subsequently binned the g, r, i, and z_s data to 300, 240, 180, and 120 s, respectively. These data are shown in Figure 7.

We also observed a full transit with MuSCAT (Narita et al. 2015), which is installed on the 188 cm telescope of the National Astronomical Observatory of Japan (NAOJ) in Okayama, Japan. MuSCAT has a similar optical design as MuSCAT3 but has three CCD cameras for the g, r, and z_s bands. On the night of 2021 July 28 we observed TOI-1696 from BJD 2459424.228358 to 2459424.30679. At that point, the r-band camera was not available due to an instrumental issue, so we observed with only the g and z_s bands, using an exposure time of 60 s for both bands.

The data reduction and differential photometry was performed using the pipeline described in Fukui et al. (2011). The optimized aperture radii were 4 and 6 pixels (144 and 216) for the g and z_s bands, respectively. Similarly to the MuSCAT3 data, we binned the g and z_s data to 300 and 120 s, respectively. These data are shown in Figure 8.

On the nights of 2020 December 03 and 2021 October 14, TOI-1696 was observed with the 'Alopeke speckle imager (Scott 2019), mounted on the 8.1 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 = 5.8 at a separation of 1'' in the 832 nm band. No secondary sources were detected. The data taken on 2021 October 14 are shown in Figure 3.

Zoom In Zoom Out Reset image size

On 2021 September 19 we conducted near-infrared high-resolution imaging using the adaptive optics instrument Palomar High Angular Resolution Observer (PHARO) mounted on the 5 m Hale telescope at Palomar Observatory (Hayward et al. 2001). We observed TOI-1696 separately in the Br γ (2.18 μm) and H_cont (2.29 μm) bands, reaching a contrast of Δmag = 8 at a separation of 1'' in both bands. The AO images and corresponding contrast curves are shown in Figure 4.

Zoom In Zoom Out Reset image size

We obtained high-resolution spectra of TOI-1696 in the near-infrared with the InfraRed Doppler (IRD; Tamura et al. 2012; Kotani et al. 2018), mounted on the 8.2 m Subaru telescope. IRD can achieve a spectral resolution of ∼70,000 in the wavelength range 930 nm to 1740 nm. The derived spectra were used for three purposes: to search for spectral companions (e.g., double-lined spectroscopic binary scenarios), to measure fundamental stellar parameters (e.g., effective temperature and metallicity), and to rule out large radial velocity (RV) variations that would indicate an eclipsing binary (EB), as well as placing a limit on the mass of the planet. From UT 2021 January 30 to 2022 January 08, we obtained 13 spectra of TOI-1696 using 1800 s exposure times, as part of a Subaru Intensive Program (Proposal IDs S20B-088I and S21B-118I). The raw data were reduced using IRAF (Tody 1993) as well as a pipeline for the detector's bias processing and wavelength calibrations developed by the IRD instrument team (Kuzuhara et al. 2018; Hirano et al. 2020). For the RV analyses and stellar parameter derivation, we computed a high-S/N coadded spectrum of the target following the procedures described in Hirano et al. (2020).

For use as a spectral template in the analysis described in Section 3.2, we also downloaded archival IRD data of GJ 699 (Barnard's Star), ⁴⁰ which were obtained on 2019 March 23 (HST). We reduced and calibrated the GJ 699 data following the same procedures as the TOI-1696 data.

We collected observations of TOI-1696 on UT 2020 December 09 using SpeX, a medium-resolution spectrograph on the NASA Infrared Telescope Facility (IRTF) on Maunakea (Rayner et al. 2003). We obtained our observations in SXD mode with a 03 × 15'' slit, providing a spectral resolution of R ≈ 2000 over a wavelength range from 700 nm to 2550 nm. In order to remove the sky background and reduce systematics, the spectra were collected using an ABBA nod pattern (with a separation of 75 between the A and B positions) and with the slit synced to the parallactic angle. We reduced our spectra using the Spextool reduction pipeline (Cushing et al. 2004) and removed telluric contamination using xtellcor (Vacca et al. 2003). The derived spectra were used to calculate the stellar metallicity.

In the next subsections we estimate the fundamental stellar parameters of TOI-1696. First, the stellar effective temperature T_eff and metallicity [Fe/H] are derived from two independent methods; one is from the IRD spectra, and the other is from the SpeX spectra and photometric relations. Second, the stellar radius R_⋆, mass M_⋆, and other related parameters are derived using empirical relations and the above T_eff and [Fe/H] values.

3.1.1. Estimation of T_eff and [Fe/H]: from IRD Spectra

3.1.2. Estimation of T_eff and [Fe/H]: from SpeX Spectra and Photometric Relations

3.1.3. Estimation of Stellar Radius and Mass

If a stellar companion orbits the target star, the observed spectra will generally be the combination of two stellar spectra with different radial velocities. To see if TOI-1696 is a double-lined spectroscopic binary, we calculated the cross-correlation function (CCF) of the TOI-1696's IRD spectra with that of the well-known single-star GJ 699 (Barnard's Star). The spectrum of TOI-1696 used for the analysis was obtained on UT 2021 January 30 08:53, which corresponds to an orbital phase of 0.247 based on the TESS ephemeris.

For the analysis, we divided the spectra into six wavelength bins that are less affected by telluric absorption: [988, 993 nm], [995, 1000 nm], [1009, 1014 nm], [1016, 1021 nm], [1023, 1028 nm], and [1030, 1033 nm]. We corrected the telluric absorption signal using the spectra of the rapid-rotator HIP 74625, which was observed on the same night. The CCF to the template spectrum was calculated for each segment, after barycentric velocity correction. Finally, we computed the median of the CCFs from each segment. As shown in Figure 5, the resulting CCF is clearly single-peaked. If the observed transit signals were actually caused by an eclipsing stellar companion, the RV difference at quadrature would be >100 km s⁻¹, which would result in a second peak in the CCF given that the flux of such a companion would be detectable. We thus conclude TOI-1696 is not an eclipsing binary.

Zoom In Zoom Out Reset image size

Because young stars are active and rapidly rotating, the stellar activity and rotation period can be used as proxy for determining its youth. We did not find any stellar rotational signal in the TESS SPOC light curve, suggesting that the star is not very active. Similarly, no strong rotational signal was found in archival photometric data from Zwicky Transient Facility (ZTF) Data Release 9 (Bellm et al. 2019; Masci et al. 2019) and the All Sky Automated Survey for SuperNovae (ASAS-SN; Kochanek et al. 2017).

GJ 699 has a rotation period of 145 days and $v\sin i$ of less than 3 km s⁻¹ (Toledo-Padrón et al. 2019), which is below the limit of IRD's resolving power (∼70000, corresponding to ∼4.5 km s⁻¹). While the CCF of TOI-1696 has an FWHM value consistent with that of GJ 699 (see Figure 5), even if we assume the rotation axis of TOI-1696 is in the plane of the sky, relatively short rotation periods cannot be ruled out, as their rotational broadening would not be resolvable with IRD. However, fast rotation would most likely be accompanied with surface magnetic activity levels that would produce detectable photometric signals. We also used banyan Σ (Gagné et al. 2018) to check if TOI-1696 is a member of any known stellar associations, using its proper motion and the parallax from Gaia EDR3. banyan Σ tool ⁴³ returned a value of 99.9% field star, suggesting it is not a member of any nearby young moving group. The nondetection by Galaxy Evolution Explorer also means that the star is not young enough to be bright in the UV. We thus conclude that TOI-1696 is most likely a relatively old, slow rotator.

We jointly fit the TESS, KeplerCam, Sinistro, MuSCAT3, and MuSCAT data sets using the PyMC3 ⁴⁴ (Salvatier et al. 2016), exoplanet ⁴⁵ (Foreman-Mackey et al. 2019), starry (Luger et al. 2019), and celerite2 (Foreman-Mackey et al. 2017; Foreman-Mackey 2018) software packages.

3.4.1. Mean Model

3.4.2. Noise Model

3.4.3. Priors

We placed Gaussian priors on the white noise scale parameters, with the center and width equal to unity, and Gaussian priors on the stellar mass and radius based on the results in Table 1. We also placed Gaussian priors on the limb darkening coefficients based on interpolation of the parameters tabulated in Claret et al. (2012) and Claret (2017), propagating the uncertainties in the stellar parameters in Table 1 via Monte Carlo simulation. ⁴⁶ We used uniform priors for all other free parameters, except for the GP hyperparameters (red noise amplitude and length scale), which we sampled in logarithmic space with wide (uninformative) Gaussian priors.

Zoom In Zoom Out Reset image size

Table 1. Main Identifiers, Equatorial Coordinates, Proper Motion, Parallax, Optical and Infrared Magnitudes, and Fundamental Parameters of TOI-1696

Parameter	Value	Source
Main identifiers
TIC	470381900	TIC v8 a
2MASS	J04210733+4849116	ExoFOP a
WISE	J042107.34+484911.5	ExoFOP a
UCAC4	695-028795	ExoFOP a
Gaia EDR3	270260649602149760	Gaia EDR3 b
Equatorial coordinates, parallax, and proper motion
R.A. (J2015.5)	04h21m0736	Gaia EDR3 b
Decl. (J2015.5)	+48 ° 49'1138	Gaia EDR3 b
π (mas)	15.4752 ± 0.0345	Gaia EDR3 b
μα (mas yr−1)	12.8726 ± 0.0345	Gaia EDR3 b
μδ (mas yr−1)	−19.0463 ± 0.0269	Gaia EDR3 b
Optical and near-infrared photometry
TESS	13.9664 ± 0.00730068	TIC v8 a
G	15.3056 ± 0.0028	Gaia EDR3 b
Bp	17.0511 ± 0.0051	Gaia EDR3 b
Rp	14.0457 ± 0.0039	Gaia EDR3 b
B	18.467 ± 0.162	ExoFOP a
V	16.82 ± 1.133	ExoFOP a
J	12.233 ± 0.023	2MASS c
H	11.604 ± 0.031	2MASS c
Ks	11.331 ± 0.023	2MASS c
W1	11.134 ± 0.023	AllWISE c
W2	10.984 ± 0.021	AllWISE c
W3	10.71 ± 0.11	AllWISE d
W4	8.748±	AllWISE d
Fundamental parameters
Teff(K)	3185 ± 76	This work
$\mathrm{log}g$(cgs)	4.959 ± 0.026	This work
[Fe/H](dex)	0.336 ± 0.060	This work
M⋆(M⊙)	0.255 ± 0.0066	This work
R⋆(R⊙)	0.2775 ± 0.0080	This work
ρ⋆(g cm−3)	${16.8}_{-1.4}^{+1.5}$	This work
distance (pc)	65.03 ± 0.36	This work
Luminosity (L⊙)	${0.00711}_{-0.00075}^{+0.00083}$	This work

Notes.

^a ExoFOP (2019). ^b Stassun & Torres (2021). ^c Skrutskie et al. (2006). ^d Cutri et al. (2021).

Download table as:  ASCIITypeset image

3.4.4. Exploration of Parameter Space

To optimize the model we used the gradient-based BFGS algorithm (Nocedal & Wright 2006) implemented in scipy.optimize to find the initial maximum a posteriori (MAP) parameter estimates. We then used these estimates to initialize an exploration of the parameter space via "no U-turn sampling" (NUTS; Hoffman & Gelman 2014), an efficient gradient-based Hamiltonian Monte Carlo (HMC) sampler implemented in PyMC3. We sampled four independent chains initialized with the MAP parameter values for 3000 tuning steps followed by an additional 1000 steps. The MAP solution found during sampling had a higher posterior probability than that found by BFGS, so we repeated the sampling procedure once more starting from the new MAP values. The Gelman–Rubin $\hat{R}$ statistic (Gelman & Rubin 1992) was less than 1.003 for all parameters, indicating that the samplers had converged, and the number of effectively independent samples was large enough to ensure negligibly small sampling errors. The results are summarized in Table 2, and detailed plots showing the model fits to the TESS and ground-based data sets are shown in Figures 6, 7, 8, and 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Results of Joint Fit to the TESS and Ground-based Transit Data Sets

Parameter	Value
Primary transit parameters
M⋆ [M⊙]	0.255 ± 0.007
R⋆ [R⊙]	0.277 ± 0.008
T0 [BJD]	2458834.20115 ± 0.00058
P [days]	2.500311 ± 0.000004
b	${0.59}_{-0.04}^{+0.03}$
RP /R⋆ (T)	0.0952 ± 0.0062
RP /R⋆ (V)	0.1021 ± 0.0057
RP /R⋆ (g)	${0.1036}_{-0.0068}^{+0.0060}$
RP /R⋆ (r)	0.1053 ± 0.0034
RP /R⋆ (i)	0.1023 ± 0.0020
RP /R⋆ (z)	0.1026 ± 0.0020
RP /R⋆	0.1025 ± 0.0014 a
Limb darkening parameters
u1 (T)	0.16 ± 0.01
u2 (T)	0.48 ± 0.01
u1 (V)	0.48 ± 0.02
u2 (V)	0.30 ± 0.01
u1 (g)	0.49 ± 0.01
u2 (g)	0.31 ± 0.01
u1 (r)	0.50 ± 0.01
u2 (r)	0.25 ± 0.01
u1 (i)	0.37 ± 0.01
u2 (i)	0.28 ± 0.01
u1 (z)	0.24 ± 0.01
u2 (z)	0.36 ± 0.01
Derived parameters
Rp [R⊕]	3.09 ± 0.11 a
a [au]	0.0229 ± 0.0002
Teq [K]	489 ± 13 b
T14 [hours]	1.00 ± 0.01

Notes.

^a Derived from achromatic transit model fit. ^b Assuming a Bond albedo of 0.3.

Download table as:  ASCIITypeset image

3.4.5. Achromatic Model

We did not detect any significant wavelength dependence of the transit depth (see Figure 10), which rules out many plausible false-positive scenarios involving eclipsing binaries (see Section 3.6 for more details). Having established that the transit depth is achromatic, we robustly determined the planet radius by conducting a second fit with a single R_P /R_⋆ parameter for all data sets. This fit resulted in a final value of R_P /R_⋆ = 0.1025 ± 0.0014, corresponding to an absolute radius of 3.09 ± 0.11 R_⊕, and all other parameters were unchanged.

Zoom In Zoom Out Reset image size

To put a limit on the mass of TOI-1696.01, we fit an RV model with a circular orbit to the RV data from Subaru/IRD. Between the H-band and the YJ-band spectra obtained with IRD, we opted to use the H-band spectra for RV analysis because of its higher S/N. ⁴⁷ The data observed on 2021 January 29 were excluded because of the possibility of an RV offset, as there was a gap of 8 months relative to the succeeding observations. We also removed any data with the clouds passing, which can cause systematic errors. The final data set consisted of 9 RV measurements from 2021 September 29 to 2022 January 8.

We used the RV model included in PyTransit(Parviainen 2015), which we simplified to have five free parameters: phase-zero epoch T₀, period, RV semi-amplitude, RV zero-point, and RV jitter term. Here, we defined the word RV jitter term' as a parameter that describes the amplitude of white-noise-like scatter due to the systematics. For T₀ and the period, we put Gaussian priors using T₀ and the period derived from the transit analysis. For the other parameters we put wide uniform priors. We ran the built-in Differential Evolution optimizer and then sampled the parameters with Markov Chain Monte Carlo (MCMC) using 30 walkers and 10⁴ steps. We use the following equation to derive the planet mass,

$$\begin{eqnarray}&&{M}_{p}={\left(\displaystyle \frac{{{PM}}_{\star }^{2}}{2\pi }\right)}^{1/3}\displaystyle \frac{K{\left(1-{e}^{2}\right)}^{1/2}}{\sin (i)},\end{eqnarray} \tag{ 1 }$$

where M pis the planet mass, M_⋆ is the star mass, P is the orbital period, K is the RV semi-amplitude, e is the eccentricity (fixed to zero), and i is the inclination (fixed to 90°). To propagate uncertainties, we use the posteriors for M_⋆ and P from previous analyses.

In Figure 11 we plot Keplerian orbital models corresponding to different masses encompassing the 68th, 95th, and 99.7th percentiles of the semi-amplitude posterior distribution. The 2σ upper limit is 48.8 M_⊕, which places the companion 2 orders of magnitude below the deuterium burning mass limit. The best-fit semi-amplitude is K = 14.4 ms⁻¹, which corresponds to a mass of M_p = 12.3M_⊕, and the best-fit jitter value is σ_K = 62 m s⁻¹.

Zoom In Zoom Out Reset image size

We calculated an expected planetary mass of ∼8 M_⊕ with MRExo, ⁴⁸ which uses a mass–radius relationship calibrated for planets around M dwarfs (Kanodia et al. 2019). This mass corresponds to a semi-amplitude of 9.4 m s⁻¹, but the observed RV data exhibit significantly larger variability (σ ≈ 52 m s⁻¹). We interpret this variability as being responsible for the large jitter value found by the fit, which suggests it is out-of-phase with TOI-1696.01. As the star appears to be quiet, one explanation for this signal is the existence of an additional (possibly nontransiting) planet, but more RV measurements would be required to determine if this is the case. Furthermore, if such a planet were dynamically interacting with TOI-1696.01, then this could help explain TOI-1696.01's location in a sparsely populated part of the period–radius plane (see Section 4).

A number of astrophysical scenarios can mimic the transit signal detected from TESS photometry, including an eclipsing binary (EB) with a grazing transit geometry, a hierarchical EB (HEB), and a diluted eclipse of a background (or foreground) EB (BEB) along the line of sight of the target. In the following, we will examine the plausibility of each scenario.

First, the renormalized unit weight error (RUWE) from Gaia EDR3 is 1.12, which suggests that TOI-1696 is single (Belokurov et al. 2020). We can also rule out the EB scenario based on the analysis of the IRD CCF in Section 3.2, and the mass constraint derived in Section 3.1.1. Finally, the absence of any wavelength dependence of the transit depth from our chromatic transit analysis (Section 3.4) is incompatible with contamination from a star of different spectral type (color) than the host star, the details of which are discussed in Appendix B. In the absence of dilution, the measured radius of 3.09 ± 0.11 R_⊕ (0.27 R_Jup) equals the true radius, which makes it significantly smaller than the lower limit of 0.8 R_Jup expected for brown dwarfs (Burrows et al. 2011).

Grazing transit geometries can also be eliminated, as the impact parameter is constrained to b < 0.7 at the 99% level based on our transit and contamination analyses. The apparent boxy shape of our follow-up light curves is in stark contrast with the V-shaped transit expected for grazing orbits. Hence, grazing EB scenario is ruled out.

Moreover, we can constrain the classes of HEBs that can reproduce the observed transit depth and shape using our multiband observations. We aim to compute the eclipse depths for a range of plausible HEBs in the bluest and reddest bandpasses where they are expected to vary significantly. We adopt the method presented in Bouma et al. (2020) to perform the calculation taking into account nonzero impact parameter, the details of which are discussed in Appendix C. Comparing the simulated eclipse depths with the observed depth in each band, we found that there is no plausible HEB configuration explored in our simulation that can reproduce the observed depths in multiple bands simultaneously. Hence, the HEB scenario is ruled out.

Although, TOI-1696's probability of being a BEB is very high a priori given its location at the galactic plane, we argue in the following that the BEB scenario is extremely unlikely. Our MuSCAT3 observation can resolve the signal down to 3'', which represents the maximum radius within which the signal must originate. Furthermore, our high-resolution speckle imaging ruled out any nearby star and blended sources down to 01 at a delta mag of 4.5. We checked archival images taken more than 60 yr apart, but the proper motion of TOI-1696 is not enough to obtain a clear view along the line of sight of the star. However, we can use statistical arguments to estimate the probability of a chance-aligned star. To do this, we use the population synthesis code Trilegal ⁴⁹ (Girardi et al. 2005), which can simulate the Galactic stellar population along any line of sight. Given the position of TOI-1696, we found a probability of 5 × 10⁻⁸ to find a star brighter than T = 16, ⁵⁰ within an area equal to the smallest MuSCAT3 photometric aperture (aperture radius = 3''). Assuming all such stars are binary and preferentially oriented edge-on to produce eclipses with period and depth consistent with the TESS detection, then this can represent a very conservative upper limit of a BEB scenario. Despite the small probability of a BEB based on the trivial star counting argument, we discuss relevant tools in the following section for a more thorough statistical modeling.

Here we quantify the false-positive probability (FPP) of TOI-1696.01 using the Python package Vespa and Triceratops (Morton 2015a; Giacalone & Dressing 2021), the details of which are discussed in Appendix D. Although we were able to rule out the classes of EB, BEB, and HEB in Section 3.6, we ran Vespa considering all these scenarios for completeness and computed a formal FPP < 1 × 10⁻⁶, which robustly quantifies TOI-1696.01 as a statistically validated planet. Additionally, we validated TOI-1696.01 using Triceratopsand found FPP = 2 × 10⁻³. Giacalone et al. (2021) noted that TOIs with FPP < 0.015 have a high enough probability of being bona fide planets to be considered validated. The low FPPs calculated using Vespa and Triceratops added further evidence to the planetary nature of TOI-1696.01. We now refer to the planet as TOI-1696 b in the remaining sections.

Here, we consider the nature of TOI-1696 b by placing it in context with the population of known exoplanets ⁵¹ (NASA Exoplanet Science Institute 2020). Figure 12 shows a radius versus period diagram, indicating that there are only a handful of planets with similar characteristics to TOI-1696 b. The measured planetary radius R_p of 3.09 ± 0.11 R_⊕ and the orbital period P of 2.50031 ± 0.00001 days, places it securely within the bounds of the Neptunian desert as defined by Mazeh et al. (2016). The region occupied by TOI-1696 remains sparsely populated despite recent discoveries of TESS planets within the Neptunian desert (e.g., Murgas et al. 2021; Brande et al. 2022).

Zoom In Zoom Out Reset image size

It should be noted that the Neptunian desert was originally determined based on a population of planets orbiting mainly solar-type stars from the Kepler mission. Because TOI-1696 is an M dwarf, the incident flux at a given orbital separation will be less than for solar-type stars. Nevertheless, we emphasize that the target exists in a sparsely populated region of the parameter space, despite the large number of planets discovered around M dwarfs since the Kepler mission (i.e., from K2 and TESS). For example, if we limit the comparison to the 279 confirmed planets around M dwarfs with T_eff below 3800 K, only 14 planets have been found so far with orbital periods shorter than 10 days and planetary radii in the range of 2.5R_⊕ < R_p < 5 R_⊕. As shown in Figure 12, TOI-1696 b is similar to K2-25 b (Mann et al. 2016), K2-320 b (Castro González et al. 2020), GJ 1214 b (Charbonneau et al. 2009), TOI-269 b (Cointepas et al. 2021), and TOI-2406 b (Wells et al. 2021) in terms of the orbital period and radius. In particular, TOI-2406b appears most similar to TOI-1696 b as it orbits around a mid-M dwarf with an effective temperature of 3100 ± 75, and has a radius of 2.94 ± 0.17 R_⊕ and orbital period of 3.077 days. TOI-2406 is also thought to be relatively old without any activity signal. As both TOI-1696 b and TOI-2406b are excellent targets for detailed characterization studies, together they may provide unique insights into this class of planet. There is also some similarity between TOI-1696 b and the Neptunian Desert planets orbiting young host stars, such as au Mic b and c, K2-25 b, K2-95 b, and K2-264 b. It has been suggested that these planets may have inflated radii and could possibly still be undergoing atmospheric mass loss (e.g., Mann et al. 2016). Further study of TOI-1696 b could reveal whether its similarity to these planets (despite being older) is only superficial, or if it is indicative of an inflated radius.

Given the rarity of this planet, it would be useful to assess its prospects for future atmospheric observations to understand its formation and evolution. In particular, the relatively large size of the planet compared to its host star makes it a good candidate for transmission spectroscopy. Using Equation (1) in Kempton et al. (2018), we calculated the transmission spectroscopy metric (TSM) of TOI-1696 b from its mass, radius, equilibrium temperature, stellar radius, and J-band magnitude. We used the values in Tables 1 and 2, and assumed a mass of 8 M_⊕ estimated by MRExo. The derived TSM value of TOI-1696 b is 105.6. For reference, Kempton et al. (2018) suggested that planets with TSM > 90 are ideal targets for atmospheric follow-up.

For comparison, we calculated the TSM for the known population of transiting M dwarf planets. We selected planets with T_eff < 3800 K, R_p < 10 R_⊕, and H < 11mag. ⁵² For planets without mass measurements, we assumed the masses predicted by MRExo. For planets without an equilibrium temperature, we estimated it from the semimajor axis and the host star's effective temperature (assuming zero albedo). Figure 13 shows the computed TSM values for the selected samples of planets. The TSM of TOI-1696 b places it in the top 10, making it one of the best targets for future atmospheric investigations.

Zoom In Zoom Out Reset image size

Up to this point in the section, the discussion has been based on the assumption that the target has an atmosphere. Usually it is thought that planets above the so-called radius gap can retain their atmospheres (Weiss & Marcy 2014; Rogers 2015). However, does TOI-1696 b actually have a H₂/He atmosphere? Here we study the atmospheric mass that TOI-1696 b can retain after ∼8 Gyr under a stellar XUV irradiation. The mass of TOI-1696 b remains poorly constrained as discussed in Section 3.5. We modeled TOI-1696 b as a rocky planet with Earth-like core compositions (MgSiO₃:Fe = 7:3) in the mass range from 0.5 M_⊕ to 20 M_⊕. The silicate mantle and iron core were described by the third-order Birch–Murnagham equation of state (EoS) for MgSiO₃ perovskite (Karki et al. 2000; Seager et al. 2007) and the Vinet EoS for -Fe (Anderson et al. 2001), respectively. The Thomas–Fermi–Dirac EoS (Salpeter & Zapolsky 1967) was applied to high-pressure EoS for MgSiO₃ at P ≥ 4.90 TPa and Fe at P ≥ 2.09 × 10⁴ GPa (Seager et al. 2007; Zeng & Sasselov 2013). The pressure and temperature in a H₂/He envelope were calculated using the SCvH EoS (Saumon et al. 1995).

We computed the thermal evolution of TOI-1696 b with a H₂/He atmosphere by calculating its interior structure in hydrostatic equilibrium for ∼8 Gyr and calculated its mass-loss process. The initial mass fraction of a H₂/He atmosphere for a rocky planet ranges from 0.001% to 30% of its core mass. The energy-limited hydrodynamic escape (Watson et al. 1981) controls the mass-loss rate given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{{\rm{p}}}}{{dt}}=-\eta \displaystyle \frac{{R}_{{\rm{p}}}^{3}{L}_{\mathrm{XUV}}(t)}{4{{GM}}_{{\rm{p}}}{a}^{2}{K}_{\mathrm{tide}}},\end{eqnarray} \tag{ 2 }$$

where η is the heating efficiency due to stellar XUV irradiation, L_XUV is the stellar XUV luminosity, G is the gravitational constant, and R_p is the planetary radius (Erkaev et al. 2007). As the heating efficiency for a hydrogen-rich upper atmosphere was lower than 20% (Shematovich et al. 2014; Ionov & Shematovich 2015), we adopted η = 0.1. K_tide is the reduction factor of a gravitational potential owing to the effect of a stellar tide:

$$\begin{eqnarray}&&{K}_{\mathrm{tide}}(\xi )=1-\displaystyle \frac{3}{2\xi }+\displaystyle \frac{1}{2{\xi }^{3}}\lt 1,\,\,\xi =\displaystyle \frac{{R}_{{\rm{H}}}}{{R}_{{\rm{p}}}},\end{eqnarray} \tag{ 3 }$$

where R_H is the Hill radius. The XUV luminosity (L_XUV) of TOI-1696 followed from the X-ray-to-bolometric luminosity relations of M-type stars (Jackson et al. 2012), where we adopted the current luminosity of TOI-1696 as its bolometric luminosity. We also considered a 10L_XUV model because of the large uncertainty in L_XUV of young M dwarfs.

Figure 14 shows the initial H₂/He atmosphere of a TOI-1696 b-like planet that reproduces the radius of 3.09 ± 0.11R_⊕ at the current location (i.e., T_eq = 489 ± 13 K) after the mass loss driven by the standard XUV radiation (L_XUV: blue) and 10 times higher one (10L_XUV: red). The gray region shows the H₂/He atmospheric mass fraction of TOI-1696 b with a rocky core that satisfies its observed radius. The observed radius of TOI-1696 b favors the existence of a H₂/He atmosphere atmosphere with ≳3 wt% (a percentage by weight) unless its core contains icy material. We find that TOI-1696 b can possess the H₂/He atmosphere for 8 Gyr if its core mass is larger than ∼1.5 M_⊕ (∼4 M_⊕ for 10L_XUV models). We rule out the possibility that TOI-1696 b with a rocky core initially had a H₂/He atmosphere of ≲3%. Also, TOI-1696 b with mass of ≳10 M_⊕ can retain almost all the H₂/He atmosphere accreted from a disk. These suggest that if TOI-1696 b has a core of ≳1.5–4 M_⊕, it is likely to be a sub-Neptune with a H₂/He atmosphere.

Zoom In Zoom Out Reset image size

We calculate the abundance of seven elements in addition to iron from IRD spectra. We used 28 lines in total caused by neutral atoms of Na, Mg, Ca, Ti, Cr, Mn, and Fe and singly ionized Sr. The detailed procedures of abundance analysis and error estimation are described in Ishikawa et al. (2020). Figure A1 shows the final values of abundance after the iteration. From the final values of the abundances of the eight elements, [M/H] was determined by calculating the average weighted by the inverse of the square of their estimated errors.

Zoom In Zoom Out Reset image size

As an independent determination of the basic stellar parameters, we performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia EDR3 parallax (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, 2018). We pulled the JHK_S magnitudes from 2MASS, the W1–W3 magnitudes from the Wide-field Infrared Survey Explorer (WISE), and the grizy magnitudes from the Panoramic Survey Telescope and Rapid Response System. Together, the available photometry spans the full stellar SED over the wavelength range 0.4–10 μm (see Figure A2).

Zoom In Zoom Out Reset image size

We performed a fit using NextGen stellar atmosphere models (Hauschildt et al. 1999), with the effective temperature (T_eff) and metallicity ([Fe/H]) constrained from the spectroscopic analysis. The remaining free parameter is the extinction A_V , which we fixed at zero due to the star's proximity. The resulting fit (Figure A2) has a reduced χ² of 1.7. Integrating the (unreddened) model SED gives the bolometric flux at Earth, F_bol = 5.20 ± 0.25 × 10⁻¹¹ erg s⁻¹ cm⁻². Taking the F_bol and T_eff together with the Gaia parallax, gives the stellar radius, R_⋆ = 0.276 ± 0.015 R_⊙. We used the T_eff and [Fe/H] values from spectroscopic results as priors for the parameter estimation.

In addition, we estimated the stellar mass from the empirical relations of Mann et al. (2019), giving M_⋆ = 0.279 ± 0.014 M_⊙. Finally, the radius and mass together imply a mean stellar density of ρ_⋆ = 18.79 ± 3.26 g cm⁻³.

In addition to the methods described above, we used the Python package isochrones, which calculates stellar parameters from the stellar evolution models. The three methods are not fully independent, as some of them use the same relations such as mass derivation from Mann et al. (2019), but comparing three results are useful to confirm the results are robust. The derived stellar parameters agreed within 1 ∼ 2σ, as shown in Table A1. We pick up the results from the empirical relations as our final stellar parameters in Table 1.

Table A1. Stellar Parameters Derived from Empirical Relations (Method 1; Section 3.1.3), SED Fitting (Method 2; Appendix A.2), and isochrones (Method 3; Appendix A.3)

Parameter	Method 1	Method 2	Method 3
Teff (K)	⋯	3130 ± 75	3159 ± 40
[Fe/H] (dex)	⋯	0.2 ± 0.3	${0.232}_{-0.038}^{+0.035}$
M⋆ (M⊙)	0.255 ± 0.0066	0.279 ± 0.014	${0.276}_{-0.005}^{+0.006}$
R⋆ (R⊙)	0.2775 ± 0.0080	0.280 ± 0.014	0.291 ± 0.005
$\mathrm{log}g$ (cgs)	4.959 ± 0.026	4.990 ± 0.049	${4.955}_{-0.007}^{+0.008}$
ρ⋆ (g cm−3)	${16.8}_{-1.4}^{+1.5}$	18.0 ± 2.9	${16.282}_{-0.617}^{+0.598}$
distance (pc)	65.03 ± 0.36	⋯	${65.390}_{-0.464}^{+0.540}$
Luminosity (L⊙)	${0.00711}_{-0.00075}^{+0.00083}$	⋯	⋯
AV (mag)	0 (fixed)	0 (fixed)	∼0 †
Fbol (cgs)	⋯	5.19(18) × 10−11	⋯
MKs (mag)	7.265 ± 0.026	⋯	⋯

Note.^† The resulting posterior is approximately zero and non-Gaussian as a result of using a tight uniform prior close to 0 for numerical reasons.

Download table as:  ASCIITypeset image