TOI-4010: A System of Three Large Short-Period Planets With a Massive Long-Period Companion

We report the confirmation of three exoplanets transiting TOI-4010 (TIC-352682207), a metal-rich K dwarf observed by TESS in Sectors 24, 25, 52, and 58. We confirm these planets with HARPS-N radial velocity observations and measure their masses with 8 - 12% precision. TOI-4010 b is a sub-Neptune ($P = 1.3$ days, $R_{p} = 3.02_{-0.08}^{+0.08}~R_{\oplus}$, $M_{p} = 11.00_{-1.27}^{+1.29}~M_{\oplus}$) in the hot Neptune desert, and is one of the few such planets with known companions. Meanwhile, TOI-4010 c ($P = 5.4$ days, $R_{p} = 5.93_{-0.12}^{+0.11}~R_{\oplus}$, $M_{p} = 20.31_{-2.11}^{+2.13}~M_{\oplus}$) and TOI-4010 d ($P = 14.7$ days, $R_{p} = 6.18_{-0.14}^{+0.15}~R_{\oplus}$, $M_{p} = 38.15_{-3.22}^{+3.27}~M_{\oplus}$) are similarly-sized sub-Saturns on short-period orbits. Radial velocity observations also reveal a super-Jupiter-mass companion called TOI-4010 e in a long-period, eccentric orbit ($P \sim 762$ days and $e \sim 0.26$ based on available observations). TOI-4010 is one of the few systems with multiple short-period sub-Saturns to be discovered so far.


INTRODUCTION
Two major classes of planetary system architectures have arisen from the census of exoplanets close to their stars: hot Jupiters (R p > 8 R ⊕ , P < 10 days), and Neptune-sized and smaller planets (R p < 4 R ⊕ ) in compact multi-planet systems ("compact multis"). Hot Jupiters tend to have no co-planar transiting companions (e.g. Steffen et al. 2012;Hord et al. 2021) with few exceptions (e.g. WASP-47, Becker et al. 2015;TOI-1130, Huang et al. 2020a, have a wide range of orbital eccentricities, and are more common around metal-rich stars (Fischer & Valenti 2005). Meanwhile, compact multis tend to involve small planets in nearly circular and co-planar orbits (e.g. Xie et al. 2016;Van Eylen et al. 2019), and typically lack massive companions in similarly packed orbits.
Super-Neptune to sub-Saturn-sized planets (R p = 4 − 8 R ⊕ ) do not fit cleanly into either of the above classes. The most massive sub-Saturns (up to ∼ 80 M ⊕ ) are similar to hot Jupiters in that they typically orbit metal-rich stars and are the only detected planets in their systems, while the least massive sub-Saturn-sized planets (down to ∼ 5 M ⊕ ) tend to exist in low-eccentricity multiplanet systems ). However, very rarely do such systems feature sub-Saturns close to their stars. Only 63 sub-Saturns with P < 20 days in multiplanet systems have been confirmed so far, compared to 1213 smaller planets in these orbits. 1 Furthermore, only five known systems feature multiple sub-Saturns with P < 20 days: Kepler-18 (Cochran et al. 2011), Kepler-223 (Mills et al. 2016, K2-19 (Armstrong et al. 2015), V1298 Tau , and TOI-942 (Carleo et al. 2021). V1298 Tau and TOI-942 are both young stars (∼ 20 − 50 Myr old), meaning their planets could be precursors to the more typical compact systems of small exoplanets found by Kepler.
In this work, we present the confirmation of TOI-4010: a sixth example of a multi-planet system with two shortperiod sub-Saturns, in addition to an inner planet in the hot Neptune desert. We discuss the photometry, highresolution imaging, and spectroscopy observations that resulted in the discovery and confirmation of the system in §2. In §3 and §4, we describe our determination of the stellar and planet parameters, as well as the detection of a fourth planet in precise radial velocity (RV) observations. In §5 we place TOI-4010 in the context of other multi-planet systems and discuss the implications of our results.

Transiting Exoplanet Survey Satellite
NASA's Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) launched in April 2018 with the goal of searching for transiting exoplanets around nearby, bright stars. While TESS targets ∼ 20, 000 high-interest stars every 27.4-day sector with 2-minute cadence observations, it also enables the detection of exoplanets around millions more stars by regularly recording images of its entire field of view. These Full-Frame Images (FFIs) were saved every 30 minutes for the TESS Prime Mission (2018 July -2020 July), every 10 minutes for the first Extended Mission (2020 July -2022 Septem-ber), and every 200 seconds for the second Extended Mission (2022 September -).
However, at T = 11.5 mag, the host star was not bright enough for the signal to be further processed by the standard QLP vetting process, which is reserved for stars brighter than T = 10.5 mag (Guerrero et al. 2021). The signal was later vetted as part of an independent review of QLP search results for planets around faint stars (Kunimoto et al. 2022), which also uncovered a similarly sized candidate with R p = 6.4 ± 0.6 R ⊕ , P = 14.687 ± 0.003 days, T 0 = 1985.976 ± 0.003 TBJD, and S/N = 17. These signals were alerted on the MIT TESS Alerts portal 2 on 2021 June 23 as TOI-4010.01 and TOI-4010.02, respectively.
A third, smaller candidate with R p = 3.3 ± 0.6 R ⊕ , P = 1.3480 ± 0.0004 days, T 0 = 2007.533 ± 0.009 TBJD, and S/N = 8 was later detected with both an independent BLS search using a custom FFI light curve extraction and a multi-planet search from QLP. This signal was alerted as TOI-4010.03 on 2021 July 06.
The target was also observed in Sectors 52 (Camera 4, CCD 3, 2022 May 18 -June 13) and 58 (Camera 2, CCD 1, 2022 October 29 -November 26) in both the FFIs and at 2-minute cadence. 2-minute cadence data products were processed by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) at the NASA Ames Research Center. All three transiting planets were re-detected by QLP and SPOC in Sector 52, and passed all diagnostic criteria in the SPOC Data Validation (DV) reports (Twicken et al. 2018;Li et al. 2019) including even-odd transits comparison, eclipsing binary discrimination tests, tests against effects associ-2 https://tev.mit.edu/data/ ated with scattered light, and tests against background eclipsing binaries.

This Work
Here, we use custom 30-minute FFI light curves from Sectors 24 and 25 and the SPOC 2-minute cadence light curves from Sectors 52 and 58. We chose SPOC data for these sectors in order to take advantage of the faster cadence, which helps better resolve transit shapes. We followed the procedure of Vanderburg et al. (2019) to produce FFI light curves. In summary, after downloading the pixels surrounding TOI-4010 using TESScut (Brasseur et al. 2019), we extracted light curves from a series of both circular and irregularly shaped photometric apertures (Vanderburg et al. 2016). We then removed systematics by decorrelating each light curve with the mean and standard deviations of TESS spacecraft quaternion time series. Decorrelation was performed via linear regression with iterative outlier exclusion. We simultaneously modeled stellar variability by including a basis spline in our linear regression with breakpoints every 1.5 days. In-transit datapoints were manually masked out from the detrending process, and an extrapolation of the trends was used to estimate corrections during transit. We selected the final FFI light curve from the aperture that minimized the photometric scatter of out-of-transit data.
From Sectors 52 and 58, we used the SPOC Pre-search Data Conditioning Simple Aperture Photometry (PD-CSAP; Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014 light curve with all data points with a nonzero quality flag removed. We further removed low-frequency trends from the PDCSAP light curve using the biweight time-windowed slider implemented in the wotan Python package (Hippke et al. 2019) with a window of 0.5 days, and masked all in-transit points to avoid distorting the transit shapes.
The final four-sector light curve contained 2391 flux measurements at 30 minute cadence and 35,293 measurements at 2 minute cadence, with a total observing baseline of 954 days. We measured the 1 hour precision of the light curve to be 513 -583 ppm depending on the sector, estimated as 1.48 times the median absolute deviation (MAD) of the light curve after binning to 1 hour bins.
While TESS light curves are the primary data product for planet search and vetting, analysis of the TESS pixels can also help assess planet candidacy. Because photometers measure all light from pixels within a specific aperture, the flux from a nearby eclipsing binary (NEB) can contaminate the aperture and cause transitlike events in the light curve. These cases can be iden-tified using difference imaging, whereby the difference of averaged in-and out-of-transit pixel images is found (Bryson et al. 2013). The location of the strongest flux variation in the difference image should correspond to the transit source.
To produce difference images, we used TESS-plots 3 , a Python package for difference image generation from TESS FFIs.
For planets in multi-planet systems, TESS-plots masks out all cadences corresponding to other planet transits. This ensures that the primary source of variability in the difference image is due to the planet transit of interest. We inspected both Sector 24 and 25 images, with Sector 25 images shown in Figure  1. All difference images indicate that the three transit sources are co-located with the target star. These results are consistent with the difference images and centroid analysis from the Sector 52 SPOC DV reports, which find that the transit sources for all three planets are well within 3σ of the expected target location from the TESS Input Catalog (TICv8.2; Paegert et al. 2021): the transit source for TOI-4010.01 is located within 1.7 ± 2.6 ′′ of the target star, that for TOI-4010.02 is within 4.6±2.6 ′′ , and TOI-4010.03 is within 1.9 ± 2.8 ′′ .

Archival Images
TOI-4010 was observed as part of the second Palomar Sky Survey (POSS-II; Reid et al. 1991) between 1989 and 1995 in blue, red, and infrared bands. We obtained images of the corresponding finder charts from the NASA/IPAC Infrared Science Archive, 4 as shown in Figure 2. The star was also observed in the first survey by the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS1; Chambers et al. 2016). We obtained processed images in g, r, i, z, and y from the Pan-STARRS1 image cutout service. 5 Only the y band cutout from 2011 is included in Figure 2 due to saturation in other bands.
TOI-4010 has moved only ∼ 3 ′′ since the oldest image in 1989, and is relatively isolated. All six Gaia sources within 21 ′′ (the size of one TESS pixel) are fainter by ∆T ≈ 7 mag or more, and are therefore too faint to reproduce the transits.

Ground-Based Photometry
We obtained ground-based photometric observations through the TESS Follow-up Observing Program (TFOP) SG1, which is specialized in seeing-limited photometry to aid in the validation of TESS planet can-didates. Ground-based photometry can identify false positives due to NEBs at smaller separations than is achievable with difference image analysis. Observations can also provide improved transit ephemerides, independent transit depth constraints, and confirmation that the event is on-target. Our observations are summarized in Table 1, and all time series are available on the Exoplanet Follow-up Observing Program website 6 (Exo-FOP; NExScI 2022). At least two on-target full transits of each planet were detected.

Las Cumbres Observatory
Transits of all three planets were observed with telescopes in the Las Cumbres Observatory global telescope network (LCOGT; Brown et al. 2013). Observations taken by the 1m telescopes at the McDonald Observatory (LCO-McD) and Teide Observatory (LCO-Teide) use Sinistro cameras with a pixel scale of 0.389 ′′ px −1 . Observations taken by the 2m telescope at the Haleakala Observatory (LCO-HAL) use a MuS-CAT3 camera (Narita et al. 2020) with a pixel scale of 0.265 ′′ px −1 . Images were calibrated using the standard LCOGT BANZAI pipeline (McCully et al. 2018), and photometric data were extracted with AstroImageJ .
TOI-4010.01: A full transit of TOI-4010.01 was observed on 2021 August 28 by LCO-McD in gp and ip bands, each giving depths and durations consistent with the TESS transits. A 10.0 px (3.90 ′′ ) aperture optimized the photometric precision in both bands, and a smaller 5.0 px (1.95 ′′ ) aperture was used to clear the field of NEBs at the transit ephemerides out to 2.5 ′ .
TOI-4010.02: A full transit of TOI-4010.02 was observed on 2021 August 30 by LCO-Teide in gp and ip bands. A 15.0 px (5.85 ′′ ) aperture optimized the photometric precision, and a 5.0 px aperture was used to clear the field out to 2.5 ′ . A second full transit was observed on 2021 October 27, optimized by an 11.0 px (4.29 ′′ ) aperture. Both transits had depths and durations consistent with the TESS transits.
TOI-4010.03: A first attempt to detect TOI-4010.03 was made on 2021 September 8 at LCO-McD. No transit on-target was observed, which was later attributed to the TESS ephemerides predicting the transit too early. Observations on 2021 October 03 at LCO-Teide tentatively detected the event in a 9.0 px (3.51 ′′ ) aperture, but clarity depended on the choice of aperture size. TOI-4010.03 was also tentatively detected on October 08 at LCO-Teide and October 16 at LCO-McD in 14.0 px (5.46 ′′ ) apertures, but a smaller 10.0 px aperture did  , and infrared (center-right) filters from the second Palomar Sky Survey, and the y band (right) from the first Pan-STARRS survey. The location of TOI-4010 at the start of TESS Sector 24 in 2020 is marked with an orange cross. TOI-4010 has moved only ∼ 3 ′′ since the oldest image in 1989, and is relatively isolated, with no nearby brighter stars. not produce valid light curves. A further attempt on 2021 October 28 at LCO-Teide in the ip band was inconclusive.
The star was later re-observed by TESS in Sector 52, allowing us to precisely constrain the ephemerides and revisit the LCO data. We determined that partial transits of TOI-4010.03 occurred in the LCO observations from October 03 and 08, and an egress was caught shortly before the transit of TOI-4010.02 on August 30.
Post-TESS Sector 52, further attempts to detect TOI-4010.03 were made at LCO-HAL on 2022 July 11 and August 11 in gp, rp, ip, and zs bands. Both attempts were affected by poor weather, though the latter recovered on-time transits in each band. Three on-time full transits were finally detected by LCO-Teide on 2022 August 06, August 24, and September 06, providing definitive ground-based detections of TOI-4010.03.

Acton Sky Portal Observatory
A full transit of TOI-4010.01 was observed on 2021 September 08 at the private observatory Acton Sky Portal (ASP), using an SBIG ST-8-XME camera mounted on a 0.36m Celestron EdgeHD 1400. The camera has a pixel scale of 0.69 ′′ px −1 , and observations were taken with an R c filter. Cloud disruption affected data near ingress, and there is a systematic noise feature in the light curve due to an image reset. However, the transit depth and egress time matched predictions from

Fred Lawrence Whipple Observatory
A detection of TOI-4010.03 in the ip band was attempted with the 1.2m telescope at the Fred Lawrence Whipple Observatory (FLWO) on 2021 September 12. The telescope uses the KeplerCam, which has a pixel scale of 0.672 ′′ px −1 . We subsequently found that the observations were too early to cover the transit. Nevertheless, these were useful in determining that future observations should attempt to detect TOI-4010.03 later than predicted by the TESS data.

Observatori Astronòmic Albanyà
A full transit of TOI-4010.01 was observed on 2021 December 30 at the Observatori Astronòmic Albanyà using a Moravian G4-9000 camera mounted on the 0.4m OAA telescope. The camera has a pixel scale of 1.44 ′′ px −1 and observations were taken with an I c filter. Photometric data were extracted with AstroImageJ. The transit was detected on-target at the expected ephemeris, with the expected transit depth and duration, in an uncontaminated 10.0 px (14.4 ′′ ) aperture.

High-Resolution Imaging
We obtained high-resolution imaging to detect nearby objects that are not resolved in the TIC or by seeinglimited photometry. Unresolved eclipsing binaries may be the source of a false positive signal, while light from nearby stars can also dilute the transit depth, thus reducing the inferred planet radius.

W. M. Keck Observatory
TOI-4010 was observed on 2021 July 19 by the 10m Keck II telescope at the W. M. Keck Observatory. We obtained adaptive optics (AO) imaging in the K p band with the Near-Infrared Camera (NIRC2), 7 as well as data for non-redundant aperture masking interferometry (NRM; Tuthill et al. 2000). The image processing and analysis was conducting following the methods described by Kraus et al. (2016). Contrast curves from both types of observations are included in Figure 3. The AO images revealed a likely neighbor at a separation of 3.174 ± 0.004 ′′ , a position angle of 73.09 ± 0.07 • , and a contrast of ∆K p = 8.58 ± 0.07 mag, placing it roughly 0.6 mag above the detection limits. This neighbor is too faint to reproduce the transits of TOI-4010, and thus we do not consider it a false positive source or significant contaminant to the light curve. No other neighbors were revealed in either the AO or NRM analysis.
We obtained additional NIRC2 AO contrast curves through TFOP SG3 based on observations from 2021 August 28, shown in Figure 3, demonstrating that there are no nearby stars or companions in the field of view down to the instrumental detection limits. Observations were made in the K s band with an estimated PSF of 0.05 ′′ . These results are available on ExoFOP.

Caucasian Observatory
TOI-4010 was also observed on 2021 September 09 with the speckle polarimeter on the 2.5m telescope at the Caucasian Observatory of Sternberg Astronomical Institute (SAI) of Lomonosov Moscow State University (Safonov et al. 2017). Observations were made in the I band with an estimated PSF of 0.083 ′′ and obtained through TFOP SG3. The sensitivity curve is shown in Figure 3, again finding no detectable companions.

High-Resolution Spectroscopy
Spectroscopic observations enable precise characterization of planet host stars, and can identify false positives caused by spectroscopic binary stars. If sufficiently precise, RV data derived from spectroscopic observations can also further characterize orbits, constrain planet masses, and subsequently lead to the confirmation of the candidates as bonafide planets.

TRES Spectroscopy
We obtained a single reconnaissance spectrum of TOI-4010 on 2021 September 20 from the Tillinghast Reflector Echelle Spectrograph (TRES) through TFOP SG2. TRES has a spectral range from 385 − 910 nm at a resolution of R ∼ 44, 000 over 51 orders (Szentgyorgyi & Furész 2007;Mink 2011) and is located on the FLWO 1.5m Tillinghast Reflector telescope on Mount Hopkins, Arizona. The spectrum suggested that TOI-4010 was a slowly rotating K-dwarf suitable for more precise RV observations.

HARPS-N Spectroscopy
We obtained 112 spectra of TOI-4010 between 2021 July -2023 January with the High Accuracy Radial velocity Planet Searcher for the Northern Hemisphere (HARPS-N; Cosentino et al. 2012) at the 3.6m Telescopio Nazionale Galileo (TNG) at the Roque de los Muchachos Observatory. HARPS-N is a high resolution (R = 115, 000) optical echelle spectrograph with a spectral range of 378−691 nm, whose long-term pressure and temperature stability enable it to reach stability suitable for precise RV work.
The spectra were reduced and had RVs extracted using the HARPS-N Data Reduction Software (DRS) v2.3.5. The RVs were determined by computing the center of a Gaussian fit to the cross correlation function (CCF) of the spectrum with a binary mask (Baranne et al. 1996;Fellgett 1955;Griffin 1967). The weight associated to each line in the mask is proportional to the RV content of the corresponding stellar line (Dumusque 2018;Pepe et al. 2002). The resulting 112 RVs have a mean uncertainty of 5.0 m s −1 and cover a total baseline of 545 days.

STELLAR PARAMETERS
TIC-352682207 is listed in the TICv8.2 with R ⋆ = 0.88 ± 0.06 R ⊙ , M ⋆ = 0.80 ± 0.10 M ⊙ , T eff = 4888 ± 138 K, and log g (cgs) = 4.45±0.09, implying a K dwarf. We further characterized the stellar properties using three different methods. A summary of stellar properties and their sources is given in Table 2.

Spectral Energy Density Analysis
As an independent determination of the basic stellar parameters, we performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia DR3 parallax (Gaia Collaboration et al. 2016) (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the  (2016); Stassun et al. (2017Stassun et al. ( , 2018. We pulled the U BV magnitudes from Mermilliod (2006), the JHK S magnitudes from 2MASS (Skrutskie et al. 2006), the W1-W3 magnitudes from WISE (Wright et al. 2010), and the G, G BP , and G RP magnitudes from Gaia. Together, the available photometry spans the full stellar SED over the wavelength range 0.2-12 µm (see Figure 4). We performed a fit using Kurucz stellar atmosphere models (Kurucz 1993), with the effective temperature (T eff ), surface gravity (log g) and metallicity ([Fe/H]) adopted from the spectroscopically determined values. The remaining free parameter is the extinction A V , which we limited to the maximum line-of-sight value from the Galactic dust maps of Schlegel et al. (1998). The resulting fit ( Figure 4) has a reduced χ 2 of 1.5, with best-fit A V = 0.20 ± 0.05. Integrating the (unreddened) model SED gives the bolometric flux at Earth, F bol = 4.09±0.14×10 −10 erg s −1 cm −2 . Taking the F bol and T eff together with the Gaia parallax, gives the stellar radius, R ⋆ = 0.846 ± 0.023 R ⊙ . In addition, we can estimate the stellar mass from the empirical relations of Torres et al. (2010), giving M ⋆ = 0.89 ± 0.03 M ⊙ , consistent with the (less precise) value obtained from R ⋆ together with the spectroscopic log g, 0.74 ± 0.17 M ⊙ .

Isochrones Fitting
To derive a self-consistent set of physical parameters for the host star, we fit the observed properties of TOI-4010 to MIST evolutionary models (Dotter 2016;Choi et al. 2016) using the stellar model grid package isochrones (Morton 2015). We defined a single star model using parallax from Gaia DR3, observed magnitudes (J, H, K s , G, G BP , G RP , W 1, W 2, and W 3), and T eff and [Fe/H] from the HARPS-N SPC fit. Following the convention of Eastman (2017), we inflated the magnitude uncertainties to 0.02 for the Gaia bands and 0.03 for the WISE bands.
Our fit parameters were stellar mass, age, metallicity, distance, and extinction. We used the default priors from isochrones, except for a broad flat prior for metallicity. 8 We ran the fit using the emcee ensemble sampler Goodman & Weare 2010) with 100 walkers for 40,000 steps, and considered the fit converged based on the chain being more than 50 times longer than the estimated autocorrelation time. We discarded the first 500 steps as burn-in.
Our isochrones fit concluded that TOI-4010 is a metal-rich ([Fe/H] = 0.37 ± 0.07 dex) K-dwarf (T eff = 4960 ± 36 K, log g (cgs) = 4.54 ± 0.02) with mass M ⋆ = 0.88 ± 0.03 M ⊙ , and radius R ⋆ = 0.83 ± 0.01 R ⊙ . All properties are within 1σ of those derived from the HARPS-N SPC and SED analysis. We defer to the isochrones-fitted parameters when requiring stellar parameters in the remainder of this paper.

Stellar Activity
The HARPS-N DRS pipeline computes chromospheric activity indicators from each spectrum, including the Mount-Wilson S-index (S mw ; Noyes et al. 1984;Vaughan et al. 1978), Hα-Index Bonfils et al. 2007), and Na-Index ), which measure emission in the cores of the Ca II H and K, Hα, and Na I D1 and D2 lines, respectively. The DRS also provides indicators that quantify changes in the mean line profile, namely the full full width at half maximum (FWHM) and bisector inverse slope (BIS; Queloz et al. 2001) of the CCF. Prior to further analysis, we removed one negative (and therefore unphysical) measurement from the S mw time series.
We computed generalized lomb-scargle (GLS) periodograms (Zechmeister & Kürster 2009) for all indicators ( Figure 5). There are peaks just above 0.1% false alarm probability (FAP) levels at P ∼ 24 days and P ∼ 33 days in Na-Index and FWHM, respectively. However, we find no significant signals at the periods of the transiting planets. There are also no clear pe-8 Our isochrones priors were as follows. Mass: power law prior with power-law index α = −2.35, bounded between 0.1 and 10 M ⊙ ; Age: flat log prior between 9 and 10.15; Metallicity: uniform prior between -2 and 10 dex; Distance: power law prior with α = 2, bounded between 0 and 3000 pc; Extinction: uniform prior between 0 and 1.
riodicities shared between periodograms, and the time series of the indicators show no long-term trends over the 545-day baseline. We also searched for correlations between indicators and the RVs, recognizing that signals from activity indicators may not necessarily be periodic. In particular, anticorrelation between BIS and RVs can be considered a signature of short-term photospheric variability (Queloz et al. 2001). We find a Pearson rank correlation coefficient of r = −0.16, indicating no significant correlation between BIS and RVs. All other indicators feature weak or non-existent correlations (|r| < 0.3) with the RVs.
Overall, TOI-4010 appears to have low levels of stellar activity, and we conclude that activity is unlikely to influence our measurements of the planet masses.

Stellar Rotation
A rotation signal with a period of P rot ∼ 30 days may be visible in the TESS photometry in Sectors 24 and 25 ( Figure 6). However, the star appears quieter in Sectors 52 and 58, and measuring rotation periods longer than ∼ 13 days using TESS is challenging due to the telescope orbit (e.g. Hedges et al. 2020;Claytor et al. 2022).
Lacking strong periodicities in the activity indicators, we also cannot use the GLS periodograms to find the stellar rotation signal. Instead, we estimate P rot by relating it to log R ′ HK , the logarithmic fraction of the flux in the H and K line cores to photometric contributions from the star. We use the S-index and B − V colour of the star to estimate log R ′ HK using the method outlined by Noyes et al. (1984). Based on the weighted mean of ⟨S mw ⟩ = 0.34 ± 0.06 from the HARPS-N observations and color of B−V = 1.09±0.07 mag from the Fourth US Naval Observatory CCD Astrograph Catalog (UCAC4; Zacharias et al. 2013), we find log R ′ HK = −4.90 ± 0.24, consistent with a low-activity star. This log R ′ HK suggests a stellar rotation period of P rot = 37.7 ± 15.0 days based on the relationship for K stars provided by Suárez Mascareño et al. (2016), which does not correspond to any of the planet orbital periods.

Transit-Only Fit
We first fit our TESS Sector 24 and 25 FFI-extracted light curves, the Sector 52 and 58 SPOC light curves, and ground-based photometry with a three-planet transit model to get precise ephemerides for each planet. We The model was parameterized by the orbital period (P ), transit epoch (T 0 ), planet-to-star radius ratio (R p /R ⋆ ), and impact parameter (b) of each planet, as well as host star mass and radius. We adopted a quadratic limb-darkening law parameterized by q 1 , q 2 from Kipping (2013a), with limb-darkening parameters fit separately for TESS, LCO, and OAA bands. We fit for flux offsets to the FFI and SPOC light curves separately as well as each night of ground-based photometry, and jitter terms added in quadrature to the errors of all observations. We assumed circular orbits, under the expectation that our subsequent fits including HARPS-N RV data would provide significantly better constraints on orbital eccentricity (e) and argument of pericenter (ω) than transit photometry alone.
To evaluate the model's joint posterior, we used PyMC3 ) within the exoplanet package (Foreman-Mackey et al. 2021). We sampled four chains for 4000 tuning steps and 4000 draw steps each, and confirmed that the chains converged according to the Gelman-Rubin convergence statistic for each parameter satisfyingr < 1.01 (Gelman & Rubin 1992).
The transit phase diagrams are plotted in Figure 7. Due to having observations spanning 954 days, the orbital periods and transit epochs are known precisely for all three planets. We improved the uncertainties in P by two orders of magnitude with respect to the Sector 24 and 25 QLP values for all three candidates. The posterior distributions of P and T 0 from this fit are indistinguishable from the results of the final joint model, so we show only the latter results (see §4.3) for simplicity.

RV-Only Fit
We separately fit the HARPS-N RV observations to ensure each planetary signal was present in the data. Visual inspection of the RV data revealed a significant long-term trend, resulting in a peak-to-peak change of ∼100 m s −1 across the observations (Figure 8). While this could be modeled by a simple quadratic trend in time, we do not neglect the possibility that it could be better modeled as a Keplerian due to a long-period giant planet companion. Furthermore, the stellar activity indicators do not show significant signals near the periods expected for the long-period planet ( §3.4), so we do not believe the signal is activity-related. Therefore, we consider five classes of models: • 0p ("no-planet") model: a quadratic trend, • 1p ("one-planet") models: a quadratic trend + a single Keplerian, • 2p ("two-planet") models: a quadratic trend + two Keplerians, • 3p ("three-planet") model: a quadratic trend + three Keplerians, and • 4p ("four-planet") model: four Keplerians.
For one-and two-planet models, we tried all possible combinations of the three transiting planets. For the four-planet models, we set wide uniform priors on the period and mid-transit time of the fourth planet. To test the sensitivity of our fit to choice of prior, we let the period vary between 500 − 600 days, 600 − 700 days, 700 − 800 days, 800 − 900 days, and 500 − 900 days in five separate fits. While the correct orbital period may be longer than 900 days, fits with priors covering longer periods often failed to converge, likely due to the complexity of the system and the lack of sufficient observations to fit an eccentric orbit with a period much longer than the timespan of the data (545 days).
We fixed the periods and transit times of each transiting planet equal to the median of the transit-only fit results, under the expectation that the transit photometry should provide orders of magnitude more precise measurements of these properties than the RV data alone could constrain. This decision significantly sped up the computational time of the model fits. We also relaxed the circular orbit assumption by fitting for the eccentricity and argument of pericenter of each planet's orbit. For eccentricity priors, we used the distribution from Van Eylen et al. (2019) for the three inner planets, and Kipping (2013b) for the long-period giant planet. The argument of pericenter was constrained to the range −π to π, and the sampling was performed in two-dimensional vector space ( √ e sin ω, √ e cos ω) to avoid the sampler seeing a discontinuity at values of π. All models fit for an RV offset and a jitter term which was added in quadrature to the uncertainties. Because the activity indicators for TOI-4010 showed low levels of stellar activity and no strong periodicities, we opted not to fit the RV observations with a more complicated noise model.
We compared the models based on the leave-oneout cross-validation information criterion (LOO; Vehtari et al. 2015) as computed from the exoplanet results using the arviz package for analysis of Bayesian models (Kumar et al. 2019a). Unlike simpler estimates for   for the model, LOO is a fully Bayesian model selection method that uses the entire posterior distribution. A model is favoured if it has a higher log(LOO) score. As shown in Table 3, the three-planet model is a better fit to the data compared to all zero-, one-, and two-planet models, with ∆ log (LOO) = 2.7 ± 2.6 over the next best model (the two-planet model featuring TOI-4010.01 and .02). Combined with the fact that ground-based photometry detected all three transiting signals on-target, high-resolution imaging follow-up did not find evidence for an unresolved eclipsing binary star, and none of the periods of the planets are co-incident with high levels of stellar activity, we consider that the transiting planet candidates from TESS are bonafide planets as confirmed by independent radial velocity measurements. We refer to TOI-4010.01, .02, and .03 as TOI-4010 c, d, and b, respectively, ordered alphabetically by increasing orbital period.
However, the four-planet models are significantly favoured over all other models, as indicated by differences in log(LOO) that are several times larger than the standard errors of the differences. The 500 − 900-day prior gives ∆ log (LOO) = 97.4 ± 10.7 over the threeplanet model. Subtracting the median RV model from the observations also results in the lowest spread in residuals, with root-mean-square (RMS) error of 6.3 m s −1 compared to 14.8 m s −1 for the three-planet case. While the RV semi-amplitude and eccentricity of the fourth planet change slightly depending on the choice of period prior (from K = 49.1 +1.2 −1.2 m s −1 and e = 0.21 +0.03 −0.03 assuming a period between 500 − 600 days, to K = 57.7 +2.5 −2.3 m s −1 and e = 0.29 +0.03 −0.03 assuming a period between 800 − 900 days), all solutions point to a long-period giant planet companion in a moderately eccentric orbit. The fitted properties of the transiting planets also remain consistent across all four-planet fits, indicating that their solutions are robust.
Using the widest prior (500−900 days) to best capture the uncertainties, we find that the fourth planet can be possibly described by P = 760 +91 −90 days, K = 54.5 +3.5 Finally, we inspected the residuals with a GLS periodogram after subtracting the median RV model of each planet from the RVs. The residuals do not show any statistically significant periodicity that could be attributed to additional planets or stellar activity (Figure 9).

Joint transit + RV Model Fit
To obtain the final parameters, we performed a joint analysis of the TESS photometry and HARPS-N RV data using exoplanet. Due to the significantly increased computational expense of the joint fit, we did not include ground-based photometry.
We fit a four-planet model to the RVs based on the results of our RV-only fit, but only a three-planet model to the light curve due to the lack of transit for TOI-4010 e. Our adopted priors are in Table A1, and histograms of the distributions of all fit parameters after sampling four chains with 4000 tuning steps and 8000 draw steps are shown in Figure A1. The final fitted values are reported in Table 4.3, along with additional fundamental transit and physical properties derived from the fits. Because insolation flux and equilibrium temperature require knowledge of the host star's effective temperature, we drew random values from T eff posterior from the isochrones results for each realization of the joint fit. Figure 10 shows the TESS photometry and HARPS-N RV observations folded on the orbits of the three transiting planets, while Figure 11 shows the full RV observations and phase diagram for the nontransiting planet.
In summary, we find the following properties for the four planets of TOI-4010: • TOI-4010 b (TOI-4010.03) is a transiting inner sub-Neptune (P = 1.3 days, R p = 3.02 +0.08 −0.08 R ⊕ , M p = 11.00 +1.29 −1.27 M ⊕ , e = 0.03 +0.03 −0.02 ), Table 3. Results from fitting the HARPS-N RV data alone, testing no-planet (0p), one-planet (1p), two-planet (2p), threeplanet (3p), and four-planet (4p) models. We set the periods of the transiting planets fixed to the transit-only fit results, and used wide uniform priors for the period of the fourth planet in the 4p model. A model is favoured if it has a higher log(LOO) score, demonstrating that the 4p models are the best match to the data. The root-mean-square (RMS) error for each scenario was computed from the difference between the data and the median RV model fit. As discussed in §4.2, we explored orbital periods for TOI-4010 e only up to 900 days, and chose a simple jitter term to capture RV noise rather than more complicated models, such as Gaussian Processes. These choices Figure 10. Plots of the final joint transit + RV model fits, with residuals after subtracting the median models provided in the lower panel of each phase diagram. Black points show observations with offsets subtracted and jitter terms added in quadrature with uncertainties, while coloured circles represent binned data. Solid coloured lines represent median model values, while shaded regions cover 16th to 84th percentiles. The transit (left) and RV (right) phase diagrams for each planet are shown in each row. In each planet's phase diagram, the best-fit models for the other three planets have been subtracted from the data.
were due to the challenges of fitting a four-planet model with incomplete data and indications that TOI-4010 is a low-activity star. We do not expect our treatment of RV noise to affect the fitted RV semi-amplitudes and thus the derived masses of the planets, though it may affect our measurements of the orbital eccentricities (Hara et al. 2019). The RMS error of the RV residuals (6.3 m s −1 ) is also slightly larger than the mean uncertainty of the RVs (5 m s −1 ), indicating that there may be excess scatter due to activity. We leave fitting longer orbital pe-riods and more complicated noise models to when more RV observations become available.

DISCUSSION
The three-transiting-planet system bears a close resemblance to Kepler-18 (Cochran et al. 2011) andK2-19 (Armstrong et al. 2015), which are the only other known systems with at least three planets at periods less than 15 days, two of which are larger than Neptune. All three systems feature two sub-Saturn-sized planets with   Table 4. Results of the four-planet joint fit of transit photometry and RV observations. The central values are the median of the posteriors, with uncertainties computed from the 16th and 84th percentiles. Figure 11. Top: HARPS-N RV time series, which has a long-term trend best explained by an outer giant planet. The four-planet model is over-plotted in red, and the residuals after subtracting the model are in the lower panel. Bottom: RV phase diagram for the long-period giant planet TOI-4010 e with the signals of the other three planets removed, in the same format as in Figure 10. The Keplerian model is plotted in purple. The gap in the phase diagram reflects the fact that the RV observing baseline is shorter than the estimated orbital period of the planet.
a smaller inner planet. However, TOI-4010 is the only of these systems with an additional long-period giant companion, and the only configuration to lack strong resonances, which may be due to dynamical interaction with the moderately eccentric companion. The transiting planets are plotted in mass-radius and period-radius diagrams in Figures 12 and 13 alongside planets listed in the Planetary Systems Composite Data Table 9 on the NASA Exoplanet Archive (2022). TOI-4010 c and d land well within the large spread in typical masses (M p ∼ 5 − 80 M ⊕ ) for sub-Saturn-sized planets (R p = 4 − 8 R ⊕ ), while TOI-4010 b is slightly low-mass for its radius. The orbits of the planets are shown in Figure 14.
5.1. TOI-4010 b: A Sub-Neptune in the Hot Neptune Desert TOI-4010 b is an addition to the sparsely populated region of planets with R p ∼ 2 − 10 R ⊕ and P ≲ 3 days known as the hot Neptune (or sub-Jovian) desert (e.g. Szabó & Kiss 2011;Mazeh et al. 2016). TOI-4010 b is especially unusual because it exists in a multi-planet system. Most planets in the desert lack planet companions (Dong et al. 2018), and TOI-4010 b is so far one of only two P < 2-day hot Neptunes with measured masses and radii that have known planet companions: the recently discovered TOI-969 b (P = 1.8 days, R p = 2.8 R ⊕ , M p = 9.1 M ⊕ ; Lillo-Box et al. 2022) also has a companion, though at a significantly wider separation (a ∼ 2.5 AU, P ∼ 1700 days) than any of the other planets in the TOI-4010 system. As shown in Figure 13, TOI-4010 b's specific location in the desert appears to be in an oasis of hot Neptunes with R p ∼ 3 − 5 R ⊕ , above and below which there are almost no members of the desert with P ≲ 2 days. It is unclear if this is a selection effect, such as a bias in selecting planets for RV follow-up, or an emerging feature of the hot Neptune population.
Both the upper and lower boundaries of the desert have been explained by tidal effects following higheccentricity migration (Matsakos & Königl 2016;Owen & Lai 2018), though this mechanism may not be able to explain planets in or near the desert that are in multiplanet systems such as TOI-4010 b. The lower boundary could instead be sculpted by atmospheric mass-loss due to photoevaporation (Owen & Lai 2018). However, TOI-4010 b still remains in the desert despite being highly irradiated (insolation flux S = 714 +28 −26 S ⊕ ). Furthermore, when compared to theoretical mass-radius curves from Zeng et al. (2018) for planets with R p = 2−4 R ⊕ at 1000 K surface temperature (a proxy for the T eq = 1441 +14 −13 K measured for TOI-4010 b), we find that TOI-4010 b must still be maintaining an atmosphere. TOI-4010 b is most consistent with an Earth-like rocky core with 1 − 2% H 2 envelope by mass, or water-rich core with 0.1 − 0.3% H 2 envelope by mass (included in Figure 12). The combination of hosting an atmosphere while still being highly irradiated may be explained by the fact that TOI-4010 b orbits a metal-rich star, and metalrich atmospheres are less susceptible to photoevaporative mass-loss (Owen & Murray-Clay 2018;Wilson et al. 2022). Figure 12. Mass-radius diagram for planets with density known to better than 50% precision (σρ p /ρp < 0.5), plotted alongside the three transiting, short-period TOI-4010. Planets are coloured by equilibrium temperature, which we re-calculated assuming an albedo of 0 for uniformity. Mass-radius curves from Zeng et al. (2018) assuming 1000 K surface temperature for Earth-like rocky cores with no atmosphere or 0.1%, 1%, and 2% H2 envelopes by mass are overlaid in green, along with water-rich cores with 0.1%, 0.3%, and 1% H2 envelopes by mass (blue) and a pure rock composition (red). Lundkvist et al. (2016) defined an alternative desert for hot super-Earths using insolation flux instead of period, with empirical bounds of R p = 2.2 − 3.8 R ⊕ and S > 650S ⊕ . TOI-4010 b shares this desert with only five other confirmed planets with measured masses and radii: K2-66 b , NGTS-4 b (West et al. 2019), TOI-132 b (Díaz et al. 2020), TOI-849 b (Armstrong et al. 2020), TOI-2196 b (Persson et al. 2022). All of these planets are remarkably different than TOI-4010 b, having two to five times the mass (M p = 21 to 39 M ⊕ ) and up to three times the density (ρ p = 3.1 to 7.6 g cm −3 ) of TOI-4010 b (M p = 11.0M ⊕ , ρ p = 2.2 g cm −3 ). The other members of the desert are also the only planets known in their systems, and their host stars are less metal-rich than TOI-4010 ([Fe/H] = −0.05 to 0.19 dex, compared to [Fe/H] = 0.37 ± 0.07 dex). These findings directly support the negative correlation between planet density and host-star metallicity found for sub-Neptunes by Wilson et al. (2022), further suggesting that high metallicity shaped their histories. TOI-4010 b will be a unique target for future studies on highly irradiated sub-Neptunes.

TOI-4010 c and d: A Pair of Similar Sub-Saturns
TOI-4010 hosts two similarly sized sub-Saturns, one of which (TOI-4010 d) is among the most massive sub-Saturns discovered in a multi-planet system to date (Figure 15). When compared to other systems with multiple short-period sub-Saturns, TOI-4010 stands out further: the combined mass of TOI-4010 c and d is ∼ 58 M ⊕ , greater than the combined masses of the two sub-Saturns orbiting Kepler-18 (∼ 34 M ⊕ ; Cochran et al. 2011) and K2-19 (∼ 43 M ⊕ ;Petigura et al. 2020).
Giant planets are well-established to be more abundant around more metal-rich stars (e.g. Fischer & Valenti 2005;Johnson et al. 2010;Petigura et al. 2018), andPetigura et al. (2017) found a positive correlation between sub-Saturn mass and host star metallicity. These correlations have been interpreted as support for core accretion theory, given that high metallicity disks should have more solid materials available to form more massive cores. TOI-4010 is fittingly one of   Table 4.3.
the most metal-rich sub-Saturn hosts known so far, with [Fe/H] = 0.37 ± 0.07 dex, and only Kepler-27 (Steffen et al. 2012) is a more metal-rich star hosting a multiplanet system with a sub-Saturn. However, as shown in Figure 15, a more updated sub-Saturn population with well-measured densities (precision better than 50%) reveals a weaker correlation between sub-Saturn mass and Figure 15. Mass-metallicity diagram for sub-Saturn-sized planets (Rp = 4 − 8 R⊕) from Figure 12, including TOI-4010 c and d. This population shows a tentative positive correlation between mass and host star metallicity. TOI-4010 d is one of the most massive sub-Saturns discovered in a multi-planet system. TOI-4010 follows the general trend of more massive sub-Saturns orbiting higher-metallicity stars.
host star metallicity compared to Petigura et al. (2017). We find a Spearman's rank correlation coefficient of r s = +0.19, indicating a weak positive correlation. This correlation is slightly stronger among sub-Saturns in multi-planet systems (r s = +0.33), and moderate among sub-Saturns with M p < 60M ⊕ (r s = +0.37). While the relatively large masses of TOI-4010 b and c follow this general trend, additional planets will be needed to better characterize the strength of the correlation.

TOI-4010 e: A Long-Period Giant Planet Companion
Our RV observations reveal an outer companion with a possible orbital period of P = 762 +90 −90 days and minimum mass of M p sin i = 692 +66 −63 M ⊕ (2.2 +0.2 −0.2 M Jup ) consistent with at least a super-Jupiter. At this period, the planet receives an insolation flux of S = 0.15 +0.03 −0.02 S ⊕ , placing it just outside of its host star's habitable zone (Kopparapu et al. 2013).
Based on the fitted next transit time of 3218 +86 −86 TBJD (BJD -2457000) from the RV observations, no transits of TOI-4010 e were expected to occur during any of the TESS sectors over which TOI-4010 was observed. We cannot yet know whether or not TOI-4010 e has a transiting orbital geometry. Without a transit, we can not know the radius, orbital inclination, and true mass of TOI-4010 e. Nevertheless, we will continue to monitor TOI-4010 with RVs to improve orbital phase coverage and better constrain its period and eccentricity.
We explored what constraints Gaia absolute astrometry can provide on the true mass of TOI-4010 e. The Gaia DR3 archive entry for TOI-4010 reports a Renormalized Unit Weight Error (RUWE) value of 0.84. A RUWE > 1.4 (Lindegren et al. 2018(Lindegren et al. , 2021 is typically adopted as a threshold to indicate significant departure from a single-star model and likely evidence of orbital motion in Gaia astrometry. The star thus appears welldescribed as a single star. We then performed numerical simulations following the approach by e.g. Belokurov et al. (2020) to constrain the regime of companions masses and orbital separations that would produce excess astrometric residuals with respect to a single-star model, and therefore a RUWE value larger than the one reported in the Gaia DR3 archive. At the currently assumed separation of TOI-4010 e (a = 1.57 +0.13 −0.12 AU), its true mass would need to be > 6 M Jup in order to record a larger RUWE with 99% probability. Gaia DR3 astrometry thus allows us to rule out quasi-face-on orbital configurations (i ≲ 20 • ) for TOI-4010 e. Significantly improved true mass constraints even in the case of a close-to-coplanar orbit for this massive companion are expected to come from Gaia DR4, slated to be made public in late 2025.

Dynamical Stability
TOI-4010 b and c are within 2% of a 4:1 (third order) mean motion resonance (period ratio 4.0158), but we do not expect there to be significant transit timing variation (TTVs) because no adjacent planet pairs are near strong mean motion resonances. We ran the three transiting, short-period planets through 1000 simulations using ttvfast-python, a Python wrapper for the TTVFast symplectic integrator for TTV computation (Deck et al. 2014). We ran the integrator for 954 days (the total timespan of our TESS photometry) with a time step of 0.067 days (1/20th the period of TOI-4010 b). The median largest TTV amplitude for any planet in the system was less than 2 minutes.
We used the SPOCK stability classifier  to assess the likely long-term stability of the system. SPOCK is an XGBoost-based machine learning classifier trained with 100,000 compact 3-planet systems to quickly predict the probability that a system will remain dynamically stable for 10 9 orbits of the innermost planet, and was directly demonstrated in  to give reliable results for systems of up to five planets. We ran the SPOCK classifier on the posterior from the joint fit under two separate assumptions for the unmeasured inclination of TOI-4010 e: one in which the mutual inclination with the inner system was 0, and another where the mutual inclination was drawn from a uniform distribution between -10 • and 10 • . For each assumption, we computed the stability probability for every draw from the posterior. Both cases showed a high level of stability for the system (with the majority of draws having stability probabilities upwards of 0.84. We find that the least dynamically stable configurations have larger eccentricities for TOI-4010 c and d or a higher mutual inclination between those two planets. We also used full N-body integrations to constrain the feasible parameter space for the orbital inclination of TOI-4010 e. We first randomly selected a draw from the posterior and chose the inclination of TOI-4010 e from a grid spanning a mutual inclination with the inner system between -85 • to 85 • . We initialized each simulation with the drawn parameters, chosen inclination for i e , and M e computed from the drawn value of M e sin i e . We then integrated each set of initial conditions for 1 Myr in rebound (Rein & Liu 2012;Rein & Spiegel 2015) using the MERCURIUS integrator and the additional general relativistic correction implemented in the REBOUNDx gr effect. Almost all of these simulations were dynamically stable for the full integration time, but the integrations with 85 • mutual inclination went dynamically unstable (signified by destruction of the inner system) within a few thousand years. This result suggests that we can use the stability of the inner system to constrain the mutual inclination of TOI-4010 e to be less than 85 • relative to the plane of transiting planets, leaving a very large range of allowable parameter space for the orbital inclination of TOI-4010 e where its presence will not disrupt the inner system.

Follow-Up Feasibility
We computed the Transmission Spectroscopy Metric (TSM; Kempton et al. 2018) for each transiting planet using their derived planet masses. The TSM measures how promising planets are for atmospheric characterization studies. We find TSM = 49, 115, and 50 for TOI-4010 b, c, and d, respectively. When compared to other planets with densities measured to better than 50% precision, TOI-4010 b has a higher TSM than 72% of planets with R p < 4 R ⊕ , and TOI-4010 c and d have higher TSM than 88% and 56% of other sub-Saturns in multi-planet systems, respectively. Importantly, TOI-4010 has multiple promising planets, which can enable atmospheric comparative planetology studies.

SUMMARY
We have presented the confirmation of three transiting exoplanets and the discovery of a fourth planet orbiting TOI-4010, a metal-rich ([Fe/H] = 0.37 ± 0.07 dex) K dwarf observed by TESS in Sectors 24, 25, 52, and 58: • TOI-4010 b (TOI-4010.03; P = 1.3 days, R p = 3.02 +0.08 −0.08 R ⊕ , M p = 11.00 +1.28 −1.27 M ⊕ ) is an addition to the census of planets in the hot Neptune desert, specifically occupying an oasis of hot Neptunes with R p ∼ 3−5 R ⊕ , which may be an emerging feature of the desert. Despite being highly irradiated, TOI-4010 b may still hold on to an atmosphere, likely due to having a metal-rich atmosphere that resists photoevaporative mass-loss. Compared to similarly irradiated planets near its size with both measured mass and radius, it is significantly less massive and dense, and the only one in a multi-planet system.
• TOI-4010 d (TOI-4010.02; P = 14.7 days, R p = 6.18 +0.14 −0.15 R ⊕ , M p = 38.15 +3.27 −3.22 M ⊕ ) is a more massive close-in sub-Saturn, and is the most massive sub-Saturn-sized planet discovered in a system with at least three planets.
• TOI-4010 e (P ∼ 762 +90 −90 days, M p sin i ∼ 692 +66 −63 M ⊕ ) is a moderately eccentric (e ∼ 0.26 +0.04 −0.04 ), long-period giant planet companion detected in our RV observations. These are preliminary physical and orbital properties based on exploring possible orbits with periods up to 900 days due to a lack of sufficient RV baseline. We will continue to obtain observations to better constrain the orbit and determine whether or not the period may be longer. TOI-4010 was followed-up with ground-based photometry, high-resolution imaging, and reconnaissance spectroscopy, and all three transiting planets were confirmed with precise RV observations from HARPS-N. This is only the sixth discovered system with two sub-Saturnsized planets with P < 20 days, and the only known system containing three planets with R p > 3R ⊕ in P < 15day orbits. TOI-4010 c and d are the most massive planets discovered among other rare multi-planet systems with multiple sub-Saturn-sized planets. Sub-Saturns as a population exhibit a large diversity in planet mass at a given size, and two such planets in the same system present an opportunity to study their properties with other variables (e.g. host star age, metallicity) held constant. The TOI-4010 system currently has no known analogues, making it an intriguing target for future work on planet formation and evolution theory.  Table A1. Priors placed on all model parameters for the joint fit of transit photometry and RV observations.