The Discovery and Follow-up of Four Transiting Short-period Sub-Neptunes Orbiting M Dwarfs

Sub-Neptunes with radii of 2–3 R ⊕ are intermediate in size between rocky planets and Neptune-sized planets. The orbital properties and bulk compositions of transiting sub-Neptunes provide clues to the formation and evolution of close-in small planets. In this paper, we present the discovery and follow-up of four sub-Neptunes orbiting M dwarfs (TOI-782, TOI-1448, TOI-2120, and TOI-2406), three of which were newly validated by ground-based follow-up observations and statistical analyses. TOI-782 b, TOI-1448 b, TOI-2120 b, and TOI-2406 b have radii of Rp=2.740−0.079+0.082R⊕ , 2.769−0.068+0.073R⊕ , 2.120 ± 0.067 R ⊕, and 2.830−0.066+0.068R⊕ and orbital periods of P = 8.02, 8.11, 5.80, and 3.08 days, respectively. Doppler monitoring with the Subaru/InfraRed Doppler instrument led to 2σ upper limits on the masses of <19.1 M ⊕, <19.5 M ⊕, <6.8 M ⊕, and <15.6 M ⊕ for TOI-782 b, TOI-1448 b, TOI-2120 b, and TOI-2406 b, respectively. The mass–radius relationship of these four sub-Neptunes testifies to the existence of volatile material in their interiors. These four sub-Neptunes, which are located above the so-called “radius valley,” are likely to retain a significant atmosphere and/or an icy mantle on the core, such as a water world. We find that at least three of the four sub-Neptunes (TOI-782 b, TOI-2120 b, and TOI-2406 b), orbiting M dwarfs older than 1 Gyr, are likely to have eccentricities of e ∼ 0.2–0.3. The fact that tidal circularization of their orbits is not achieved over 1 Gyr suggests inefficient tidal dissipation in their interiors.


INTRODUCTION
About 70-80% of stars in the solar neighborhood are M dwarfs (e.g.Winters et al. 2019Winters et al. , 2021)).Nearby M dwarfs are suitable targets to search for low-mass planets.Recent discoveries of temperate planets orbiting nearby M dwarfs, such as TRAPPIST-1 (Gillon et al. 2017), provide a unique opportunity for atmospheric characterization of terrestrial planets beyond the solar system.Near-infrared spectrographs designed for highprecision radial velocity (RV) measurements, such as the Very Large Telescope/CRIRES+ (Dorn et al. 2016), the Hobby-Eberly Telescope/Habitable Planet Finder (Mahadevan et al. 2018), Subaru/InfraRed Doppler (Kotani et al. 2018), CARMENES (Quirrenbach et al. 2014), PARVI (Gibson et al. 2020), NIRPS (Wildi et al. 2017), CFHT/SPIRou (Donati et al. 2020), and Gemini/MAROON-X (Seifahrt et al. 2018), have begun to find planets down to an Earth mass around M dwarfs, such as Teegarden's Star (Zechmeister et al. 2019).Also, space-based planet surveys have accelerated the detection of planets orbiting M dwarfs.The successor to Kepler, the Transiting Exoplanets Survey Satellite (TESS), which launched in 2018, has conducted the first-ever allsky photometric survey for planets around bright stars (Ricker et al. 2015).The TESS bandpass spanning from 600 to 1000 nm is well suited for detecting planets around cool stars.About one-fourth of confirmed TESS planets have been discovered around cool stars (Guerrero et al. 2021).
The measured densities of almost all of sub-Neptunes are too low to be consistent with pure rock and metal, and thereby suggest a large mass fraction of volatile material.Many sub-Neptunes with densities lower than those of planets composed of pure H 2 O retain atmospheres up to the present day.The atmospheres of planets in the proximity of a central star reflect the competing processes of the accumulation of the disk gas (e.g.Ikoma & Hori 2012;Lee et al. 2014;Ormel et al. 2015;Kurokawa & Tanigawa 2018;Kuwahara et al. 2019) and secondary gases, as well as atmospheric escape by photoevaporation (e.g.Lopez & Fortney 2013;Owen & Wu 2013, 2017) and core-powered mass loss (Ginzburg et al. 2016(Ginzburg et al. , 2018;;Gupta & Schlichting 2019;Gupta et al. 2022).Understanding the prevalence of sub-Neptunes with significant atmospheres helps to disentangle the accretion histories of low-mass planets.
We report four sub-Neptunes orbiting M dwarfs from TESS (TOI-782 b, TOI-1448 b, TOI-2120 b, and TOI-2406 b) validated by ground-based follow-up observations using multicolor photometry, high-resolution imaging, and Subaru/IRD RV measurements and statistical analyses, three of which are newly confirmed by this study; the discovery of TOI-2406 b was first reported in Wells et al. (2021).We also obtain their (upper limit) masses at the 2σ confidence level from high-precision RV measurements using Subaru/IRD.We find that at least three sub-Neptunes around a mature M dwarf with age ≳ 1 Gyr are likely to have eccentricities of e ∼ 0.2 − 0.3.Such a nonzero eccentricity of a short-period planet invokes inefficient tidal dissipation in the interior (e.g.Wang & Ford 2011).
This paper is structured as follows.We describe the TESS photometry and our ground-based follow-up observations in Section 2. We present data analyses and the properties of four M dwarfs and their planets in Section 3. In Section 4, we discuss the interior structures of short-period sub-Neptunes with nonzero eccentricities.In the last section, we summarize our results.
For further analyses, we downloaded the Presearch Data Conditioning Simple Aperture Photometry (PDC-SAP; Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014) ) where available, otherwise, the TESS-SPOC light curves that are extracted from the FFIs (Caldwell et al. 2020) in the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute.The light curves were nor-malized to unity sector by sector after removing data points with bad-quality flags.

Multicolor Transit Photometry from the Ground
We conducted multiband transit photometry of TOI-782 b, TOI-1448 b, TOI-2120 b, and TOI-2406 b from the ground as part of the TESS Follow-up Observing Program (TFOP; Collins 2019)1 to confirm the transit signals detected by TESS.We checked for chromaticity in the transit depths, that is, a sign of false positives due to contamination from eclipsing binaries and refined the physical parameters of the planetary candidates.We utilized a customized version of the Tapir software package (Jensen 2013), the TESS Transit Finder, to schedule our transit observations.All the observations that are used in our analyses are summarized in Table 1.
Below we briefly describe the observations and data reductions conducted at each facility.

TCS 1.52m/MuSCAT2
We observed one transit of TOI-782.01 and one transit of TOI-1448 b with the multiband imager MuSCAT2 (Narita et al. 2019), which is mounted on the 1.52m TCS telescope at the Teide Observatory in the Canary Islands, Spain.MuSCAT2 has four optical channels each equipped with a 1k × 1k CCD camera with a pixel scale of 0. ′′ 44 pixel −1 , providing a field of view (FOV) of 7. ′ 4 × 7. ′ 4, and is capable of simultaneous imaging in the g, r, i, and z s bands.The exposure times were set independently for each band as shown in Table 1.The airmass range and the median FWHM of the stellar point spread function in each band of each observation are also given in Table 1.We performed image calibration and aperture photometry to extract relative light curves using a pipeline described in Fukui et al. (2011).

FTN 2m/MuSCAT3
We observed one transit of TOI-782 b, one transit of TOI-1448 b, three transits of TOI-2120 b, and one transit of TOI-2406 b with the multiband imager MuSCAT3 (Narita et al. 2020), which is mounted on the 2 m Faulkes Telescope North (FTN), owned and operated by the Las Cumbres Observatory (LCO), located at the Haleakala Observatory in Hawaii, USA.As with MuSCAT2, MuS-CAT3 is capable of simultaneous imaging in the g, r, i, and z s bands, but is equipped with a 2k × 2k CCD camera with a pixel scale of 0. ′′ 27 pixels −1 , providing an FOV of 9. ′ 1 × 9. ′ 1, at each channel.The exposure times, the median FWHM values, and the airmass range for each observation are reported in Table 1.We performed image calibration and aperture photometry on the obtained images in the same way as for the MuSCAT2 data.

LCO 1m/Sinistro
We observed six transits of TOI-782 b, one transit of TOI-1448 b, and one transit of TOI-2120 b using the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network.The telescopes are equipped with 4096 × 4096 Sinistro cameras having an image scale of 0. ′′ 389 per pixel, resulting in an FOV of 26 ′ × 26 ′ .The exposure times, the typical FWHM values, and the airmass range for each observation are given in Table 1.The images were calibrated using the standard LCOGT BANZAI pipeline (McCully et al. 2018) and differential photometric data were extracted using AstroImageJ (Collins et al. 2017).

TRAPPIST-North
We observed one full transit of TOI-2120 b with the 0.6 m TRAPPIST-North telescope located at the Oukaimeden Observatory in Morocco (Jehin et al. 2011;Gillon et al. 2013;Barkaoui et al. 2019) on UT 2021 July 14.TRAPPIST-North is equipped with a thermoelectrically cooled 2K × 2K Andor iKon-L BEX2-DD CCD camera with a pixel scale of 0. ′′ 64, offering an FOV of 20'× 20'.We collected 186 images in the Sloan z ′ band with an exposure time of 70 s.The typical FWHM values and airmass range for this observation are shown in Table 1.We performed data reduction and differential aperture photometry using prose2 (Garcia et al. 2022), which automatically selected the optimum apertures for the photometric data extraction to be 5.05 pixels (3.′′ 23).

High-dispersion Spectroscopy with Subaru/IRD
We obtained time series of high-resolution spectra of TOI-782, TOI-1448, TOI-2120, and TOI-2406 with the IRD instrument on the 8.2 m Subaru telescope (Tamura et al. 2012;Kotani et al. 2018) between UT 2019 June 17 and 2022 October 10, as part of the Subaru/IRD TESS intensive follow-up program (ID: S20B-088I).The IRD is a fiber-fed spectrograph covering near-infrared wavelengths from 930 nm to 1740 nm with a spectral resolution of ≈ 70, 000.The integration time per exposure was set to 600-1200 s depending on the stellar brightness and observing conditions on each night.In total, we obtained 25, 21, 24, and 22 exposures for TOI-782, TOI-1448, TOI-2120, and TOI-2406, respectively.We also observed at least one telluric standard star (A0 or A1 star) on each night to correct the telluric lines when extracting the template spectrum for the RV analysis.
The raw IRD data were reduced by the procedure described in Kuzuhara et al. (2018) and Hirano et al. (2020), where the wavelengths were calibrated by spectra from a laser-frequency comb.The reduced onedimensional spectra have typical signal-to-noise ratios (S/Ns) of 30-60 pixel −1 at 1000 nm, except for TOI-2406, whose S/NR was 13-28 due to its faintness.We discarded data with S/Ns lower than 18, which could be affected by detector persistence.We also discarded the data that are logged as being affected by passing clouds or instrumental troubles (e.g., guiding errors) during the exposures, leaving 24, 16, 17, and 11 spectra for TOI-782, TOI-1448, TOI-2120, and TOI-2406, respectively.The RVs for each target were extracted following the procedure of Hirano et al. (2020) using the forward modeling technique: for each target, the RV-analysis pipeline derives a telluric-free template spectrum of the target, which is used to model and fit the individual observed spectra taking into account the instantaneous variations in the instrumental profiles.The measured RV values and 1σ uncertainties along with the S/N ratios at 1000 nm of the original spectra are reported in Table 2.The median RV internal errors are 5.3, 5.7, 4.6, and 14.5 m s −1 for TOI-782, TOI-1448, TOI-2120, and TOI-2406, respectively.The RV uncertainties vary slightly from frame to frame, reflecting the sky conditions and resulting S/N ratio for each frame.About half of the spectra were obtained under good observing conditions, in which case we achieved the expected S/N ratios and RV internal errors.

High-resolution Imaging
To investigate whether there are companion or background stars close to the host stars, we obtained adap- tive optics (AO) and speckle images of TOI-782, TOI-1448, and TOI-2120 with various facilities as part of TFOP.The observational dates, instruments, filters, and achieved contrasts are summarized in Table 3.The observations for individual targets are briefly described below.

TOI-2120
We conducted speckle imaging observations with Speckle Polarimeter (Safonov et al. 2017) on the SAI-2.5mtelescope on UT 2020 October 27 in the I band, and with 'Alopeke 3 (Scott et al. 2021) on the Gemini-N telescope on UT 2020 December 3 with the 562nm (∆λ = 54nm) and 832nm (∆λ = 40nm) band filters.We achieved a 5σ contrast at 0. ′′ 5 of up to 5.6 mag.We did not detect a nearby star within 2 ′′ in any of the observations.

Spectroscopic Properties
We estimated the metallicity ([Fe/H]) and the effective temperature (T eff ) of four host stars using the IRD spectra following the procedure described in Mori et al. (2022).For TOI-782b, TOI-1448b, and TOI-2120b, we used the template spectra identical to those used for the RV analysis, but for TOI-2406b, we used a combined spectrum before deconvolving the instrument profile to avoid amplifying the additional noise due to the low S/N.The basic concept of our analysis is an equivalent width (EW) comparison of individual atomic and molecular absorption lines between synthetic and observed spectra.To determine the T eff value, 47 FeH molecular lines at 990 − 1012 nm were used (see Ishikawa et al. (2022) for details).We employed a total of 25 atomic lines corresponding to Na I, Mg I, Ca I, Ti I, Cr I, Mn I, Fe I, and Sr II as a measure of the elemental abundances (Ishikawa et al. 2020).We selected only atomic lines with a high S/N.We obtained the EWs of the atomic lines by fitting the synthetic spectra line by line, while the EWs of the FeH lines were obtained from a Gaussian fitting to each line.The logarithm of the surface gravity log g and the microturbulent velocity were fixed at 5.0 and 0.5 km s −1 , respectively, for all the targets.We checked that this assumption causes a negligible change in the final T eff and [Fe/H] estimates.
3 https://www.gemini.edu/sciops/instruments/alopeke-zorro/First, we derived a T eff from the EWs of the FeH lines, assuming solar abundances.Subsequently, we determined the individual abundances of the eight elements, i.e., [X/H], using the T eff value.Then, we carried out these analyses for [X/H] and T eff in an iterative process until the T eff and [Fe/H] matched the values assumed in the previous iteration step within their uncertainties.We obtained [Fe/H] = 0.34 ± 0.15, 0.37 ± 0.14, 0.39±0.22,and −0.26±0.24dex, and T eff = 3379±100 K, 3390 ± 100 K, 3183 ± 100 K, and 3308 ± 100 K for TOI-782, TOI-1448, TOI-2120, and TOI-2406, respectively.Note that our [Fe/H] value for TOI-2406 is consistent with that derived by Wells et al. (2021) within the uncertainties.
Using the IRD spectra, we also derived the systemic RVs of our targets.We fitted Gaussian functions to relatively deep atomic lines (30 − 40 lines) at near-infrared wavelengths and compared the line centers with their vacuum wavelengths.The systemic RVs and their uncertainties were then derived from the mean values and the scatter of the measurements as listed in Table 4. Inputting the systemic RVs together with the coordinates, parallaxes, and proper motions, we computed the Galactic space velocities (U, V, W ) for the four stars using the Python script uvw errs (Rodriguez 2016).The derived U V W velocities are also listed in Table 4.For all of our targets, we found that the sky-projected rotational velocity (v sin i s ) of each star is too small (smaller than the spectral resolution, or ∼4.3 km s −1 ) to detect from IRD spectra in a well-calibrated way.
Next, we updated T eff , R s , and d by fitting a stellaratmospheric model to the SED constructed from the photometric magnitudes.For the photometric magnitudes, we used the catalog magnitudes in the BP , G, RP bands from Gaia DR3, J, H, K s bands from 2MASS, and W1, W2 from AllWISE (Wright et al. 2010;Mainzer et al. 2011) for all host stars except for TOI-2120.For TOI-2120, because there is a nearby background star (TIC 389900766) that was separated from TOI-2120 by ∼9 ′′ , ∼6 ′′ , and ∼4 ′′ at the times when the data were collected by 2MASS, AllWISE, and Gaia, respectively, we did not use the W1 and W2 magnitudes from AllWISE which could be contaminated by the flux from this nearby star.For the fitting, the tabulated BT-Settl model (Allard et al. 2011) (2019). 4As a result, we obtained the values of T eff to be 3370 +29 −32 K, 3412 +28 −32 K, 3131 ± 30 K, and 3100 +31 −24 K, and values of R s of 0.413 ± 0.008 R ⊙ , 0.380 ± 0.007 R ⊙ , 0.245 ± 0.006 R ⊙ , and 0.204 ± 0.004 R ⊙ for TOI-782, TOI-1448, TOI-2120, and TOI-2406, respectively.The updated values of T eff are consistent with those derived from the IRD spectra (Section 3.1.1)within 1σ for TOI-782, TOI-1448, and TOI-2120 and within 2σ for TOI-2406.The updated values of R s are all within 1σ of those derived from the empirical relation.Our values of T eff , R s , and d are listed in Table 4.The prior and posterior SED models along with the fitted SED data are shown in Figure 12.

Stellar Activity and Rotation
In order to investigate the magnetic activity and rotation periods of the host stars, we searched for long-term photometric variability of the host stars in the TESS light curves.We found no significant variability in any of the four stars.We also searched for periodic variability in the g-and r-band archival light curves of the Zwicky Transient Facility (ZTF) DR165 which are available for TOI-1448, TOI-2120, and TOI-2406, and in the g-and V -band light curves from ASAS-SN6 for TOI-782, by performing a Lomb-Scargle periodogram analy-  (Bakos et al. 2013) during two observing seasons between 2010 January 26 and 2013 July 23, with a gap in the observations between 2010 August 10 and 2012 December 18.A total of 25,247 photometric observations were obtained in the r-band with an exposure time of 4 minutes.The observations were reduced to trend-filtered light curves following the procedures described in Bakos et al. (2013).We used the generalized lomb-scargle periodogram (GLS) (Zechmeister & Kürster 2009) to detect a significant periodic signal in the light curve at a period of P = 64 days, and with a semiamplitude of 7.28 ± 0.23 mmag, assuming a strictly periodic sinusoidal signal.The signal has a formal bootstrap-calibrated formal FAP of less than 10 −200 , with the number of independent trials calibrated through a bootstrap procedure.This indicates a vanishingly small chance of white noise data producing a periodogram peak of this height.However, the signal also appears to vary over time.When we analyze the two seasons separately, the peak period in the GLS periodogram for the first season is at a period ∼ 30 days, with a lower significance, while the second season shows a strong peak at ∼ 60 days.The semi-amplitude of the P = 64 days signal is also lower in the first season, with a value ∼ 2 mmag.Because of the apparent evolution of the signal with time, we also calculated the Discrete Autocorrelation Function (DACF; Edelson & Krolik 1988) of the light curve, and find a clear cyclic variation with a period of P = 65.653± 0.051 days based on fitting two parabolas to the first and second prominent peaks in the DACF.Some less prominent peaks also appear at shorter time lags.The analysis discussed so far is for a partial detrending of the light curve against a set of auxiliary instrumental parameters (henceforth referred to as external parameter decorrelation (EPD)), which yields a point-to-point RMS of 25 mmag.When we apply the more aggressive trend-filtering algorithm (TFA; Kovács et al. 2005)  the HATSouth light curve, which reduces the RMS to 20 mmag, the signal remains at a high level of significance, but the amplitude of this long-period variation is reduced to 4.18 ± 0.22 mmag.The DACF of the TFA light curve appears to show a cleaner cyclic variation, with the lag between the first and second peaks indicating a period of P = 65.257± 0.063 days.This suggests that the additional less pronounced, shorter-lag peaks present in the EPD DACF are likely due to instrumental artifacts.The periodogram, DACF and light curve are shown in Figure 2.

Stellar Age
TOI-782 and TOI-1448 have a slow rotation with P ∼ 65.7 days and 42.3 days, respectively.Neither TOI-782 nor TOI-2120 show any activity variations due to stellar spots.The low-velocity stars, TOI-782, TOI-1448, and TOI-2120, are probably thin disk stars (see Table 4), whereas the total velocity of the metal-poor TOI-2406, , suggests that it is a member of the thick disk (see also Wells et al. 2021).In fact, the thick disk to thin disk probability ratios of TOI-782, TOI-1448, TOI-2120, and TOI-2406 are estimated to be 0.0125, 0.0144, 0.0213, and 38.9, respectively, from their Galactic kinematics (Bensby et al. 2014), respectively.Also, the BANYAN Σ analysis (Gagné et al. 2018) confirm that four these M dwarfs are not bona fide members of the young association within ∼100 pc.

Validation of the Planetary Candidates
In the following paragraphs, we qualitatively eliminate the classes of false positive scenarios that can mimic the transit signal including an EB with a grazing transit geometry, a hierarchical eclipsing binary (HEB), and a diluted eclipse of a background (or foreground) eclipsing binary (BEB) along the line of sight of the target.Although the planetary nature of TOI-2406 b has already been validated by Wells et al. (2021), we will still validate all four planetary candidates for completeness.
First, the renormalised unit weight error (RUWE) from Gaia DR3 7 is 1.25, 1.08, 1.15, and 1.06 for TOI-782, TOI-1448, TOI-2120, and TOI-2406, respectively, which are all low enough to be consistent with a single star being the source of the astrometric motion (Belokurov et al. 2020).We also rule out the EB scenario based on the mass constraint derived in Section 3.4.3.Moreover, we do not detect any wavelength dependence of the transit depth from our chromatic transit analysis (Section 3.4.2),which is consistent with the absence of contamination from a star of a spectral type (color) different from the host star.In the absence of dilution, the measured planetary radii become significantly smaller than the lower limit of 0.8R Jup expected for brown dwarfs (Burrows et al. 2011).Grazing transit geometries can be eliminated, as the impact parameter 7 https://gaia.ari.uni-heidelberg.de/singlesource.html is constrained to b < 0.6 at the 99% level based on our transit and contamination analyses.The apparent boxy shapes of the transit light curves are in stark contrast with the V-shaped transit expected for grazing orbits.Hence, the grazing EB scenario is ruled out.
Second, we constrain the classes of HEBs that 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 Mori et al. (2022), which is based on Bouma et al. (2020), to perform the calculation over a grid of eclipsing companions. 8Comparing the simulated eclipse depth with the observed one (Section 3.4.2) in each band, we find that no plausible HEB configuration explored in our simulation can reproduce the observed transit depths in the multiple bands simultaneously as well as the transit shape.To check the possibility of stellar companions spectroscopically, we visually searched for secondary peaks in the stellar line profiles.For all the four targets, the mean line profile, a cross-correlation function of the IRD spectra against a template spectrum, exhibited no sign of a stellar companion with a contrast higher than ∼ 0.1 in the near-infrared wavelengths.Hence, the HEB scenario is ruled out.
Third, the probability of being a BEB is low a priori because our targets are far from the Galactic plane, except for TOI-1448 and TOI-2120,9 and we argue below that the BEB scenario for TOI-1448 and 2120 is very unlikely.First, there is a background star (TIC 389900766; ∆G=5.1)near TOI-2120 with a separation of 2. ′′ 2-2.′′ 5 at the times of the ground-based transit observations, which contaminates the photometric apertures for these observations (3 ′′ or larger in radius).However, we did not detect any significant chromatic difference in the transit depth in the analysis of the ground-based transit light curves (Section 3.4.2),which cannot be explained by the BEB scenario given the chromatic flux ratios of the background star to TOI-2120 of 1.3%, 1.6%, 0.78%, and 0.79% for the g, r, i, and z s bands, respectively, measured from the Pan-STARRS1 catalog (Chambers et al. 2016).Next, our high-resolution speckle and AO imaging ruled out any other nearby star and blended sources down to 0. ′′ 5 at a ∆mag of 7.7 (2166 nm), 7.0 (2166 nm), 5.6 (832 nm), and 5.2 (832 nm) for TOI-782, TOI-1448, TOI-2120 and TOI-2406, respectively.We statistically estimate the probability of a chance-aligned star, using the population synthesis code TRILEGAL 10 ( Girardi et al. 2005), which simulates the Galactic stellar population along any line of sight.Given the positions of our targets, we found probabilities of P BEB < 10 −6 10 http://stev.oapd.inaf.it/cgi-bin/trilegal to find a star brighter than T mag ≃16,11 within 0. ′′ 5. Assuming that all such stars are binaries and preferentially oriented edge on to produce eclipses with the period and depth consistent with the TESS detection, this can represent a conservative upper limit to the BEB scenario.
Finally, we quantify the false positive probability (FPP) of our targets using the Python package TRICERATOPS (Giacalone & Dressing 2020).We used the ephemerides and depths reported in Table 5, the contrast curves from the high-resolution images with the highest contrast achieved, and the MuSCAT3 i-band follow-up light curves as inputs to TRICERATOPS.Although we were able to rule out the classes of EB, BEB, and HEB, we ran TRICERATOPS considering all of these scenarios for completeness.We obtained formal FPPs of the order of 10 −5 , 10 −4 , 10 −7 , and 10 −5 for TOI-782 b, TOI-1448 b, TOI-2120 b, and TOI-2406 b, respectively.These values are well below the threshold of FPP< 1.5% prescribed by Giacalone et al. (2021) which statistically validates them as bona fide planets.

Properties of the Planets
To derive the physical properties of the planets, we analyzed the TESS light curves, ground-based transit light curves, and RV data of each planetary system in the following steps.First, we modeled the TESS light curves with a transit + Gaussian Process (GP) model to determine the hyperparameters of the GP model.Next, we modeled the ground-based transit light curves with a transit + GP model using the posteriors of the transit parameters obtained from the above analysis as priors and then determined the hyperparameters of the GP model for the ground-based data.Finally, we jointly modeled all the light curves and RV data with transit and RV models along with the predetermined hyperparameters of the GP models to derive the final values and uncertainties of the planetary parameters.We describe each step in more detail in the following subsections.

TESS Light Curves
For the transit model, we adopted the Mandel & Agol model (Mandel & Agol 2002) implemented by PyTransit (Parviainen 2015) with the following parameters: scaled semimajor axis a R (= a/R s , where a is the semimajor axis), impact parameter b, planet-to-star radius ratio k (= R p /R s , where R p is the planetary radius), eccentricity e, the argument of periastron ω, orbital period P , individual transit times T c,i , and two coefficients u 1 and u 2 of a quadratic limb-darkening law.Note that the model is supersampled, i.e., averaged over the exposure time.At this stage, we set both e and ω to zero and P to the value provided by the TOI catalog.For u 1 and u 2 , we applied Gaussian priors with the values and uncertainties of the stellar parameters given by LDTk (Parviainen & Aigrain 2015).Note that we increased the uncertainties of u 1 and u 2 given by LDTk by a factor of 3 to account for the possible systematics in the stellar models.We assumed uniform priors for the other parameters.
Simultaneously with the transit model, we also modeled the time-correlated noise in the light curves using a GP model implemented in celerite (Foreman-Mackey et al. 2017) with a stochastically driven, damped simple harmonic oscillator (SHO) kernel function, where the model parameters are the frequency of the undamped oscillation ω 0 , the scale factor to the amplitude of the kernel function S 0 , and the quality factor Q. We set Q to unity in all sectors to avoid overfitting, while leaving ω 0 and S 0 free for each sector to account for the different noise properties from sector to sector.We also fit a white jitter term in the flux for each sector.We ran an MCMC analysis for each target using emcee to estimate the posterior distributions of the parameters.

Ground-based Light Curves
We simultaneously modeled the ground-based, multiband transit light curves of each target with a transit + GP model.We used the same transit model as described in Section 3.4.1,but left k free for each band in order to search for chromaticity in k, which can be a sign of flux contamination in the photometric aperture.The limb-darkening parameters u 1 and u 2 were fitted with Gaussian priors for each band in the same way as for the TESS light curves.Note that all the measured fluxes of TOI-2120 are contaminated by a nearby faint background star (TIC 389900766), as mentioned in Section 3.3.We, therefore, corrected the contamination to TOI-2120 b by replacing the transit model F tr with (F tr +f cont )/(1+f cont ), where f cont is the flux ratio of the target to the contaminating source (1.3%, 1.6%, 0.78%, and 0.79% for the g, r, i, and z s bands, respectively, assuming that all the fluxes were contaminated in the photometric apertures).The time-correlated noise in the ground-based light curves was computed by a GP model as a function of time in the celerite package whose covariance function is an approximate Matérn 3/2 kernel with signal standard deviation σ and length scale ρ.The kernel function has the same ρ in all light curves of each transit event but different σ.Note that this kernel is a specific version of the SHO kernel (Section 3.4.1)with no oscillation term.We applied this simplified form to the ground-based light curves considering that their observational duratios are shorter than the stellar variability timescale.The MuSCAT2 data exhibit nonnegligible systematics arising from telescope drifts.Therefore, for the time-correlated noise model of the MuSCAT2 data, we multiplied the GP model by a linear function of ∆x and ∆y, where ∆x and ∆y denote stellar displacements in the detector in the x-and y-directions, respectively.We performed an MCMC analysis for each target using emcee with uniform priors for all parameters except for u 1 and u 2 .

Curves
From a preliminary analysis of the RV data, we did not detect significant (> 3σ) planetary signals in any of the systems.Nonetheless, the RV data are still useful for placing upper bounds on the planetary masses and orbital eccentricities.Because both the RV and transit light-curve models depend on eccentricity, the best constraint on the mass and eccentricity can be obtained by a joint analysis of the RV and transit light-curve data.We therefore simultaneously modeled the RV data and all the available transit light curves (from both TESS and the ground) of each system.The parameters of the RV model include RV amplitude K, orbital inclination i, eccentricity parameters e and ω (fitted as √ e sin ω and √ e cos ω), orbital period P , time of periastron passage t 0 , RV offset v 0 , and RV jitter σ v,jit .We additionally included a parameter for a linear slope in the RV, γ, for the data of TOI-1448 b and TOI-2120 b, for which a preliminary fit of the RV data only indicated that a model with a linear slope was preferred over a model without a slope, determined by corresponding values of the Bayesian information criterion of 3.8 and 21.8, respectively.The parameters for the transit model are the same as in Section 3.4.2,but this time the assumption of a circular orbit was removed.In addition, we assumed a constant orbital period, which introduces a reference transit time T c,0 for each system, because we did not detect significant variations in the individual transit timings (T c,i ) measured in the TESS and ground-based light curves from a constant period in any system.We also adopted a constant radius ratio k across the bands for each system because we did not detect any significant chromaticity in k in the analyses of Sections 3.4.1 and 3.4.2.The hyperparameters of the GP for the TESS and ground-based light curves were set to the best-fit values derived in Sections 3.4.1 and 3.4.2,respectively, to suppress the number of free parameters.Note that the correlations between the hyperparameters and other physical parameters are small enough that fixing the hyperparameters at the best-fit values does not have much effect on the posterior distributions of the hyperparameters.
With the above parameterization, we performed an MCMC analysis using emcee.Note that in this analysis i and t 0 are not fitted but i is converted from b, a R , e, and ω by the equation of b = a R cos i(1−e 2 )/(1+e sin ω) and t 0 is converted from T c,0 , P and Kepler's equation.In addition, we use log stellar density log ρ s as the fitting parameter instead of a R , which is converted from ρ s and P by the equation a R = (ρ s GP 2 /3π) 1/3 (assuming k 3 ≪ 1), where G is the gravitational constant.
For the priors, we applied a two-sided Gaussian distribution for ρ s using the values derived in Section 3.1.2.We also applied the joint probability distribution for e and ω of Kipping (2014), Equation ( 23) in that paper), which is a conditional probability given that the planet has a transit geometry adopting a beta function as the underlying prior for e.For the two coefficients of the beta function, we adopted α = 1.58 and β = 4.4 from Van Eylen et al. (2019), which were derived from a sample of Kepler single-transiting systems.Uniform priors were applied to all the other parameters.
The derived posterior values (median and 1σ boundaries or 2σ upper limits) are listed in Table 5.The transit light curves detrended and phase folded by the best-fit parameters, and the time series and phase-folded RV data along with the posterior RV models, for TOI-782b, TOI-1448b, TOI-2120b, and TOI-2406b are shown in Figures 3, 4, 5, and 6, respectively.

Results
Combining the results of the joint analysis in Section 3.4.3with the stellar properties derived in Section 3.1.2,we derive the physical properties of the planets as listed in Table 5.We find that TOI-782 b, TOI-1448 b, TOI-2120 b, and TOI-2406 b have radii of 2.734 ± 0.085 R ⊕ , 2.749 +0.067 −0.063 R ⊕ , 2.122 ± 0.065 R ⊕ , and 2.860 +0.063 −0.069 R ⊕ , respectively, all of which are within the sub-Neptune regime.Our value for TOI-2406 b is consistent with that measured by Wells et al. (2021) within the uncertainties but is more precise thanks to the additional TESS data taken after their work and our high-precision multiband transit light curves taken with MuSCAT3.
Two-dimensional posterior distributions between M p , e, and ω obtained from the joint analysis are shown in blue in Figure 7.Although we are unable to measure the masses of any of the planets with a high enough significance, we place 2σ upper limits on the masses of 19.1, 19.5, 6.8, and 15.6 M ⊕ for TOI-782 b, TOI-1448 b, TOI-2120 b, and TOI-2406 b, respectively.We find tentative linear trends in the RV data of TOI-1448 and TOI-2120 with γ = 0.041 ± 0.019 and −0.044 ± 0.007 m s −1 day −1 , respectively, which could be due to an additional outer planet or companion, or some unrecognized systematics.Further RV observations are needed to determine the definitive masses of these planets and to confirm or refute the presence of outer planets/companions to TOI-1448 b and TOI-2120 b.
While we cannot unambiguously determine the eccentricity due to the lack of significant planetary signals in the RV data, our joint analyses of RV data and transit light curves with stellar density priors show that three of the four planets, TOI-782 b, TOI-2120 b, and TOI-2406 b, have nonzero eccentricity with significances of ∼2-4σ (0.19 +0.09 −0.11 , 0.32 ± 0.08, and 0.17 +0.11 −0.05 , respectively).These eccentricities are mostly constrained by the transit-light curves with the stellar density priors.To illustrate this, in Figure 7 we overplot the posterior distributions of M p , e, and ω obtained from analyses of the RV data alone (where T c,0 , P and i are fixed to the values measured in the transit light curves) in yellow and the posterior distributions of e and ω from analyses of the transit light curve data alone (both TESS and ground-based ones) in red.We note that our conclusion that three of the four planets are likely to have nonzero eccentricity does not depend on the prior on the eccentricity; we also detect nonzero eccentricities for TOI-782 b, TOI-2120 b, and TOI-2406 b at similar significances even with a uniform prior on e.We further note that a nonzero eccentricity of TOI-2406 b was originally found by Wells et al. (2021), with a uniform prior on e, and we confirmed it with almost independent transit data (except for TESS Sectors 3 and 30) and new constraints from the RVs.

DISCUSSION
In this study, we discovered and followed up four single sub-Neptunes transiting nearby (<100 pc) M dwarfs.Although these four planets are not outliers in the radius-period diagram among known close-in sub-Neptunes around M dwarfs (see Figure 9), we found that at least three of them are likely to have eccentric orbits with e ∼ 0.2 − 0.3, which are slightly larger than the majority of known close-in sub-Neptunes with orbital periods of ≲ 10 days (see Figure 10).One possible explanation for the nonzero eccentricities of these planets is that there is an unseen (long-period) planet (e.g.Limbach & Turner 2015) or a substellar companion (e.g.Ngo et al. 2016) in each system.Once the disk gas has dissipated, the eccentricity of a planet can 13.0 +4.0 −2.9 3.5 ± 2.0 4.7 +5.9 Note-The reported values, uncertainties, and upper limits represent the median, 1σ lower and upper boundaries, and 2σ upper boundaries of the posterior probability distributions, respectively.a Barycentric TESS Julian Date, which corresponds to Barycentric Julian Date (BJD) -2457000.be excited by other planets and (sub-)stellar companions through gravitational perturbations, such as planetplanet scattering (e.g.Weidenschilling & Marzari 1996) and the von Zeipel-Lidov-Kozai mechanism (e.g.Fabrycky & Tremaine 2007).In fact, TOI-1448 b and TOI-2120 b may have an unseen siblings as suggested by the linear trend in the RV data (see Figure 5).Another dynamical process that influences the eccentricity of a close-in planet is tidal interaction between the planet and its host star.As an eccentric planet ap-proaches to its host star, tidal forces circularize its orbit.Nonzero eccentricities of the close-in sub-Neptunes may suggest that tidal circularization of the planet's orbit is not achieved because the host star is too young or because inefficient tidal dissipation occurs in the planetary interior.The estimated ages of the four M dwarfs are all older than 1 Gyr (see Section 3.2).Therefore, all four sub-Neptunes should experience tidal forces from their host stars over 1 Gyr.The nonzero eccentricities of at least three of the sub-Neptunes are suggestive of less efficient tidal dissipation in their interiors.We consider the tidal evolution of thr three shortperiod sub-Neptunes in eccentric orbits around relatively old M dwarfs (see Table 6).A tidal interaction between a short-period planet and an M dwarf circularizes the orbit of the planet over gigayears.The eccentricity damping of a planet in equilibrium tide (e.g.Murray & Dermott 1999) is written as where e is the eccentricity of the planet, M p is the planetary mass, R p is the planetary radius, a is the semimajor axis, μp is the effective bulk modulus, Q p is the tidal dissipation factor of the planet, M ⋆ is the mass of the host star, and Ω is the Keplerian angular velocity.Fig- ure 11 shows the eccentricity damping timescales of the three sub-Neptunes.We adopt the upper limit masses of TOI-782 b, TOI-2120 b, and 2406 b.A nonzero eccentricity for each planet can persist over gigayears if the Q value is ≳ 10 3 for TOI-782 b and 2120 b, and ≳ 10 5 for TOI-2406 b, assuming that μp ∼ 1 and Q p is constant for 1 Gyr.A typical Q value for a terrestrial planet and a gas giant is 10 − 100 and 10 5 − 10 6 , respectively.We consider here that massive rocky planets may have a Q-value greater than 10 − 100.Massive rocky planets have magma oceans for a long period of time due to high pressures and temperatures in their interiors.Tidal damping in rocky super-Earths can be less efficient due to the size effect (Efroimsky 2012).Tidal dissipation in massive rocky planets may be weaker than that in terrestrial planets.On the other hand, the observed eccentricity of super-Earths suggests that their typical Q value is ∼ 10 − 100 (Hansen & Murray 2015), for example, Q p ∼ 100 for CoRoT-7 b (Clausen & Tilgner 2015).The interiors of these three eccentric sub-Neptunes with large Q values are different from those of close-in super-Earths like CoRoT-7 b.In fact, the mass-radius relationship of these four sub-Neptunes supports that they are not composed of pure rock (see Figure 8).The four validated planets, as well as K2-18 b (e.g.Madhusudhan et al.Other known eccentric, sub-Neptunes close to M dwarfs, K2-25 b with R p = 3.44 ± 0.12R ⊕ , M p = 24.7 +5.7 −5.2 M ⊕ , and e = 0.43 ± 0.05 in the Hyades star cluster (600-800Myr; (Mann et al. 2016;Stefansson et al. 2020)) and TOI-269 b with R p = 2.77 ± 0.12R ⊕ , M p = 8.8±1.4M⊕ , and e = 0.425 +0.082 −0.086 (Cointepas et al. 2021) around an old M dwarf with an age of a few gigayears, have drawn attention to the dynamical histories of close-in planets.K2-25 b and TOI-269 b can remain in an elliptical orbit until the present day if Q p ≳ 10 4 −10 5 (see Figure 11).These Q values larger than those of terrestrial planets are consistent with low bulk densities of both planets, indicating the presence of volatile material in their interiors.

SUMMARY
We report the discovery and follow-up of four sub-Neptunes with radii of 2 − 3R ⊕ and P ≲ 8 days orbiting M dwarfs (TOI-782, TOI-1448, TOI-2120, and TOI-2406), three of which were newly validated by groundbased follow-up observations and statistical analyses.Our RV follow-up observations with Subaru/IRD suggest that all of the four planets are less massive than 20M ⊕ at 2σ significance.The orbits of at least three short-period sub-Neptunes orbiting relatively old M dwarfs with ages ≳ 1 Gyr are not yet tidally circularized.The data suggest that these planets have eccentricities of e ∼ 0.2 − 0.3, which would require inefficient tidal dissipation in their interior.The slow damping of the eccentricities due to tidal interactions requires a large Q value of 10 3 − 10 5 for these three sub-Neptunes.The mass-radius relationship of all four sub-Neptunes suggests that they are unlikely to be pure rocky planets, that is, they have a significant atmosphere and/or an icy mantle on the core.The validated eccentric sub-Neptunes give us an unprecedented opportunity to infer the elusive interiors and formation histories of sub-Neptunes close to a central star.Wavelength (Å)

Figure 1 .
Figure 1.(a) Power spectra for the g-(top) and r-band (bottom) ZTF light curves of TOI-1448.The gray dashed lines indicate FAPs of 0.1, 1, and 10% from top to bottom.Significant, coincident signals as shown by the orange vertical lines) can be seen at a period of 42.3 days.(b) The ZTF light curves of TOI-1448 in the g band (top) and the r band (bottom) folded with the peak periods.Light dots and dark circles are unbinned and binned data, respectively.Black curves show the best-fit sine function.

Figure 2 .
Figure2.Top left: the HATSouth light curve of TOI-782 plotted against time.The grey points show the unbinned data, while the red squares show the data binned by 1 day.For clarity, we cut out a gap in the data as shown.The light curve shown in this figure has had an EPD detrending applied to it, but not the more aggressive TFA detrending that is typically also applied to HATSouth light curves.Top right: The GLS unnormalized periodogram of the HATSouth light curve of TOI-782 shows a strong peak at a period of P = 64 days.Middle left: The DACF of the HATSouth EPD light curve.The first and second prominent peaks are fitted with parabolas to determine the period.Middle right: same as on the left, but here we show the DACF of the HATsouth TFA light curve.Some shorter-lag variations in the EPD DACF are not present, suggesting that they are likely due to instrumental systematics.Bottom left: the HATSouth EPD light curve phase folded at the detected period.The gray-scale points show the unbinned data, while the red squares show the phase-binned data with a bin size of 0.02.Bottom tight: same as on the left, but here we show only the phase-binned data and zoom in to see the detected photometric variation more clearly.

Figure 3 .
Figure3.(a) Detrended and phase-folded light curves of TOI-782b (upper five curves) and their residuals from the best joint-fitting model (lower five curves).The black, blue, green, orange, and red points (from top to bottom in each five data sets) indicate 5 minutes binned data from the TESS, g, r, i, and zs bands, respectively.Gray lines indicate the best-fit models.The plots are vertically shifted for display.(b) Top: time series of RV data from Subaru/IRD (blue) along with the 1σ, 2σ, and 3σ credible regions calculated from the posteriors of the joint-fitting analysis (from thick to thin orange).Error bars represent the estimated 1σ uncertainties without a jitter term.Middle: residuals of the RV data from the best joint-fitting model.Bottom: same as the top panel, but phase folded with the orbital period of TOI-782 b.Areas from dark to light orange indicate 1, 2, and 3σ credible regions, respectively.

Figure 5 .
Figure 5. Same as Figure 3 but for TOI-2120 b.A linear trend in the RV data was injected in the joint analysis of TOI-2120.

Figure 7 .
Figure 7. (a) Corner plots for the posterior samples of Mp, e, and ω for TOI-782 b.Yellow, red, and blue contours indicate the posteriors from the RV-only fit, the transit-only fit, and the RV+transit joint fit, respectively.Contours with thick to thin solid lines indicate 68%, 95%, and 99% confidence regions, respectively.The density of the posterior samples from the joint analysis is indicated by the color map.(b), (c), (d) Same as (a) but for TOI-1448 b, TOI-2120 b, and TOI-2406 b, respectively.

Figure 8 .Figure 9 .
Figure8.Mass-radius relationship of planets with measured radii of 2 − 4R⊕ and orbital periods of ≲ 50 days around Sun-like stars and M dwarfs (data from http://www.exoplanet.eu).The curves show the interior models calculated for planets composed of pure H2O, MgSiO3, and Fe.We adopted the third-order Birch-Murnagham equation of state (EOS) for MgSiO3 perovskite(Karki et al. 2000;Seager et al. 2007), the Vinet EOS for ϵ-Fe(Anderson et al. 2001), the Thomas-Fermi-Dirac EOS(Salpeter & Zapolsky 1967) at high pressure(Seager et al. 2007;Zeng & Sasselov 2013), and AQUA EOS for H2O(Haldemann et al. 2020).We also plotted the contours of mass ranges derived from joint analyses of the four sub-Neptunes validated in this study.

Figure 10 .
Figure 10.Eccentricity distribution of planets with measured radii of 1.5 − 4R⊕ and orbital periods of ≲ 50 days around M dwarfs with T eff = 3000 − 3900 K.Only planets with measured eccentricities are shown.The color contour represents the effective temperature of a planet-hosting M dwarf.The four close-in, sub-Neptunes validated in this study are shown as red stars.

Figure 11 .
Figure 11.The eccentricity damping timescale of the three sub-Neptunes: TOI-782 b (red), TOI-2120 b (blue), and TOI-2406 b (orange).We adopted the upper limit values of the planetary masses for TOI-782 b, TOI-2120 b, and TOI-2406 b.The expected tidal dissipation factors of two other eccentric sub-Neptunes, K2-25b (black) and TOI-269 b (gray), are also shown for comparison.The vertical dashed lines correspond to μpQp of planets with τe = 1 Gyr.

APPENDIXA.
SPECTRAL ENERGY DISTRIBUTIONSThe SEDs of the four host stars are shown in Figure12.
B. TESS LIGHT CURVES OF THE INDIVIDUAL SECTORSThe TESS PDC-SAP light curves of each sector for TOI-782, TOI-1448, TOI-2120, and TOI-2406 are shown from Figures 13 to 13 .
C. INDIVIDUAL TRANSIT LIGHT CURVES FROM GROUND-BASED OBSERVATIONSThe individual transit light curves of the four planets obtained with the ground-based instruments are shown in Figure17.

Figure 12 .
Figure 12.(a) SED of TOI-782.Panels (b), (c), and (d) are the same as (a), but for TOI-1448, TOI-2120, and TOI-2406, respectively.Yellow and blue curves show the prior and posterior SED models.Red diamonds and blue squares are observational data and best-fit results, respectively.

Table 1 .
Observation Log of Transit Photometry from the Ground

Table 2 .
Radial Velocities of the Host Stars Measured with Subaru/IRD

Table 3 .
Observation Log of High-contrast Imaging Note-SAI: Sternberg Astronomical Institute adopted these values as our fiducial values for M s , which are listed in Table

Table 4 .
Properties of the Four M Dwarf Hosts .As a result, there were no plausible periodic signals in any of the light curves except for the g-and r-band ZTF light curves of TOI-1448, in which we found significant and coincident periodic signals with false alarm probabilities (FAPs) of less than 0.1%, as shown in Figure1.The strongest signals in the power spectra for the g-and r-bands are present at periods of 42.32 ± 0.11 d and 42.22 ± 0.13 days, having variability amplitudes of 10.8 ± 1.8 mmag and 9.2 ± 1.5 mmag, respectively.Considering the significance and period coincidence of the two signals, we attribute these signals to stellar rotation.The weighted mean of the two periods is calculated to be 42.28 ± 0.085 days, which is listed in Table4.TOI-782 was observed by the HATSouth groundbased telescope network a From Gaia DR3.The values are given for epoch=J2016.0.sis

Table 5 .
Posterior Values of the Parameters Used in the Joint Fitting

Table 6 .
Properties of the Four Sub-Neptunes around M dwarfsNote-The reported values, uncertainties, and upper limits represent the median, 1σ lower and upper boundaries, and 2σ upper limits of the posterior probability distributions, respectively.
Hori et al.