TOI-199 b: A Well-characterized 100 day Transiting Warm Giant Planet with TTVs Seen from Antarctica

We present the spectroscopic confirmation and precise mass measurement of the warm giant planet TOI-199 b. This planet was first identified in TESS photometry and confirmed using ground-based photometry from ASTEP in Antarctica including a full 6.5 hr long transit, PEST, Hazelwood, and LCO; space photometry from NEOSSat; and radial velocities (RVs) from FEROS, HARPS, CORALIE, and CHIRON. Orbiting a late G-type star, TOI-199 b has a 104.854−0.002+0.001day period, a mass of 0.17 ± 0.02 M J, and a radius of 0.810 ± 0.005 R J. It is the first warm exo-Saturn with a precisely determined mass and radius. The TESS and ASTEP transits show strong transit timing variations (TTVs), pointing to the existence of a second planet in the system. The joint analysis of the RVs and TTVs provides a unique solution for the nontransiting companion TOI-199 c, which has a period of 273.69−0.22+0.26days and an estimated mass of 0.28−0.01+0.02MJ . This period places it within the conservative habitable zone.


INTRODUCTION
In the quest to understand how planets form and evolve, warm giant planets are a vital piece of the puzzle.Warm giants are generally defined as planets with sizes R p ≥ 4R ⊕ and periods 10 d ≲ P ≲ 300 d.Unlike their closer-in hot Jupiter cousins, they are far enough from their host stars to not be strongly irradiated, and thus are not expected to be affected by radius inflation (see, e.g.Sarkis et al. 2020).Likewise, tidal circularisation is not expected to impact their orbital parameters, meaning these can potentially be used to disentangle their migration history (see Dawson & Johnson 2018, for a review).
Identifying and characterizing a large sample of warm giants is the necessary first step towards understanding this population.Transiting warm giants around bright stars are especially valuable, as both their radii (via transits) and their masses (via radial velocities) can be measured, and thus their density computed and their internal structure modelled.The Transiting Exoplanet Survey Satellite mission (TESS, Ricker et al. 2015) is proving invaluable for the detection of these warm giants, with some 60 already published (according to the table of TESS planets on the NASA Exo-planet Archive1 , as of 6 February 2023) and hundreds more expected from yield simulations of the prime and extended missions (Sullivan et al. 2015; Barclay et al. 2018;Kunimoto et al. 2022).Thanks primarily to TESS and its predecessor Kepler (Borucki et al. 2010), some 80 warm giants have had their radii and masses characterized to better than 25%, according to the TEPCAT catalogue (Southworth 2011) 2 .The majority of these are on the shorter end of the 10-300 d period range, with half having periods shorter than 20 days.Longerperiod transiting warm giants remain rare.The Warm gIaNts with tEss collaboration (WINE, e.g.Brahm et al. 2019;Jordán et al. 2020;Schlecker et al. 2020;Hobson et al. 2021;Trifonov et al. 2021) seeks to confirm and characterize warm giant candidates from TESS, using a variety of photometric and spectroscopic facilities.
While the transit technique has been immensely fruitful, leading the list of planet discoveries with over 3900 of the over 5200 known exoplanets 3 , it is of necessity limited by orbital geometry -the planet must cross be-tween the star and the observer -and biased towards short-period planets whose transits are more likely to be observed.However, for multi-planetary systems where the planets exert strong gravitational influences on each other, the transit timing variations (TTVs, Agol et al. 2005;Holman & Murray 2005) technique can help to overcome these limitations in some cases: by measuring the TTVs for a transiting planet, we can detect nontransiting planets, and measure the masses of all bodies involved.The first detection of the TTV effect was for Kepler-19 b by Ballard et al. (2011), but they could not unambiguously determine the mass and period of the non-transiting companion.The first planet with a unique solution determined from TTVs was KOI-872 c (Nesvorný et al. 2012).Since then, this technique has enabled the detection of 24 planets4 , many of which are warm giants (e.g.Trifonov et al. 2021, a pair of warm giants where only the inner one transits).
In this paper, we present the TOI-199 system, consisting of two giant planets.The inner planet is transiting, with a period of 104.854 +0.001  −0.002 d, and was found in the TESS photometry.The outer planet does not transit, but was revealed by the TTVs it induces on the transits of the inner planet.The paper is organized as follows.We present the observational data in Section 2. In Section 3 we analyse the data, and characterize the star and planet.Our results are discussed in Section 4, and finally, we provide our conclusions in Section 5.
The 2-minute cadence data for TOI-199 were processed by the TESS Science Processing Operation Center (SPOC, Jenkins et al. 2016) at NASA Ames Research Center.A potential transit signal was identified in the SPOC transit search (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2020) ) of the light curve, and designated as a TESS Object of Interest (TOI) by the TESS Science Office based on SPOC data validation results (Twicken et al. 2018;Li et al. 2019) indicating that it was consistent with a transiting planet.The SPOC difference image centroiding analysis locates the source of the transit signal to within 1.0 ± 2.5" of the target star's catalog position (Twicken et al. 2018).The planetary candidate , so named because it is the first transit signal detected in this light curve, is listed in the ExoFOP-TESS archive 5 , with a period of P 01 = 104.87075d.Prior to its release as a TOI, the WINE team had independently identified the initial single transit event in Sector 2 and begun follow-up.
Unfortunately, one of the transits, close to the start of sector 10, is cut off by quality checks in the 2-minute cadence Presearch Data Conditioning (PDCSAP, Stumpe et al. 2012Stumpe et al. , 2014;;Smith et al. 2012) light curve, and another at the end of sector 32 is distorted.However, both are clearly seen in the light curves obtained from the full frame images (FFIs), which we extracted with our own tesseract pipeline (Rojas et al. in prep. 6).We attempted to re-reduce the 2-minute cadence data to recover them, but were unable to retrieve the transit at the edge of sector 10, where the light curve drops off sharply.Therefore, we chose to use the FFI light curves from tesseract for the whole of our analysis.The PD-CSAP and FFI light curves for the seven sectors with transits are shown in Fig. 1; note that the cadence for the FFIs changed from 30 minutes for the prime mission (top row) to 10 minutes for the first extended mission (middle row), and to 200s for the second extended mission (bottom row).

Follow-up photometry
At 21 ′′ per pixel, the TESS pixels are relatively large, meaning that nearby companions can contaminate its photometry.Follow-up photometry is used to identify such cases of false positives, and to monitor transits not observed by TESS.TOI-199 was observed by five facilities, described in the following sub-sections.

ASTEP
5 Located at https://exofop.ipac.caltech.edu/tess/target.php?id= 309792357 Antarctica Search for Transiting ExoPlanets (ASTEP, Guillot et al. 2015;Mékarnia et al. 2016) is a 0.4 m telescope equipped with a Wynne Newtonian coma corrector, located on the East Antarctic Plateau.Until December 2021 it was equipped with a 4k × 4k frontilluminated FLI Proline KAF-16801E CCD with an image scale of 0."93 pixel −1 resulting in a 1 o × 1 o corrected field of view.The focal instrument dichroic plate split the beam into a blue wavelength channel for guiding, and a non-filtered red science channel roughly matching an Rc transmission curve (Abe et al. 2013).
In January 2022, the focal box was replaced with a new one with two high sensitivity cameras including an Andor iKon−L 936 at red wavelengths.The image scale is 1.39" pixel −1 with a transmission curve centred on 850 ±138 nm (Schmider et al. 2022).
The telescope is automated or remotely operated when needed.Due to the extremely low data transmis-sion rate at the Concordia Station, the data are processed on-site using IDL (Mékarnia et al. 2016) and Python (Dransfield et al. 2022) aperture photometry pipelines.
Three observations of TOI-199 were made with ASTEP, all of them with clear sky, winds between 3 and 5 m s −1 and temperature ranging between −62 and −69 o C.An egress was observed on 31st March 2021 with a FWHM of 4.3", corresponding to the transit observed by TESS in Sector 36.Full transits were observed on 13th July 2021 and 6th September 2022 with a FWHM of 4.3 and 5.9" respectively, which were not observed by TESS.The Moon, 80% illuminated, was present during the 31st March 2021 and 6th September 2022 observations.The light curves are shown in Fig. 2 (top row).

PEST
The Perth Exoplanet Survey Telescope (PEST) is a backyard observatory in Perth, Australia, operated by Thiam-Guan Tan.PEST was used to observe a transit egress of TOI-199.01 on 16th December 2021, corresponding to the transit observed by TESS in Sector 32.The observation was conducted with an Rc filter and 30s integration times.At the time, the 0.3 m telescope was equipped with a 1530 × 1020 SBIG ST-8XME camera with an image scale of 1.2 ′′ /pixel resulting in a 31 ′ × 21 ′ field of view.A custom pipeline based on C-Munipack was used to calibrate the images and extract the differential photometry using an aperture with radius 8.6 ′′ .The light curve is shown in Fig. 2 (middle row, left panel).

NEOSSat
The Near Earth Object Surveillance Satellite (NEOSSat, Laurin et al. 2008) is a Canadian microsatellite with a 0.15 m F/6 telescope which has a 0.86 • ×0.86 • field of view.Although designed for detecting and tracking near-Sun asteroids, it is also suitable for transit observations of bright stars (Fox & Wiegert 2022).NEOSSat observed a full transit of TOI-199.01 on 31st March 2021 with no filter, corresponding to the one observed by TESS in Sector 36.The light curve is shown in Fig. 2 (middle row, centre panel).

Hazelwood
The Hazelwood Observatory is a backyard observatory operated by Chris Stockdale in Victoria, Australia, with a 0.32 m Planewave CDK telescope working at f/8, a SBIG STT3200 2.2k × 1.5k CCD, giving a 20 ′ × 13 ′ field of view and 0.55 ′′ per pixel.The camera is equipped with B, V, Rc, Ic, g', r', i' and z' filters (Astrodon Interference).Typical FWHM is 2.2 ′′ to 3.0 ′′ .The Hazelwood Observatory observed an egress of TOI-199.01 in Rc on 31st March 2021, corresponding to the transit observed by TESS in Sector 36.Time series images were collected, bias, dark and flat fielded and the transit analysed with AstroImageJ (Collins et al. 2017).The reduced data was then uploaded to ExoFOP.The light curve is shown in Fig. 2 (middle row, right panel).

LCO
The Las Cumbres Observatory global telescope network (LCOGT, Brown et al. 2013) is a globally distributed network of telescopes.The 1 m telescopes are equipped with 4096 × 4096 SINISTRO cameras having an image scale of 0.389.′′ per pixel, resulting in a 26 ′ ×26 ′ field of view.LCO observed TOI-199 several times with the Sinistro instrument: a pre-ingress flat curve on 16th December 2021, corresponding to the transit observed by TESS in Sector 32, from the Cerro Tololo Inter-American Observatory (CTIO) site in zs filter; an ingress on 8th February 2022, from the South African Astronomical Observatory (SAAO) site, in alternating gp and ip filters; the egress of the same transit, from the CTIO site, also in alternating gp and ip filters; and a full transit on 20th Decemeber 2022, from the Siding Springs Observatory (SSO) site.The images were calibrated by the standard LCOGT BANZAI pipeline (Mc-Cully et al. 2018), and photometric data were extracted using AstroImageJ.The light curves are shown in Fig. 2 (bottom row).

High-resolution imaging
Given the relatively large pixel size of TESS, highresolution imaging from the ground is a valuable tool to assess possible blends and contamination from nearby sources.The SOAR TESS survey (Ziegler et al. 2020) observes TESS exoplanet candidate hosts with speckle imaging using the high-resolution camera (HRCam) imager on the 4.1-m Southern Astrophysical Research (SOAR) telescope at Cerro Pachón, Chile (Tokovinin 2018), in order to detect such nearby sources.TOI-199 was observed on the night of 18 February 2019, with no nearby sources detected within 3 ′′ .The contrast curve and speckle auto-correlation function are shown in Fig. 3.

Spectroscopic data
The WINE consortium carried out spectroscopic follow-up of TOI-199 with the HARPS, FEROS, CORALIE, and CHIRON spectrographs.We also received additional CORALIE and CHIRON RVs from other teams, which we incorporated into the analysis.All the resulting radial velocities (RVs) are shown in Fig. 12 (top panel).Figure 4 shows the GLS periodograms of the joint RVs from all four spectrographs (top panel); of the HARPS RVs alone (second panel); of the joint RV residuals to a fit of the inner planet alone (third panel); and of the H α and log(R ′ HK ) activity indicators computed from the HARPS spectra (fourth and fifth panels).

HARPS
We observed TOI-199 with the HARPS spectrograph (Mayor et al. 2003) at the 3.6m telescope at La Silla, resolving power R = 120 000, between 12th December 2018 and 9th March 2021.We obtained 46 spectra, under Program IDs 0101.C-0510, 0102.C-0451, 0104.C-0413, and 106.21ER.001.For the first three programs (December 2018 through February 2020), we operated in simultaneous sky mode, obtaining 26 spectra with a 900s exposure duration.For the last (December 2020 through March 2021), we switched to simultaneous Fabry-Perot mode as part of an overall revision of our observing program, obtaining 20 spectra with an increased 1200s exposure duration.The spectra have a median SNR of 62.We processed the spectra using the CERES pipeline (Brahm et al. 2017a), which provides both RVs and several activity indicators: the CCF bisector (BIS, e.g.Queloz et al. 2001a), and the H α (Boisse et al. 2009), log(R ′ HK ) (Duncan et al. 1991;Noyes et al. 1984), Na II, and He I (Gomes da Silva et al. 2011) activity indices.The resulting RVs and activity indicators are listed in Table A1; the RVs have a median error of 2.4 m/s.

FEROS
We observed TOI-199 with the FEROS spectrograph (Kaufer et al. 1999) at the MPG 2.2m telescope at La Silla, resolving power R = 50 000, between 27th November 2018 and 4th March 2020.We obtained 48 spectra, under Program IDs 0102.A-9003, 0102.A-9006, 0102.A-9011, 0102.A-9029, 0103.A-9008, and 0104.A-9007 in Object-Calibration mode, with a 500s exposure duration (save three spectra with increased exposure time due to poor observing conditions, two of 700s and one of 1000s).The spectra have a median SNR of 90.As with the HARPS data, we processed the spectra using the CERES pipeline.The resulting RVs and activity indicators are listed in Table A2; the RVs have a median error of 6.3 m/s.

CORALIE
TOI-199 was observed with the CORALIE spectrograph (Queloz et al. 2001b) at the Swiss 1.2m Euler telescope at La Silla by both the WINE and Swiss teams.In total, 19 spectra were obtained between 27th December 2018 and 2nd December 2019, with a median exposure duration of 1200s.CORALIE is a fibre-fed echelle spectrograph with a 2 ′′ science fibre, and a secondary fibre with a Fabry-Perot for simultaneous wavelength calibration.It has a spectral resolution of R = 60 000.RVs are extracted using the standard CORALIE DRS by cross-correlating the spectra with a binary G2V template (Baranne et al. 1996;Pepe et al. 2002).The BIS, FWHM, and other line-profile diagnostics were also computed via the CORALIE DRS, as was the H α index for each spectrum to check for possible variation in stellar activity.The resulting RVs and activity indicators are listed in Table A3; the RVs have a median error of 13.7 m/s.2.4.4.CHIRON TOI-199 was observed with the CHIRON spectrograph (Tokovinin et al. 2013) at the SMARTS 1.5-meter telescope at CTIO by both the WINE team and Dr Carleo.In total, 26 spectra were obtained between 11th March 2019 and 13th December 2020; the WINE team observed with single 1500s exposures, while Dr Carleo observed sets of three 600s exposures.RVs were extracted following Jones et al. (2019).The resulting RVs and activity indicators are listed in table A4; the RVs have a median error of 11.3 m/s.  ) activity indicator (bottom panel).The grey solid, dashed, and dotted horizontal lines indicate false alarm probability levels of 10%, 1% and 0.1% respectively.The periods of the two planets are indicated with vertical orange and green lines.We note that the peak at ∼81d is a one-year alias of the ∼104d planet.

Stellar parameters
We use a two-part process to characterise the host star.First, we employed the co-added FEROS spectra to determine the atmospheric parameters of  To do this, we used the ZASPE code (Brahm et al. 2017b), which compares the observed spectrum to a grid of synthetic models (generated from the ATLAS9 model atmospheres, Castelli & Kurucz 2003) in order to determine the effective temperature T eff , surface gravity log g, metallicity [Fe/H], and projected rotational velocity v sin i.
Second, we followed the second procedure described in Brahm et al. (2019) to determine the physical parameters of the host star.To summarize, we compare the broadband photometric measurements (converted to absolute magnitudes via the Gaia DR3 (Gaia Collaboration et al. 2016Collaboration et al. , 2022) ) parallax) with the stellar evolutionary models of Bressan et al. (2012).We use the emcee package (Foreman-Mackey et al. 2013) to sample the posterior distributions.Through this procedure, we determine the age, mass, radius, luminosity, density, and extinction.This also enables us to compute new values for T eff and log g, the latter of which is more precise than the one previously determined through ZASPE.Therefore, we iterate the entire procedure, using the new log g as an additional input parameter for ZASPE.One iteration was sufficient for the log g value to converge.
The final stellar parameters computed, and other observational properties of TOI-199 are presented in Table 1.We find that TOI-199 is a late G-type star.The quoted error bars for the stellar parameters are internal errors, computed as the ±1σ interval around the median of the posterior for each parameter; they do not take into account systematic errors in the stellar models.
Recently, Tayar et al. (2022) examined the uncertainty floors on fundamental stellar parameters.Specifically, they find fundamental floors of 2.4% for T eff , 2% for L ⋆ , and 4.2% for R ⋆ due to the uncertainties on input observables such as the bolometric flux and angular diameter.We also report these values in Table 1.For M ⋆ and the stellar age, they perform a comparison between four model grids.We use the online interface provided to carry out the same comparison for TOI-199, and report the standard deviation of the resulting values in Table 1.

TTV extraction
We used the juliet7 software (Espinoza et al. 2019) to extract the TTVs.We focussed on the TESS, ASTEP, and LCO transits, as the second and third  (Høg et al. 2000); PM13: using the tables of Pecaut & Mamajek (2013).Values in parentheses correspond to the fundamental uncertainty floors following Tayar et al. (2022).
ASTEP transits and second and third LCO transits are the only ones not also observed by TESS.The other light curves are consistent with the TESS observations, as shown in Figs. 8 and 9.The tesseract FFI light curves are uncorrected for systematics and other effects.We used a two-step method to fit a Gaussian Process (GP) to each TESS sector, first optimizing the hyperparameters on the outof transit data, with broad priors of J (10 −6 , 10 6 ) for σ GP,TESS,i , where J (a, b) is a Jeffreys or log-uniform distribution between a and b, and J (10 −3 , 10 3 ) for ρ GP,TESS,i for each TESS sector i.We then used the re-sulting posterior parameters as priors on the GP hyperparameters in the full fit including the transits, adopting the upper and lower 1σ bounds from the posteriors as limits on new Jeffreys distributions.Each transit is fit independently, with individual priors for the flux offsets m flux,instrument,transit and jitters σ instrument,transit .juliet also requires priors for the dilution factors m d,instrument,transit , which we have fixed to 1 for all instruments and transits.The limb-darkening parameters q 1,instrument and q 2,instrument are common to all transits observed with each instrument (save for LCO, where the observations were made at different sites or/and with different filters).For each q i,instrument we adopted truncated normal priors, where the mean of each prior was based on the derivation of limiting quadratic coefficients from the nonlinear coefficients described by Espinoza & Jordán (2015) 8 .We used ATLAS model atmospheres interpolated over 100 µ-points, as in (Claret & Bloemen 2011), and instrument response functions from the instrument website (for LCO, which provides profiles for all available filters) or from the SVO Filter Profile Service (Rodrigo et al. 2012;Rodrigo & Solano 2020; for TESS there is a specific profile; we adopted a generic Johnson V filter profile for NEOSSat, which observed without filter, and a Cousins R filter for the rest of the instruments, which observed with R filters.For the transit times T b,instrument,transit , we adopt uniform priors with a width of 6h, and midpoints determined from the orbital period listed in ExoFOP.The ExoFOP values for the planet-to-star radius ratio p b and the impact parameter b b were also used as priors.The eccentricity e b and ω b were fixed to 0 and 90 • , respectively.While we may expect planets at ∼ 100 d to have non-circular orbits, the shape of the transit will be noticeably altered only for large eccentricities, so a fixed eccentricity suffices for the TTV extraction (we note the eccentricity is a free parameter in the full RV+TTV model, see Section 3.3).The priors and posteriors for the full fit are listed in Table 2 for the planetary parameters, and in Table B1 for the instrumental and GP parameters.The fitted TESS, ASTEP, and LCO transits are shown in Figs. 5,  6, and 7 respectively, and the O-C plot is shown in Fig. 10.

Joint TTV and RV model
We used the Exo-Striker9 software (Trifonov 2019) to jointly fit the RV and the TTV data of TOI-199.We use a very similar approach to the TTV+RV N-body modelling scheme analysis done for the TOI-2202 system in Trifonov et al. (2021); we refer to this paper for more details on the TTV+RV modeling scheme used here.Briefly, the N-body RV modelling used is native to Exo-Striker, whereas the TTV modelling uses the TTVFast package (Deck et al. 2014).The fitted parameters for each planet in our joint model are the semiamplitude K, the orbital period P , eccentricity e, argument of periastron ω, and mean-anomaly M 0 .For our fitting analysis, we assumed coplanar, edge-on and prograde two-planet system (i.e., i b ,i c = 90 • and ∆i = 0 • ), and for the stellar mass of TOI-199 we adopt our best estimate of 0.936 +0.003 −0.005 M ⊙ (note that the stellar mass uncertainties are not included in the modelling, but are incorporated into the final uncertainties of the derived planetary parameters).The time step in the dynamical model was set to 1 day, assuring adequate model precision.Additionally, we optimise two parameters for each Doppler data set, the offset and the RV "jitter".This latter parameter is added in quadrature to the nominal uncertainties of the RV data (Baluev 2009).
We constructed parameter posteriors by the nested sampling algorithm (Skilling 2004), implemented via the dynesty package (Speagle 2020).For TOI-199 b, we used the previous light curve characterization with juliet and the RV period search results to constrain the priors on the parameters.For the period of the perturber, we first ran initial joint fits with broad priors in P c between 150 and 330 days, e c between 0 and 0.4, and ω c and M0 c between 0 • and 360 • covering first-and second-order mean-motion resonance (MMR) period ratios.From these, we could constrain the orbital period to a range of 265-290 d for the final fits.
We tested three models: two planets, two planets plus a linear trend, and two planets plus a quadratic trend.A visual hint of a long-term RV variation prompted the addition of trends.The comparison of the Bayesian likelihood factors favours the model with two planets and no trend (∆ ln Z = 15 between the no-trend and quadratictrend models).Therefore, we adopted the model with two planets and no additional trends, although we intend to continue monitoring this system for potential long-term RV variations that could point to the existence of more distant companions.The arguments of pericenter ω b and ω c required careful handling of the priors to avoid circular multimodality and ensure con- vergence within a reasonable time frame; we constrained them through preliminary fits.The final priors adopted are listed in Table 3.
The joint TTV and RV model results in TOI-199 b having a period of 104.854 +0.001 −0.002 d, an eccentricity of 0.09 +0.01 −0.02 , a radius of 0.888 ± 0.006 R J , and a mass of 0.17 Given the stellar radius, the planet's semi-major axis, and its predicted radii, this sets an upper limit on its inclination of i c ⪅ 89.7.Our model assumption that the system is edge-on and coplanar (i.e., i b , i c = 90 • , and ∆i = 0 • ) is still reasonable, since large mutual inclinations are unlikely, whereas small deviations from an edge-on configuration will not lead to a significant discrepancy from the derived (minimum) dynamical masses.
The phase-folded RVs for TOI-199 b and c are shown in Fig. 12 (bottom centre and right panels respectively), and the TTVs and fitted model in Fig. 12, bottom left panel.The full list of posterior parameters is given in Table 3, and the posterior probability distribution plots are shown in Figs.C1, C2, C3, and C4.The osculating orbital parameters, which evolve dynamically, are given for epoch BJD = 2458256.14.

Dynamical stability
To examine the dynamical stability of the TOI-199 system, we first considered the classical Hill stability criterion (Gladman 1993) and the Angular Momentum   Deficiency criterion (AMD, Laskar & Petit 2017).In terms of Hill stability, we found that the two planets are separated by ≈ 7.7R H , well above the ≈ 2.4R H limit given by Gladman (1993) (for initially circular orbits; however, the eccentricities of these planets are small).Likewise, the system is AMD-stable, with values of β = 0.3 and β s = 0.02 respectively (following Laskar & Petit 2017, Sect.5), well below the β = 1 limit for collisions between the planets or the innermost planet and the star, respectively.Therefore, we do expect the TOI-199 Saturn-mass pair to be stable at the estimated orbits.Nonetheless we aimed to analyse the dynamical evolution of the the best-fit, and study its properties.We performed a numerical N-body integration of the bestfit model with the Wisdom-Holman N-body algorithm (Wisdom & Holman 1991) for 10 Myr, adopting a small time step of dt =0.5 d.Fig. 13 shows the resulting evolution of the orbital semi-major axes, eccentricities, and ∆ω for an extent of 30 000 yr.We find the system to be well separated with no significant changes in the semi-major axes during the 10 Myr N-body simulation.The eccentricities, however, osculate with notable amplitudes of ≈ 0.06 for e b and ≈ 0.02 for e c .We find that the planetary pair osculates in aligned geometry, with the secular apsidal angle ∆ω = ω b − ω c librating around 0 deg, with a semi-amplitude of ∼ 20 deg.The libration of the angle between the periapses ∆ω is interesting, since it suggests that the planetary orbits must have been locked in secular apsidal alignment during the planet migration.Therefore, this libration of ∆ω is an important remnant evidence of the planets' formation and orbital evolution during the disk phase.

Possible exomoons in the system
The possibility of stable exomoons orbiting the planets of the TOI-199 system is interesting for two reasons.On the one hand, exomoons could trigger TTVs variations in TOI-199 b, which might be misinterpreted as a TTV signal caused by a second planetary companion (Kipping & Yahalomi 2022).The relatively small Hill-radius of TOI-199 b (R Hill ∼ 0.014 au), and the orbital analysis of the RV and TTV data we performed in Sec.3.3, strongly point to the existence of the outer Saturn-mass planet TOI-199 c, rather than an exomoon.Nonetheless, eliminating the exomoon hypothesis would strengthen the evidence for the existence of TOI-199 c even further.On the other hand, given the observational evidence of the existence of the Saturn-mass giant TOI-199 c, it is worth testing the possibility of exomoons around this companion.This planet orbits further out, so the Hill-radii of dynamical influence is somewhat larger (R Hill ∼ 0.03 au), positively affecting potential exomoons' stability.Further, TOI-199 c resides in the Habitable Zone (HZ), which is intriguing for the existence of stable exomoon bodies.The gas giant was likely formed even further out beyond the ice line around TOI-199.Thus, initially icy exomoons, similar to those of Saturn, are possibly in the current warm orbit, and could become potentially habitable, ocean-like, or early Mars-like worlds (see Trifonov et al. 2020, for further discussion).
We study the possibility of stable exomoons in the system following the same numerical setup as in Tri-fonov et al. (2020), who tested the exomoon stability of the Saturn-mass pair of planets around the M-dwarf GJ 1148.The N-body test was done using a Wisdom-Holman algorithm (Wisdom & Holman 1991) modified to handle the evolution of test particles in the Jacobi coordinate system (Lee & Peale 2003), which allows for invertibility of the system's hierarchy (i.e., making either of the planets a central body in the system, enabling the testing of stable orbits around each planet).Our Nbody simulations were performed with a very small time step of dt=0.01 d and for a maximum of 1 Myr, which we find to be a good balance between CPU resources and the numerical accuracy of the run.We randomly seed an arbitrary number of 4 000 "exomoon" mass-less test particles on circular prograde orbits around TOI-199 b & c.The semi-major axes of the test exomoons ranged randomly between the planetary Roche limit (∼ 0.0003 au, for both planets), and the planetary Hill-radius, where exomoons are expected to be stable.In addition, we studied the tidal interactions between the planets and the exomoons.We adopt Mars-like mass and radius for the exomoons (m = 0.107 M ⊕ , R = 0.53 R ⊕ ).In contrast, for the massive planets TOI-199 b & c, we adopt masses and radii from our best-fit analysis (the radius of TOI-199 c is unknown as it does not transit; thus, we assumed the same radius as for TOI-199 b).We integrated their orbits with the EqTide10 code (Barnes 2017), which calculates the tidal evolution of two bodies based on the models by Ferraz-Mello et al. (2008).For the rest of the numerical setup, we strictly follow the same EqTide prescription as in Trifonov et al. (2020).
The results from the test particle simulations are illustrated in Figure D1 for TOI-199 b, and Figure D2 for TOI-199 c, respectively.In the case of TOI-199 b, we found that only ∼ 48 % of the test exomoons remained long-term stable with small eccentricities, usually those near the stability limit, which we found to be at ∼ 0.0055 au (∼ 13 planet radii, ∼ 11.5 R J ).This stability limit is much smaller than the Hill-radius of the planet.As was shown in Trifonov et al. (2020) the expected stability limit for prograde exomoons is ∼ 0.5 R Hill , due to the dynamical perturbations of the second planet (see also, Grishin et al. 2017).Similarly, for TOI-199 c the stability region was found to be < 0.5 R Hill , where ∼ 60 % of the test particles survived for 1 M yr.As expected, the stable region around TOI-199 c is larger, close to 0.013 au (∼ 27 R J ) from the planet.The tidal interactions, however, further limit the possible region where the exomoons could exist around both planets.The Saturn- mass planets likely have comparable rotational periods to Saturn (∼ 10.5 h), which will not be affected significantly over the age of the system due to tidal interaction with the star.Therefore, closer Mars-like exomoons will drift away to longer orbits over time.Finally, we concluded that TOI-199 b is unlikely to host large exomoons because of the very limited range of semi-major axes of 0.0045-0.0055au where such bodies can exist.For TOI-199 c habitable exomoons could exist in the range 0.0045-0.0125au, comparable to the semi-major axes of the Galilean moons.

Interior models
We modelled the interior evolution of both planets in the system using CEPAM (Guillot & Morel 1995) and a non-grey atmosphere (Parmentier et al. 2015).We assumed simple structures consisting of a central dense core, composed of 50% ices and 50% rocks, surrounded by a hydrogen and helium envelope of solar composition.
We use the P-ρ relationships from Hubbard & Marley (1989) to calculate the density in the core.The equation of state used for hydrogen and helium accounts for nonideal mixing effects (Howard & Guillot 2023;Chabrier et al. 2019).
Figure 14 shows the resulting evolution models and the observational constraints.Defining an approximate bulk metallicity as Z ≈ M core /M tot and assuming Z ⊙ = 0.015, we find Z/Z ⊙ = 16 to 28 for TOI-199b.This is somewhat larger than what we obtain for a similarly simple model of Saturn although we stress that Saturn's enrichment corresponds to a narrower range of possibilities (between 12 and 15 according to Mankovich & Fuller 2021).However, most of the uncertainty comes from the poor constraint on the age of the system.An accurate determination of its age would thus greatly help to constrain the bulk metallicity of TOI-199 b (see also Müller & Helled 2023).We unfortunately do not have a measurement of the radius of TOI-199 c, as it does not Figure 13.N-body orbital evolution of the TOI-199 system for an extent of 30 000 yr.The initial orbit is taken from our TTV+RV best fit (see Table 3 and is   transit.Using the same evolution models and a composition with cores between 5 and 60 M ⊕ (bulk metallicity between 4 and 45), we obtain expected radii of TOI-199c between 0.6 and 1.03 R J .

Atmospheric characterization potential
TOI-199 b is a very interesting target for atmospheric characterization through transmission spectroscopy.We computed the transmission spectroscopy metric (TSM, Kempton et al. 2018) for TOI-199 b, finding a value of 107, above the TSM = 90 threshold recommended by Kempton et al. (2018).Compared to other well-characterised warm giants with similar masses and radii, it has a comparable TSM but is notably cooler (Fig. 15), making it a unique target.
The intrinsic luminosities of both TOI-199 b and c are between a third to one times the present-day luminosity of Jupiter, implying that their atmospheric structure is mostly governed by the irradiation that they receive (see Parmentier et al. 2015).With a photospheric temperature expected to range between 250 and 350 K, TOI-199 b is an ideal candidate for the observation of the consequences of condensation of water in giant planet atmospheres (see Guillot et al. 2022).

The TOI-199 system in a population context
The TOI-199 system is composed of two giant planets.The inner planet, TOI-199 b, is a transiting warm giant at a 104.854 +0.001 −0.002 d period with a 0.810 ± 0.005 R J radius, and a 0.17 ± 0.02 M J mass, comparable to those of Saturn (R S = 0.83 R J , M S = 0.30 M J ).To place TOI-199 b in context within the exoplanet populations, we plot the radius-period diagram (Fig. 16) for giant planets (R p ≥ 0.8 R J ) with well-constrained masses and radii from the TEPCAT catalogue (Southworth 2011) 11 .TOI-199 b helps populate a so far very sparse region in the period-radius space.It is noticeably less massive, less dense and less eccentric than its closest neighbours, and is the first warm exo-Saturn (i.e., a warm giant with mass and radius similar to those of Saturn) with mass measured to better than 20% and radius measured to better than 10%.
The outer planet, TOI-199 c, also falls within the warm giant parameter space at a period of 273.69 +0.26  −0.22 d.With a minimum mass of m sin i = 0.28 +0.02 −0.01 M J , it is significantly more massive than the inner planet; as it does not transit, we cannot measure its radius.Given  its orbital period, TOI-199 c falls into the conservative HZ, as defined by Kopparapu et al. (2014).The orbits of the two planets, and the conservative and optimistic HZ regions for TOI-199, are shown in Fig. 17.However, given its minimum mass, TOI-199 c is presumably a gas giant, and as such does not possess a surface on which water could be found.We use the stellar effective temperatures and luminosities listed in the NASA Exoplanet Archive 12 to compute the Habitable Zones defined by Kopparapu et al. (2014) for all stars with temperatures in the 2600 − 7200 K range and catalogued luminosities.While most exoplanet host stars (4092 out of the 5035 catalogued) fall in this temperature range, relatively few of those (897) have catalogued luminosities.In these, we find 159 planets in the conservative or optimistic Habitable Zones of their host stars, of which 138 have measured masses.Out of the planets with measured masses, the vast majority (111) have M p > 10M ⊕ , placing them in the Neptune-or-larger regime.True potentially habitable planets, with surfaces on which water could exist in a liquid state, are still rare.However, these giant planets, among which TOI-199 c falls, could still host potentially habitable exomoons.

CONCLUSIONS
We have presented the discovery and characterization of the TOI-199 system, composed of two warm giant planets.The inner planet, TOI-199 b, was first identified as a single transiter in TESS photometry, and confirmed using ASTEP, LCO, PEST, Hazelwood, and NEOSSat follow-up photometry, and HARPS, FEROS, CORALIE, and CHIRON RVs.The transits of TOI-199 b showed strong TTVs, pointing to the existence of a second planet in the system.From the joint analysis of the TTVs and RVs, we found that TOI-199 b has a 104.854 +0.001 −0.002 d period, a mass of 0.17 ± 0.02 M J , and a radius of 0.810 ± 0.005 R J , making it the first precisely characterized warm exo-Saturn.Meanwhile, the outer non-transiting planet TOI-199 c has a period of 273.69 +0.26  −0.22 d, placing it in the Habitable Zone, and a minimum mass of 0.28 +0.02 −0.01 M J .We studied the dynamical stability and the potential for these planets to host exomoons.Our N-body simulations show that the system is stable, with no significant changes in the semi-major axes.The secular apsidal angle ∆ω, however, librates, suggesting the planets' orbits were locked in secular apsidal alignment during their migration.TOI-199 b is extremely unlikely to host large exomoons.TOI-199 c, on the other hand, could potentially host habitable exomoons, though as it does not transit these would be exceedingly difficult to detect.In this appendix, we show the RVs and activity indicators (where applicable) obtained from the spectroscopical data.Tables A1 and A2 show the HARPS and FEROS results respectively, both obtained from the ceres pipeline.Table A3 shows the CORALIE results, obtained using the standard CORALIE DRS.Table A4 shows the CHIRON results, obtained following Jones et al. (2019).

C. POSTERIOR PROBABILITY DISTRIBUTIONS
In this appendix, we show the posterior probability distribution of the joint TTV+RV modelling with Exo-Striker.For ease of viewing, we have separated the cornerplots into fitted planetary parameters (Figs.C1 and C2), derived planetary parameters (Fig. C3) and instrumental RV parameters (Fig. C4).

D. EXOMOON DYNAMICAL EVOLUTION
In this appendix, we show the dynamical evolution of the test particles used in the exomoon analysis described in Sec.4.1.

Figure 1 .
Figure 1.TESS light curves for TOI-199, for the seven sectors with transits.The short-cadence PDCSAP light curves are shown in blue, and the FFI light curves extracted with tesseract in orange.For Sectors 10 and 32 the transits could not be correctly recovered from the PDCSAP light curves.

Figure 2 .
Figure 2. Follow-up photometry for TOI-199.Top row: Light curves for ASTEP observations -an egress on 31st March 2021 (left), a full transit on 13th July 2021 (centre), and a full transit on 6th September 2022 (right).Middle Row: Light curves for PEST (left, egress on 16th December 2020), NEOSSat (centre, full transit on 31st March 2021), and Hazelwood (right, egress on 31st March 2021) observations.Bottom row: light curves for LCO observations -a pre-ingress flat curve on 16th December 2020 (left), two ingresses and egresses on 8th February 2022 (centre), for which the points are colour-coded by filter and site, and a full transit on 20th December 2022 (right).In all panels, the dashed vertical line indicates the transit midpoint, and the dotted vertical lines the egress and ingress.

Figure 3 .
Figure 3. Contrast curve and speckle auto-correlation function from the HRCam at SOAR for TOI-199.The black points and solid line indicate the 5σ contrast curve; the inset shows the speckle auto-correlation function.No nearby sources are detected.

Figure 4 .
Figure 4. Periodograms of: the joint RVs (top panel); the HARPS RVs alone (second panel); the joint RV residuals to a one-planet fit for TOI-199 b (third panel); the HARPS Hα indicator (fourth panel); and the HARPS log(R ′HK ) activity indicator (bottom panel).The grey solid, dashed, and dotted horizontal lines indicate false alarm probability levels of 10%, 1% and 0.1% respectively.The periods of the two planets are indicated with vertical orange and green lines.We note that the peak at ∼81d is a one-year alias of the ∼104d planet.
T0 b b [BJD] -2458256.1286± 0.0006 R b . . . .[RJ] -0.810 ± 0.005 i b . . . .[deg] -89.81 ± 0.02 a U(a, b) indicates a uniform distribution between a and b; N (a, b) a normal distribution with mean a and standard deviation b; J (a, b) a Jeffreys or log-uniform distribution between a and b. b juliet computes the P and T0 in a TTV fit as the slope and intercept, respectively, of a least-squares fit to the transit times.We adopt the value from the full TTV+RV model (Sect.3.3) as our final period.

Figure 5 .
Figure 5. TESS FFI light curves, extracted with tesseract (blue points), and fitted models (black lines) for the seven transits observed with TESS.The GP components of each model have been subtracted.
± 0.02 M J .Meanwhile, TOI-199 c has a period of 273.69 +0.26 −0.22 d, an eccentricity of 0.096 +0.008 −0.009 , and a minimum dynamical mass (since in the dynamical model we fix i c = 90 • ) of m c = 0.28 +0.02 −0.01 M J .If it were a transiting planet, TESS should have seen transits in Sectors 5 and 35, but the light curves are flat around the predicted times of transit.In Fig. 11 we show the light curves together with model transits for the expected range of radii from interior models (Sec.4.2) and impact parameters of b = 0 and b = 0.75, any of which should have been easily detectable in the TESS light curves.

Figure 6 .
Figure 6.ASTEP light curves (blue points) and fitted models (black lines) for the three transits observed with ASTEP.

Figure 7 .
Figure 7. LCO light curves (coloured points) and fitted models (black lines) for the three transits observed with LCO.The data for transit 13 are coloured by filter and site, as in Fig. 2.

Figure 8 .
Figure 8. TESS FFI light curve for sector 32 extracted with tesseract (blue points), PEST light curve (orange points), LCO light curve (green points) and fitted model to the TESS data (black line) for transit 9.

Figure 9 .
Figure 9. TESS FFI light curve for sector 36 extracted with tesseract (blue points), ASTEP light curve (orange points), NEOSSat light curve (green points), Hazelwood light curve (red points), and fitted model to the TESS data (black line) for transit 10.

Figure 10 .
Figure 10.O-C plot (difference between predicted and observed transit midpoint times) for the TESS, ASTEP, and LCO transits.The error bars are generally smaller than the marker size.

Figure 11 .
Figure 11.Predicted transits for TOI-199 c, which does not transit.Blue and orange points indicate the PDC-SAP and tesseract light curves respectively, as in Fig. 1.The vertical dashed line indicates the expected time of midtransit.Model transits from batman are shown for two radii, Rc = 1.03 RJ (solid lines) and Rc = 0.6 RJ (dotted lines), and for two impact parameters, b = 0 (green) and b = 0.75 (purple).

Figure 12 .
Figure 12.Top: Radial velocity measurements (FEROS: blue, HARPS: red, CHIRON: purple, CORALIE: green) and best-fit two-planet model (grey line) for the combined data.Bottom left: TTVs (blue circles) and fitted model with Exo-Striker (grey line) for TOI-199 b (top panel) and residuals to the model (bottom panel).The model has been smoothed with a quadratic spline.Bottom centre and right: phase-folded representation of the two planetary signals after the RV signal of the other companion was subtracted.The respective RV residuals are shown under each panel, accordingly.The more precise HARPS RVs, which are the main driver of the fit, are highlighted.The black squares show the binned phase-folded RVs.
Figure13.N-body orbital evolution of the TOI-199 system for an extent of 30 000 yr.The initial orbit is taken from our TTV+RV best fit (see Table3and is stable for at least 10 Myr.Left to right: evolution of the semi-major axes a b (blue) and ac (red), the eccentricities e b (blue) and ec (red), and the evolution of ∆ω = ω b -ωc, respectively.See text for details.

Figure 14 .
Figure 14.Evolution models of TOI-199 b and TOI-199 c compared to a simple model for Saturn.All models assume a central ice-rock core surrounded by a hydrogen and helium envelope of solar composition.Mcore corresponds to the mass of the ice/rock core while Menv corresponds to the mass of the solar-composition envelope.The range of Mcore and Menv compatible with the observational constraints is shown for TOI-199 b.The blue error bar corresponds to observational constraints on the age and the radius of TOI-199 b.We also show the range of radii expected for likely extreme compositions of TOI-199 c in red.

Figure 15 .
Figure15.Mass-radius diagram for warm (Teq < 1000 K) giant planets with masses measured to better than 20% and radii measured to better than 10%, similar to our accuracies for TOI-199 b.The markers are colour-coded by equilibrium temperature and scaled by transmission spectroscopy metric.TOI-199 b is significantly cooler than other well-characterized planets with similar masses and radii, and has a high TSM.

Figure 16 .
Figure16.Radius-period diagram for giant planets with masses measured to better than 20% and radii measured to better than 10%, similar to our accuracies for TOI-199 b.The markers are colour-coded by eccentricity and scaled by planet mass.TOI-199 occupies a scarcely-populated region of this space.

Figure 17 .
Figure 17.Orbital configuration of the TOI-199 system.The conservative and optimistic Habitable Zones are shown in dark and light green respectively.TOI-199 c is in the Habitable Zone.

a
U(a, b) indicates a uniform distribution between a and b; T N (a, b) a normal distribution with mean a and standard deviation b N (a, b, c, d) a normal distribution with mean a, standard deviation b, limited to the [c, d] range; J (a, b) a Jeffreys or loguniform distribution between a and b.

Figure C1 .
Figure C1.Cornerplot of the posterior distributions of the fitted planet parameters for TOI-199 b for the joint modelling of the TESS and ASTEP TTVs, and the HARPS, FEROS, CHIRON, and CORALIE RVs.The distributions are explored using nested sampling.The red crosses indicate the median values, and the black contour lines the 1σ, 2σ, and 3σ confidence levels.For ω b , note that it is a circular argument and thus, e.g., −10 • = 350 • .

Figure C2 .
Figure C2.Cornerplot of the posterior distributions of the planet parameters for TOI-199 c the joint modelling of the TESS and ASTEP TTVs, and the HARPS, FEROS, CHIRON, and CORALIE RVs.The distributions are explored using nested sampling.The red crosses indicate the median values, and the black contour lines the 1σ, 2σ, and 3σ confidence levels.For ωc, note that it is a circular argument and thus, e.g., 370 • = 10 • .

Figure C3 .
Figure C3.Cornerplot of the posterior distributions of the derived planet parameters for the joint modelling of the TESS and ASTEP TTVs, and the HARPS, FEROS, CHIRON, and CORALIE RVs.The distributions are explored using nested sampling.The red crosses indicate the median values, and the black contour lines the 1σ, 2σ, and 3σ confidence levels.

Figure C4 .
Figure C4.Cornerplot of the posterior distributions of the instrumental RV parameters for the joint modelling of the TESS and ASTEP TTVs, and the HARPS, FEROS, CHIRON, and CORALIE RVs.The distributions are explored using nested sampling.The red crosses indicate the median values, and the black contour lines the 1σ, 2σ, and 3σ confidence levels.Subindex H refers to HARPS, F to FEROS, CH to CHIRON, and CO to CORALIE.

Figure D1 .
Figure D1.Evolution of mass-less test particles (i.e., exomoons) around TOI-199 b, under the gravitational perturbation of the outer planet.Shown are the position of TOI-199 b (blue dot), the planetary Roche limit (gray dashed line), and the planetary 0.5R Hill (blue dashed line), which scales with (1-e b ) due to the dynamical perturbations of TOI-199 c.

Figure D2 .
Figure D2.Same as in Figure D1, but for TOI-199 c.

Table 2 .
Prior and posterior planetary parameter distributions for the TTV extraction with juliet.Top: Fitted parameters.Bottom: derived orbital parameters and physical parameters.

Table A1 .
RV and activity indices obtained from the HARPS spectra.

Table A2 .
RV and activity indices obtained from the FEROS spectra.
B. INSTRUMENTAL PRIORS AND POSTERIORS FOR THE TTV EXTRACTIONIn this appendix, we list the instrumental and GP prior and posterior distributions for the TTV extraction with juliet.

Table A3 .
RV and activity indices obtained from the CORALIE spectra.

Table B1 .
Prior and posterior instrumental parameter distributions for the TTV extraction with juliet.