The TESS–Keck Survey. XIX. A Warm Transiting Sub-Saturn-mass Planet and a Nontransiting Saturn-mass Planet Orbiting a Solar Analog

The Transiting Exoplanet Survey Satellite (TESS) continues to increase dramatically the number of known transiting exoplanets, and is optimal for monitoring bright stars amenable to radial velocity (RV) and atmospheric follow-up observations. TOI-1386 is a solar-type (G5V) star that was detected via TESS photometry to exhibit transit signatures in three sectors with a period of 25.84 days. We conducted follow-up RV observations using Keck/High Resolution Echelle Spectrometer (HIRES) as part of the TESS–Keck Survey, collecting 64 RV measurements of TOI-1386 with the HIRES spectrograph over 2.5 yr. Our combined fit of the TOI-1386 photometry and RV data confirm the planetary nature of the detected TESS signal, and provide a mass and radius for planet b of 0.148 ± 0.019 M J and 0.540 ± 0.017 R J, respectively, marking TOI-1386 b as a warm sub-Saturn planet. Our RV data further reveal an additional outer companion, TOI-1386 c, with an estimated orbital period of 227.6 days and a minimum mass of 0.309 ± 0.038 M J. The dynamical modeling of the system shows that the measured system architecture is long-term stable, although there may be substantial eccentricity oscillations of the inner planet due to the dynamical influence of the outer planet.


INTRODUCTION
The plethora of exoplanet discoveries is primarily a result of the transit and radial velocity (RV) methods, the combination of which account for ∼95% of the currently known inventory, according to the NASA Exoplanet Archive (NEA) (Akeson et al. 2013).The Transiting Exoplanet Survey Satellite (TESS) is an all-sky photometric survey searching for transiting planets around the nearest and brightest stars (Ricker et al. 2015).At the time of writing there are 415 confirmed TESS planets, with 7027 candidates awaiting confirmation 1 .The observational resources required for follow-up RV observations are significant (Kane et al. 2009;Cloutier et al. 2018;Dalba et al. 2019;Dragomir et al. 2020;Guerrero et al. 2021;Kane et al. 2021b) and a coordinated RV campaign has been undertaken since TESS launched in 2018.The TESS-Keck Survey (TKS) is a multiinstitutional collaboration focused on planetary occurrence rates, formation, evolution, and dynamics (Chontos et al. 2022) and has directly confirmed numerous TESS candidates through precise mass measurements (Dalba et al. 2020b;Weiss et al. 2021;Dai et al. 2020;Rubenzahl et al. 2021;Scarsdale et al. 2021;MacDougall et al. 2021;Dalba et al. 2022;Lubin et al. 2022;Dai et al. 2021;Turtelboom et al. 2022).TKS primarily utilises the HIRES spectrometer on the Keck I Telescope at the W.M. Keck Observatory on Maunakea (Vogt et al. 1994) and the Levy spectrometer on the Automated Planet Finder (APF) (Radovan et al. 2014;Vogt et al. 2014) to confirm and characterize TESS objects of interest (TOIs).As TESS continues to survey the sky and TESS planets are confirmed, the planets discovered will provide answers to some of the remaining questions regarding demographics of planet populations.By surveying our nearest bright stars, the planets discovered will give a greater understanding of planet populations for a variety of stars compared to other surveys, such as Kepler, which focused on sun-like stars (Borucki et al. 2010).This will allow for validation (or invalidation) of the demographics found in the planet populations across all stellar types, such as the existence of the sub-Saturn valley (Ida & Lin 2004) and as a function of galactic latitude (Zink et al. 2023).
TOI-1386 (TIC 343019899) is a relatively bright (V = 10.51)solar-type star, for which a transit signal was detected during Sector 16 of TESS observations.This planet, designated TOI-1386 b, was adopted into the TKS single-transit program, and subsequent RV monitoring of the star commenced soon afterward.Early modeling of the transit and RV data indicated a ∼30 day orbital period for the planet, and a size and mass that places it within the category of a warm sub-Saturn.Such planets are an important component of exoplanet demographical studies, and can fall within important gaps in the known exoplanet population (Ford 2014;Winn & Fabrycky 2015;Fulton et al. 2017;Dulz et al. 2020;Hill et al. 2023;Ostberg et al. 2023).Study of the system was further motivated by the prospect of additional, non-transiting, giant planets that may lie in the system, the demographics of which can inform the early migration history and the influence on planetary architectures (Trilling et al. 1998;Alibert et al. 2005;Bonomo et al. 2017).Additional RV monitoring of the system successfully revealed another giant planet in the system, designated TOI-1386 c, whose orbital period of ∼227 days and moderate eccentricity creates an interesting dynamical environment.Multi-planet giant systems, such as TOI-1386, continue to contribute toward the inventory of exoplanet architectures that diverge significantly from that seen in the solar system (Horner et al. 2020;Kane et al. 2021c;Mishra et al. 2023a,b).
Here, we report on a detailed analysis of new measurements for the TOI-1386 system, confirming the planetary nature of the transit signals for the inner planet, detecting the presence of a non-transiting outer planet, and measuring the masses and orbital characteristics for both planets.Section 2 details the various data sources used in this analysis, including TESS photometry, high resolution imaging, and Keck RVs.Section 3 describes the properties of the host star, based on a combination of literature data and analysis of Keck spectra.Section 4 presents the results of our data modeling, including the Keplerian orbits for the two discovered planets and a dynamical simulation of their gravitational interactions.Section 5 discusses the implications of the results, particularly for understanding exoplanet demographics in the sub-Saturn valley, and we provide concluding remarks in Section 6.

Photometry
TESS observed TOI-1386 at a 30-minute cadence in Sectors 16 and 17 from 2019-Sep-11 to 2019-Nov-02, and then again at 2-minute cadence in Sectors 56 and 57 from 2022-Sep-01 to 2022-Oct-29.Thus, the star was observed for a total of 4 Sectors, providing photometric data over a period of ∼ 3.1 years.We initially downloaded data processed through the Quick Look Pipeline (QLP; Huang et al. 2020a,b) via the Mikulski Archive for Space Telescopes (MAST) portal (STScI 2018).2Only a single transit event was initially detected in Sector 16 through the QLP, though when looking at the raw data a second transit was detected at the beginning of Sector 17.We accessed the Pre-search Data Conditioning Simple Aperture Photometry (PDC-SAP; Stumpe et al. 2012;Smith et al. 2012) through MAST, cleaned the TESS photometry by keeping only points with quality flag = 0, indicating that there is no known problem or condition with the data, and stitched together and detrended the light curves from TESS sectors 16, 17, 56, and 57 into a single time-series using the Lightkurve software package (Lightkurve Collaboration et al. 2018).A single transit each was detected in TESS sectors 16, 17, and 56, resulting in an estimated orbital period for TOI-1386 b of ∼26 days.No transits were detected in Sector 57, as the predicted transit falls within the data gap for that sector.The top four panels of Figure 1 show the TESS photometry for the sectors during which the star was observed.The bottom panels show the detected transit signatures (bottom left) and folded transit (bottom right), where the red line indicates the best fit model to the transit.

Adaptive Optics and Speckle Imaging
We acquired adaptive optics (AO) imaging of TOI-1386 to search for nearby neighboring stars that may have affected or diluted the transit signal discovered by TESS.AO observations of TOI-1386 occurred on 8 Nov 2019 with Gemini/NIRI (Hodapp et al. 2003) 3 .At 0.2 ′′ , which at a distance of 146.86 pc (Gaia Collaboration et al. 2021) corresponds to a projected separation of ∼ 29 AU, we achieve a detection sensitivity of ∆mag > 4, and at 1 ′′ of separation we achieve a ∆mag > 7. A visual neighbor with a magnitude difference of ∆Kmag = 1.25 ± 0.01 and separation of ∼10.5 ′′ can be seen in the AO image (Figure 2).An additional visual companion was detected with a separation of ∼7 ′′ with a ∆Kmag = 7.23 ± 0.13.Extracting the details of these stars from Gaia Collaboration et al. ( 2021) we find the ∆Kmag = 1.25 ± 0.01 companion is a 15.2 R ⊙ star at a distance of 2924.1 pc.There are two Gaia sources ∼7 ′′ from TOI 1386, one is 10 magnitudes fainter in Gaia G magnitude and is likely a background source.The other, with a Gaia G magnitude difference of ∼ 8, has a parallax difference of 7.6 mas.Using the method outlined in Ciardi et al. (2015) to compensate for the flux dilution of neighbouring stars, we calculated that the ∆Kmag = 1.25 ± 0.01 companion could cause an error in the measurement of the planet radius of 8.7%.The ∆Kmag = 7.23 ± 0.13 companion made only a negligible difference of 0.1% of the measured radius of the transiting planet.
Speckle imaging was also acquired with the NESSI speckle imager (Scott et al. 2018;Scott 2019) mounted on the 3.5 m WIYN Telescope at Kitt Peak National Observatory on November 17, 2019.NESSI 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 follow- ing the procedures described in Howell et al. (2011).We were able to determine the limiting magnitude difference (∆mag) in both the blue (562 nm) and red (832 nm) band as a function of angular separation and found the speckle imaging from WIYN confirmed that no companions are visible from 0.05 -1.2 ′′ .Both the AO and speckle images rule out the possibility that the transit signal observed by TESS was caused by a blended stellar companion or eclipsing binary down to the detection limit of the imaging data.

Radial Velocities
We collected 64 spectra of TOI-1386 with the High Resolution Echelle Spectrometer (HIRES) instrument at the Keck Observatory (Vogt et al. 1994) between UT 2019-12-15 andUT 2022-06-11, including one high signal-noise-ratio (SNR) iodine-free template spectrum.The remaining observations included a gaseous iodine cell which was placed in the light path of the spectrometer.This allowed molecular absorption lines to be combined with the observed stellar spectrum.These lines were used as a wavelength reference for measuring the relative Doppler shift of each spectrum taken during our observations (Valenti et al. 1995;Butler et al. 1996).We reduced the spectra via the procedure outlined by the California Planet Search (CPS; Howard et al. 2010).The time series RVs are provided in Table 1 and shown in the top panel of Figure 4. Our RVs have a median nightly binned uncertainty of 1.71 m s −1 and a median SNR of ∼203 per pixel at the iodine wavelength region of ∼500 nm.Our 2.5 years of observations found both the transiting planet signal of 25.8 days along with an additional longer period, 227.6 day eccentric (e = 0.27) Keplerian signal.We also detect a long-term acceleration of the star.These discoveries are discussed further in Sections 4.1 and 4.3.

HOST STAR PROPERTIES
We analyzed our high S/N iodine-free HIRES spectrum of TOI-1386 with the SpecMatch code (Petigura et al. 2017b) to derive the T eff (5828.39±100.00K), log g (4.39±0.10cm s −2 ), and metallicity [Fe/H] (0.16±0.06 dex) of the host star.These results were used as input parameters for the EXOFASTv2 fit described in this section and in Section 4.3.
We derived the stellar mass, radius, and age using the isoclassify package (Huber et al. 2017;Berger et al. 2020).The analysis by isoclassify incorporated Gaia We used EXOFASTv2 to model the spectral energy distribution (SED) of the star together with the Gaia EDR3 parallax in order to measure the stellar radius.EXOFASTv2 computes SEDs from MIST models and Gaia parallax measurements and fits them to archival photometry to infer the properties of the host star.We used archival broadband photometry of TOI-1386 from Gaia DR2.We placed a upper bound on extinction (A V ≤ 1.841) from the galactic dust maps of Schlafly & Finkbeiner (2011).We used the SpecMatch analysis of TOI-1386 for the [Fe/H] and T eff priors.This provided a calculated stellar radius of R * = 1.015 +0.028 −0.026 R ⊙ .The values derived by isoclassify agree with those from EXOFASTv2 fit which are provided in Table 2 Reconnaissance spectra of TOI-1386 was taken on December 4, 2021 UT and December 5, 2021 UT with the Tillinghast Reflector Echelle Spectrograph (TRES; Furesz 2008) located at the Fred Lawrence Whipple Observatory (FLWO) in Arizona, USA.We compared the stellar parameters based on the Stellar Parameter Classification (SPC; Buchhave et al. 2012) analysis of six TRES spectra with our results from EXOFASTv2 and found the values agreed within the uncertainties.TOI-1386 is a chromospherically inactive star with log R ′ HK = -5.00± 0.05.Analysis of the TESS photometry found no significant periodicities that would imply the transit signal was caused by stellar activity.We computed the Ca II H&K index (S HK ) as described in Isaacson & Fischer (2010) for our Keck/HIRES time series data (Table 1).We found there was a correlation between the S HK values and RVs.Section 4.2 outlines our methods to test the robustness of the RV planet detection and rule out stellar activity as the cause of the RV signal.

Initial Spectroscopy Model with RadVel
The RV data were initially fit using the RV modeling toolkit RadVel (Fulton et al. 2018).RadVel enables users to model Keplerian orbits in radial velocity time series.RadVel uses an iterative approach to determine the best-fit orbital parameters for the observed RV curve.It then employs modern Markov chain Monte Carlo (MCMC) sampling techniques (Metropolis et al. 1953;Hastings 1970;Foreman-Mackey et al. 2013) to infer orbital parameters along with their associated uncertainties.
The RV data was fit with RadVel in order to model the ∼28 day period that was predicted from the initial single transit duration via the method outlined in Seager (2010).Early fits detected both the transiting planet period of 25.8 days along with a linear trend, which later resolved into the outer planet period of 227.6 days.The final preferred fit with RadVel included both planets b & c and a 4-sigma linear trend.This preferred model was confirmed in the final combined model fit we completed with EXOFASTv2 (Table 3 and Figure 4).

Testing the RV Model
As mentioned in Section 3, we found a Pearson correlation coefficient of 0.62 between the S HK values and RVs, with the S-value periodogram peaking at 234 days.To determine if the outer planet detection was robust we conducted a series of tests that would rule out stellar activity as the cause of the signal.We decorrelated the RV data by fitting a line of best fit to the S-values, removing this signal from the RVs and remodeling our best fit solution.We fit the decorrelated RVs and found the best fit had negligible changes from the original.We then added a planet to the RV fit with the period set to be equal to the value of the peak in the S-values and ran the model again.The best fit parameters for TOI-1386 b and TOI-1386 c did not change.We also ruled out sidereal day, lunation period, and sidereal year aliases as potential causes of the RV signal.We conclude from these tests that the RV signals are not an alias, nor are they caused by stellar activity.
As a typical stellar activity cycle is much longer than 234 days and stellar rotation rates tend to be shorter, we are unsure as to the cause of this peak in the s-values.We encourage further study into the possible causes of a peak at this period.

Combined System Modeling with EXOFASTv2
We modeled the stellar and planetary parameters for the TOI-1386 system using the EXOFASTv2 modeling suite (Eastman et al. 2013(Eastman et al. , 2019)).We included the TESS transit photometry along with the RVs from Keck-HIRES.The 30 minute cadence data from Sectors 16 & 17 with the 120 second cadence data from Sectors 56 & 57 were detrended using the spline fitting tool keplersplinev24 .For the 30-minute cadence we use the oversampling feature of EXOFASTv2, as 30 minutes is long relative to the ingress and egress times (τ ∼ 18 minutes).The oversampling feature samples the model at 10 points between ± 15 minutes of that timestamp, then numerically integrates the model over that range to get the single model point.This is important to avoid misfitting the transit shape.
When running the joint fit we placed normal priors on T eff (5828.39±100.00 K) and [Fe/H] (0.16±0.06 dex) from SpecMatch.The SED model discussed in Section 3 was simultaneously fit along with the transit and RV models.We found the preferred model included 2 planets as well as a 4-sigma linear trend ( γ) (Table 2).The best-fit transit model is shown in Figure 1 along with the data taken by TESS.Only three transits are depicted, as the transit that was set to occur in Sector 57 coincided with a data gap (Figure 1, bottom left).All RV data and the best-fit RV model are shown in Figure 4.The RV and the transit fits agreed, giving the interior planet, TOI-1386 b, a period of 25.8384 days and eccentricity of 0.061 (consistent with 0).The transiting planet has a radius of 0.540 R J and a measured mass of 0.148 M J , equating to a bulk density of 1.16 g cm −3 .This places TOI-1386 b in the sub-Saturn regime.The outer non-transiting planet, TOI-1386 c, was fit with a 227.6 day period and eccentricity of 0.27.This planet has a minimum mass (M p sin i) of 0.309 M J , making it a warm Saturn mass planet with an equilibrium temperature (T eq ) of 327.4 K (assuming no albedo and perfect temperature redistribution).The preferred fit also included a 4-sigma linear trend ( γ) which is an indication of an additional orbiting body, and further observations of TOI-1386 are needed to determine if the cause of this trend is due to an outer sub-stellar companion.The values included in this discussion and their uncertainties, along with many additional derived stellar and planet parameters provided by EXOFASTv2, are listed in Tables 2 and 3.

System Dynamics
The orbital eccentricity of planet c raises interest concerning the dynamical history of the system.To investigate this, we conducted a dynamical analysis using the N-body simulation package REBOUND (Rein & Liu 2012) with the symplectic integrator WHFast (Rein & Tamayo 2015).Planetary orbital parameters from Table 3 are taken as the starting condition of the simulation with the assumption that both planets are coplanar.Note that mutual inclinations between planets can result in non-linear trends regarding stable locations within the system, particularly for cases of mean motion resonance (Barnes et al. 2015), but can also resolve cases of orbital instability that result from coplanar assumptions (Kane 2016).Given the age of the system, the mutual inclinations present within this system would present a useful case study on the long-term stability of non-coplanar orbital configurations.
Under the assumption that both planets are coplanar, the minimum mass derived for planet c can be considered to be its true mass for the purpose of this study.The simulation was integrated for 10 7 years with a time step size of 1/20 of the planet b orbital period, or ∼1.3 days, in accordance with the recommended time step from Duncan et al. (1998) to ensure proper orbital sampling.We recorded the results every 100 years and show the time series eccentricity evolution of the two planets in Figure 5.The system is dynamically stable for the duration of the simulation.However, both planets experience eccentricity variations over time through the transfer of angular momentum (Chambers et al. 1996;Rasio & Ford 1996;Laughlin & Chambers 2001;Jurić & Tremaine 2008;Kane & Gelino 2013;Kane et al. 2021a).This is especially the case for the orbit of planet b, which varies from nearly circular to almost 0.2 eccentricity, suggesting that present observations happen to be occurring during a period of low eccentricity for the inner planet.Such large variation appears to be due to the influence of the outer planet, which has double the mass of the inner planet, driving the eccentricity evolution of both planets that exhibit in-sync periodic cycles every ∼25,800 years.Such an orbital configuration may remain longterm stable, or produce a planet-planet scattering event that leaves one remaining planet in a highly eccentric orbit (Chatterjee et al. 2008;Carrera et al. 2019;Kane et al. 2023).The relatively large separation of the two known planets makes a planet-planet scattering scenario unlikely, but the configuration of these planets and their dynamical interactions may exclude the presence of ad-ditional, smaller, planets within the inner part of the planetary system.

DISCUSSION
The size and mass of the inner transiting planet place it in the sub-Saturn valley (Figure 6).As described in Hill et al. (2023), the existence of the sub-Saturn valley is still a topic of debate.Some suggest the sub-Saturn gap could be attributed to core accretion theory (Ida & Lin 2004;Mayor et al. 2011;Emsenhuber et al. 2021), where planets that reach 10 M ⊕ enter a runaway accretion period and rapidly grow to ≥ 100 M ⊕ , provided there are sufficient materials available.Others, including Suzuki et al. (2018), argue the apparent gap is merely credited to the difficulty of detecting sub-Saturn mass planets and hence the discovery of TOI-1386 b is an important contribution to this demographic of planets.A study from Bennett et al. (2021) found that the sub-Saturn valley was missing from their analysis of RV planets observed using CORALIE/HARPS.Amongst those Sub-Saturns planets that have been detected there is a large density diversity, with sub-Saturn planet masses ranging from 6 -60 M ⊕ regardless of size Petigura et al. (2017a).TESS will help to fill out the mass-radius diagram by allowing for RV follow up of many transiting planets around the brightest stars in the sky.The inner transiting planet TOI 1386 b will contribute to the data used to ultimately verify the existence of the sub-Saturn valley and give clues as to the formation processes of these elusive planets.
As TOI-1386 b is a relatively rare sub-Saturn planet we calculated the transmission spectroscopy metric (TSM) value (Kempton et al. 2018) to determine if it is a good JWST target (Gardner et al. 2006).The TSM is proportional to the expected transmission spectroscopy signal-to-noise, based on the strength of spectral features, brightness of the host star, and mass and radius of the planet assuming cloud-free atmospheres.We calculate a TSM value of 42.5 for TOI-1386 b.We compared the parameters of this target to the population of known sub-Saturn planets (Figure 7).We limit the sub-Saturn valley planets to those with masses between 6 -60 M ⊕ in accordance with Petigura et al. (2017a).We find that compared to the rest of the sub-Saturn population TOI-1386 b is a massive, cool, moderately eccentric planet orbiting a metal rich sun-like star.In investigating the uniqueness of TOI-1386 b we found that a large number of detected sub-Saturn planets are a part of multiplanet systems that include multiple sub-Saturn planets within the system.When these rarely found planets are discovered, they are often found in duplicate.Of the 621 sub-Saturns listed on the NEA with a mass between 6-60 M ⊕ , there are 400 that are multi-planet systems.Of those systems we found 275 sub-Saturn planets were in systems where there were multiple sub-Saturn planets.This indicates that ∼ 44% of all sub-Saturn planets detected have additional sub-Saturn planets orbiting in their system.This finding agrees with the Weiss et al. (2018) study that found that in systems with multiple planets, the planets are likely to be of similar sizes.In the absence of any correction for detection biases and completeness, we are unable to provide a firm conclusion regarding this finding.More investigation into the sub-Saturn valley population needs to occur to determine whether this is a true phenomena or due to observational biases.
Originally chosen as part of the TKS single transit target list, TOI-1386 is a prime example of the importance of examining TESS raw data (Dragomir et al. 2020;Dalba et al. 2020a;Dalba et al. 2022).Investigation of the raw light curves from TESS revealed additional transits that may have initially been missed due to their proximity to data gaps.The discovery of these additional transits plus early RV data collection allowed quick identification of the ∼ 26 day orbit of the transiting planet.Although TOI-1386 was already earmarked for further investigation by the TKS team as a single transit target, the process of scrutinizing the raw data and detecting these extra transits reinforced the likelihood of the candidate being a planet and justified further study.This may be particularly important for planets that belong to key demographics, such as the sub-Saturn valley.The apparent scarcity of planets in this group may be due to the inherent challenges in detecting them, underscoring the importance of thorough data analysis.

CONCLUSIONS
TESS observed TOI-1386 in 4 Sectors (16, 17, 56 & 57) over a span of 3.2 years.Though initially only a single transit was detected, analysis of TESS data found 3 transits across the 4 Sectors.RV observations of TOI-1386 confirmed the existence of TOI-1386 b with a period of 25.83839 days and gave a mass measurement of M p = 0.148 M J , placing the 0.540 R J planet in the sub-Saturn regime.A nearby visual neighbor with a separation of ∼10.5 ′′ and magnitude difference of ∆ K mag= 1.25 ± 0.01 was calculated to have a dilution factor that could cause the radius of the planet to be underestimated by 8.7%.Even with this factor considered TOI-1386 b remains in the sub-Saturn regime and is an important addition to this demographic of planets due to the ongoing debate regarding the origin of the sub-Saturn gap.With a TSM of 42.  is a relatively massive, cool, moderately eccentric sub-Saturn planet orbiting a metal rich sun-like star.
RV analysis also detected an additional planet signal of an outer, non-transiting planet.TOI-1386 c is a eccentric (0.27) planet with a 227.6 day orbit and a minimum mass of 0.309 M J .Even with an eccentric outer planet, dynamical simulations of the system found that the two planets remained stable for the entire simulated period of 10 7 years.Both planets do, however, experience large changes in eccentricity over the simulation pe-riod with the inner planet eccentricity fluctuating from 0 to 0.2 and back every ∼ 25, 800 years.The best fit model of TOI-1386 also included a 4-sigma linear trend.Further observations of the star will be required to determine whether the cause of this trend is planetary in nature.
In our investigation of TOI-1386 b and the sub-Saturn valley population, we found that ∼ 44% of all identified sub-Saturn planets are part of systems containing other sub-Saturn planets.This finding suggests a tendency for planets of this size to coexist with similarly sized planets in their orbits.While this statistic lacks correction for detection biases and completeness, we encourage further research into the sub-Saturn valley to determine the validity of this phenomenon.
TOI 1386 is scheduled to be observed again by TESS in Sectors 76, 77 and 83.These observations, along with any additional RV observations, will help to refine the orbital parameters of the TOI 1386 system and determine the cause of the linear trend detected.We thank the time assignment committees of the University of California, the California Institute of Technol-ogy, NASA, and the University of Hawaii for supporting the TESS-Keck Survey with observing time at Keck Observatory and on the Automated Planet Finder.We thank NASA for funding associated with our Key Strategic Mission Support project.We gratefully acknowledge the efforts and dedication of the Keck Observatory staff for support of HIRES and remote observing.We recognize and acknowledge the cultural role and reverence that the summit of Maunakea has within the indigenous Hawaiian community.We are deeply grateful to have the opportunity to conduct observations from this mountain.We thank Ken and Gloria Levy, who supported the construction of the Levy Spectrometer on the Automated Planet Finder.We thank the University of California and Google for supporting Lick Observatory and the UCO staff for their dedicated work scheduling and operating the telescopes of Lick Observatory.This paper is based on data collected by the TESS mission.Funding for the TESS mission is provided by the NASA Explorer Program.This research has made use of the Exoplanet Follow-up Observation Program (ExoFOP; DOI: 10.26134/ExoFOP5) website, which is operated by the California Institute of Technol- Figure 1.Top: TESS full lightcurve for Sectors 16, 17, 56 and 57.Sectors 16 and 17 are 30-min cadence, processed from the FFIs and sectors 56 and 57 are 2-min cadence processed by the SPOC pipeline.TESS discovered 3 transits of TOI-1386 in (from top) Sectors 16, 17 and 56.The best fit period for the transiting planet positions a transit in the data gap of Sector 57 and hence a 4th transit was not recovered.Bottom Left: The 3 transits of TOI-1386 discovered by TESS in (from top) Sectors 16, 17 and 56.The red lines show the best fit model to each individual transit by EXOFASTv2.Bottom Right: EXOFASTv2 combined these transits with the RV data to determine the best fit model (red line).The residuals to the fit are shown below the transit model.

Figure 2 .
Figure 2. Left: The Gemini/NIRI Contrast curve and AO image (inset).The curve shows the limiting magnitudes (∆mag) for the non-detection of a neighboring star.At 1" of separation we achieve a ∆mag > 7. Two visual neighbors can be seen in the AO image, one with a magnitude difference of ∆ K mag= 1.25 ± 0.01 and separation of ∼10.5 ′′ .The other has a separation of ∼7 ′′ and a ∆Kmag = 7.23 ± 0.13 (Shown with a different color scaling).Both of these neighboring stars have a parallax difference to TOI 1386 of ≥ 7 mas (Gaia Collaboration et al. 2021).Right: Limiting magnitudes for the nondetection of a neighboring star based on the speckle imaging from the NESSI speckle imager on the WIYN Telescope.
DR2 parallaxes (Gaia Collaboration et al. 2018), 2MASS apparent K magnitude, the MESA Isochrones and Stellar Tracks models (MIST; Choi et al. 2016) along with priors from the preferred SpecMatch model to measure the best fitting solution of stellar properties.

Figure 3 .
Figure 3. Top: The S-Value and Velocity correlation.We found a Pearson correlation coefficient of 0.62 between the SHK values and RVs.Bottom Left: Radial Velocity periodogram, with a maximum peak at 227 days.Note the secondary peak at 787 days was investigated and ruled out as the solution to the fit.Bottom Right: The S-value periodogram, with a maximum peak at ∼ 1300 days and another peak at 234 days.We tested our model for the outer planet via the methods outlined in Section 4.2 and concluded that the RV signals were not caused by stellar activity.

Figure 4 .
Figure 4.The RV measurements of TOI-1386.Panel a: Our complete RV time series with our preferred model shown in red.Panel b: residuals between the RVs and best-fit 2 planet model.Panels c and e: Phase-folded RV time series for planets b and c respectively.Panels d and f: The residuals to each of the planet fits.

Figure 5 .
Figure 5. Eccentricity evolution of the TOI-1386 b and c planets.Shown here are the first 10 6 years from the full 10 7 year dynamical simulation, described in Section 4.4.

Figure 6 .
Figure 6.Heat map of the mass and radius distribution of exoplanets showing the position of TOI-1386 b (black X) compared to the full catalog of known exoplanets (gray).TOI-1386 b is positioned in the sub-Saturn valley between the cluster of Jupiter-sized planets and that of the Super Earths/Neptunes.

Figure 7 .
Figure 7.The attributes of planet TOI-1386 b (Panels A-D) and star TOI-1386 (Panels E & F) compared to the population of sub-Saturn planets.TOI-1386 b's position is indicated by the cross.Planets with a measured mass between 6 -60 M⊕ are included in our sub-Saturn population.The colorbar for each plot indicates the TSM value of the planet, with TOI-1386 b having a TSM of 42.5.Compared to the rest of the sub-Saturn population, TOI-1386 b is a massive, cool, moderately eccentric planet orbiting a metal rich sun-like star.edges support from a 51 Pegasi b Postdoctoral Fellowship from the Heising-Simons Foundation.A.B. is supported by the NSF Graduate Research Fellowship, grant No. DGE 1745301.R.A.R. is supported by the NSF Graduate Research Fellowship, grant No. DGE 1745301.C.D.D. acknowledges the support of the Hellman Family Faculty Fund, the Alfred P. Sloan Foundation, the David & Lucile Packard Foundation, and the National Aeronautics and Space Administration via the TESS Guest Investigator Program (80NSSC18K1583).J.M.A.M. is supported by the NSF Graduate Research Fellowship, grant No. DGE-1842400.J.M.A.M. also acknowledges the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant No. 1829740, the Brinson Foundation, and the Moore Foundation; his participation in the program has benefited this work.T.F.acknowledges support from the University of California President's Postdoctoral Fellowship Program.J.V.Z.acknowledges support from the Future Investigators in NASA Earth and Space Science and Technology (FINESST) grant 80NSSC22K1606.We thank the time assignment committees of the University of California, the California Institute of Technol-

Table 1 .
Radial Velocity Time Series