The GAPS programme at TNG XLIX. TOI-5398, the youngest compact multi-planet system composed of an inner sub-Neptune and an outer warm Saturn

Short-period giant planets are frequently found to be solitary compared to other classes of exoplanets. Small inner companions to giant planets with $P \lesssim$ 15 days are known only in five compact systems: WASP-47, Kepler-730, WASP-132, TOI-1130, and TOI-2000. Here, we report the confirmation of TOI-5398, the youngest compact multi-planet system composed of a hot sub-Neptune (TOI-5398 c, $P_{\rm c}$ = 4.77271 days) orbiting interior to a short-period Saturn (TOI-5398 b, $P_{\rm b}$ = 10.590547 days) planet, both transiting around a 650 $\pm$ 150 Myr G-type star. As part of the GAPS Young Object project, we confirmed and characterised this compact system, measuring the radius and mass of both planets, thus constraining their bulk composition. Using multidimensional Gaussian processes, we simultaneously modelled stellar activity and planetary signals from TESS Sector 48 light curve and our HARPS-N radial velocity time series. We have confirmed the planetary nature of both planets, TOI-5398 b and TOI-5398 c, alongside a precise estimation of stellar parameters. Through the use of astrometric, photometric, and spectroscopic observations, our findings indicate that TOI-5398 is a young, active G dwarf star (650 $\pm$ 150 Myr), with a rotational period of $P_{\rm rot}$ = 7.34 days. The transit photometry and radial velocity measurements enabled us to measure both the radius and mass of planets b, $R_b = 10.30\pm0.40 R_{\oplus}$, $M_b = 58.7\pm5.7 M_{\oplus}$, and c, $R_c = 3.52 \pm 0.19 R_{\oplus}$, $M_c = 11.8\pm4.8 M_{\oplus}$. TESS observed TOI-5398 during sector 48 and no further observations are planned in the current Extended Mission, making our ground-based light curves crucial for ephemeris improvement. With a Transmission Spectroscopy Metric value of around 300, TOI-5398 b is the most amenable warm giant (10<$P$<100 days) for JWST atmospheric characterisation.


Introduction
Multi-planet systems give the unique opportunity to investigate comparative planetary science and understand the interactions and processes between their planets (e.g., Dragomir et al. 2019;Kostov et al. 2019;Lacedelli et al. 2021;Trifonov et al. 2023).By studying the relative planet sizes and orbital separations, the obliquity between the planetary orbital planes and the stellar ro-⋆ Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundación Galileo Galilei (FGG) of the Istituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain).tation axis, and other parameters, we can constrain the formation and evolution processes of these systems (see, e.g., Mancini et al. 2022;Grouffal et al. 2022).The precise measurement of the orbital architecture and the bulk compositions of the planets is essential to fully characterise similar systems.In this context, high-precision radial velocities (RVs) are critical to deriving masses and eccentricities, which, combined with the radii from transit, allow measuring precise inner bulk densities and, finally, exploring the differences in planetary structure and evolution.
The architecture of multi-planet systems exhibits a significant amount of diversity, with planets spanning the whole range of masses and grouped in a wide variety of dynamical configura-tions.However, short-period giant planets (P ≲ 10 days) are typically isolated planets compared to other classes of exoplanets (Huang et al. 2016), and usually, their companions, if present, are massive long-period planets with P > 200 days (Knutson et al. 2014;Schlaufman & Winn 2016).Compactness is a rare feature among known multi-planet systems with inner giant planets.Moreover, there is a scarcity of known compact systems composed of a close-orbit outer jovian-size planet and an inner orbit small-size planet (e.g., Hord et al. 2022).Within this family of unique planetary systems, we have WASP-47 (Hellier et al. 2012), Kepler-730 (Zhu et al. 2018;Cañas et al. 2019), TOI-1130 (Huang et al. 2020;Korth et al. 2023), TOI-2000(Sha et al. 2022), and WASP-132 (Hord et al. 2022).
Studying planets younger than 1 Gyr is crucial for understanding the mechanisms at play during the early stages of planetary formation and evolution, such as orbital migration, atmospheric evaporation, planetary impacts, and so on (see, e.g., Terquem & Papaloizou 2007;Bonomo et al. 2019).However, it is challenging to identify and model such systems, because of the magnetic activity of the host star.In fact, both the photometric and spectroscopic time series are affected by the variations in star flux caused by the numerous star spots, faculae, faculae network, and strong magnetic activity on the stellar surface.Despite the stellar activity that hides planetary signals, researchers have discovered planets orbiting members of young stellar clusters (Malavolta et al. 2016;Newton et al. 2019;Mann et al. 2020;Bouma et al. 2020), young stellar associations and moving groups (Benatti et al. 2019;Damasso et al. 2023), and even young field stars (Desidera et al. 2022).There is an increasing number of young exoplanets with well-constrained ages, radii, and masses.Within the Global Architecture of Planetary Systems (GAPS, Covino et al. 2013;Carleo et al. 2020) Young Objects (YO) long-term program at Telescopio Nazionale Galileo (TNG), many systems with young transiting planets, first identified by space missions such as Kepler (Borucki et al. 2010), K2 (Howell et al. 2014), and TESS (Transiting Exoplanet Survey Satellite, Ricker et al. 2015), are being characterised through intense RV monitoring (e.g., Carleo et al. 2021;Nardiello et al. 2022).
Among the systems under intensive scrutiny by GAPS, the moderately young (∼ 650 Myr) solar-analogue star TOI-5398 (BD+37 2118) is of special interest.In this paper, we characterise the compact multi-planet system orbiting this star using a combination of TESS and ground-based photometry, and RVs collected with the High Accuracy Radial velocity Planet Searcher (HARPS-N, Cosentino et al. 2012) spectrograph (Section 2).The TESS pipeline identified two candidate exoplanets: a transiting sub-Neptune (TOI-5398.02,P ∼ 4.77 d) and a giant planet (P ∼ 10.59 d) that has been thoroughly validated by Mantovan et al. (2022) and labelled TOI-5398 b.In Section 3, we report the stellar properties determined using two independent methods.Section 4 reports the procedures used to identify and confirm the two planets in the system by outlining the detailed modelling of photometry and RV.In Section 5, we discuss our results, provide suggestions for follow-up observations, highlight the rarity of our system, and call attention to the exquisite suitability of TOI-5398 b for future atmospheric characterisation.Concluding remarks are in Section 6.

TESS Photometry
The TESS mission observed TOI-5398 (TIC 8260536) at 2 min cadence in Sector 48 from 28 January 2022 to 26 February 2022.In particular, we considered the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016) Simple Aperture Photometry (SAP; Twicken et al. 2010;Morris et al. 2020) flux light curve and corrected it for time-correlated instrumental signatures.We corrected the light curve using Cotrending Basis Vectors (CBVs) extracted through the algorithm by Nardiello et al. (2020) as well as all SAP light curves in the same Camera and CCD as TOI-5398.We did this rather than using the Presearch Data Conditioning Simple Aperture Photometry (PDC-SAP) light curve from the SPOC because, with the latter, the target experienced numerous systematic effects (Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014) ) as a result of over-corrections and/or injection of spurious signals.
In this study, we used light curve points flagged with cadence quality flag values equal to 0 and associated with local background values < 4σ above its mean1 value along the light curve.The only exceptions to this rule are the points before BTJD 2= 2609, which are flagged with a quality flag value equal to 32768 3 and had local background values up to 7σ above the mean.We did this to preserve one of the TOI-5398.02'sfour genuine transits.The corrected light curve shown in Fig. 1 displays a clear modulation that we attribute to the stellar rotation.Therefore, to identify a rotation period of TOI-5398, we computed the Generalised Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) of the aforementioned light curve.The periodogram reveals the most powerful peak at 7.18 ± 0.21 days, see Fig. 2. The link between this peak and the rotational period of the star will be discussed in Sec.3.5.The light curve also exhibits two clear dips at BTJD ∼ 2616 and 2627, which have been shown to be caused by a sub-stellar companion by Mantovan et al. (2022).

ASAS-SN
We downloaded almost three years of archival data from ASAS-SN (Shappee et al. 2014;Kochanek et al. 2017), spanning from October 2020 to May 2023.The ASAS-SN images have a resolution of 8 arcsec/pixel (∼15 ′′ FWHM PSF), and the observations of TOI-5398 were conducted in the Sloan g−band.We checked the light curve to further constrain the stellar rotation period identified from the TESS photometry.We extracted the GLS periodogram, which reveals the most powerful peak with a period of 7.34±0.01days, see Fig. 3.

Photometric Follow-up (TASTE)
A partial transit of TOI-5398 b was observed on January 31st, 2023 by The Asiago Search for Transit timing variations of Exoplanets (TASTE) programme, a long-term campaign to monitor transiting planets (Nascimbeni et al. 2011).To prevent saturation and increase the photometric accuracy, the AFOSC camera, mounted on the Copernico 1.82-m telescope at the Asiago Astrophysical Observatory in northern Italy, was purposely defocused up to 8 ′′ FWHM.A Sloan r ′ filter was used to acquire 2916 frames with a constant exposure time of 5 s.A pre-selected set of suitable reference stars was always imaged in the same field of view in order to apply precise differential photometry.At the end of the series, the sky transparency was highly variable, and thick clouds passed by during the off-transit, leaving us with a subset of 2721 images with a good S/N ratio.We reduced the Asiago data with the STARSKY code (Nascimbeni et al. 2013), a software pipeline specifically designed to perform differential photometry on defocused images.After a standard correction through bias and flat-field frames, the size of the circular apertures and the weights assigned to each reference star were automatically chosen by the code to minimise the photometric scatter of our target.The time stamps were consistently converted to the BJD TDB standard following Eastman et al. (2010), as done for the following photometric data sets as well.

Photometric Follow-up (TFOP)
The TESS pixel scale is ∼ 21 ′′ / pixel and photometric apertures typically extend out to roughly 1 arcminute, generally causing multiple stars to blend in the TESS aperture.To determine the true source of transit signals in the TESS data, improve the transit ephemerides, monitor for transit timing variations, and check the SPOC pipeline transit depth after accounting for the crowding metric, we conducted ground-based lightcurve follow-up observations of the field around TOI-5398 as part of the TESS Followup Observing Program4 Sub Group 1 (TFOP; Collins 2019).We used the TESS Transit Finder, which is a customised version of the Tapir software package (Jensen 2013), to schedule our transit observations.All the image data were calibrated, and all photometric data were extracted using AstroImageJ unless otherwise noted.We used circular photometric apertures centred on TOI-5398, and checked also the flux from the nearest known neighbour in the Gaia DR3 (Gaia Collaboration et al. 2023) and TICv8 catalogues (TIC 8260534), which is ∼ 37 ′′ north-east of TOI-5398.

LCOGT
We observed four and two partial transit windows of TOI-5398 b and TOI-5398.02,respectively, in Pan-STARRS z-short band using the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network nodes at Teide Observatory (TEID) on the island of Tenerife, McDonald Observatory (MCD) near Fort Davis, Texas, United States, and Cerro Tololo Inter-American Observatory in Chile (CTIO).We observed an ingress and egress of TOI-5398 b on UTC 2022 April 21 from TEID and MCD, respectively, an egress on UTC 2023 March 04 from TEID, and an ingress on UTC 2023 March 15 from MCD.We observed an egress of TOI-5398.02on UTC 2022 April 02 from CTIO, and an ingress on UTC 2022 May 10 from MCD.The 1 m telescopes are equipped with 4096 × 4096 Sinistro cameras having an image scale of 0.389 ′′ / pixel, resulting in a 26 ′ × 26 ′ field of view.The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018).We used circular photometric apertures having radii in the range 5.1 ′′ to 7.4 ′′ , which excluded all of the flux from TIC 8260534.

Acton Sky Portal
A full transit window of TOI-5398 b was observed in Sloan i ′ band on UTC 2022 April 21 from the Acton Sky Portal private observatory in Acton, MA, USA.The 0.36 m telescope is equipped with a SBIG Alumna CCD4710 camera having an image scale of 1.0 arcsec/ pixel, resulting in a 17.1 × 17.1 arcmin field of view.We used a circular photometric aperture with radius 9 ′′ , which excluded all of the flux from TIC 8260534.

Whitin
A full transit window of TOI-5398 b was observed in Sloan r ′ band on UTC 2022 April 21 using the Whitin observatory 0.7 m telescope in Wellesley, MA, USA.The 2048 × 2048 FLI ProLine PL23042 detector has an image scale of 0.68 ′′ /pixel, resulting in a 23.2×23.2arcmin field of view.We used a circular photometric aperture with a radius of 8.2 arcsec, which excluded all of the flux from TIC 8260534.

KeplerCam
We obtained photometric observations using the KeplerCam CCD on the 1.2m telescope at the Fred Lawrence Whipple Observatory (FLWO) at Mount Hopkins, Arizona to observe an egress of TOI-5398 b.Observations were taken in the Sloan z ′ band on UT April 21, 2022.KeplerCam is a 4096 × 4096 Fairchild detector which has a field of view of 23 ′ x23 ′ and an image scale of 0.672 ′′ /pixel when binned by 2. We used circular photometric apertures with radius 6.7 ′′ , which excluded all of the flux from TIC 8260534.

MuSCAT2
TOI-5398 was observed by MuSCAT2 on the night of 30 March 2022.MuSCAT2 is a multi-band imager (Narita et al. 2019) mounted on the Telescopio Carlos Sánchez (TCS, 1.52 m) at Teide Observatory, Spain.The instrument is capable of obtaining simultaneous images in Sloan g ′ , r ′ , i ′ , and z s bands with little readout time.Each camera has a field of view of 7.4 ′ × 7.4 ′ and a pixel scale of 0.44 ′′ / pixel.
The observation was made with the telescope defocused to avoid the saturation of the target.On the night of 30 March 2022, the i ′ -band camera presented technical issues and could not be used; the exposure times were set to 10, 5, and 10 seconds in g ′ , r ′ , and z s , respectively.Standard data reduction, aperture photometry, and transit model fit including systematic noise was done using the MuSCAT2 pipeline (Parviainen 2015;Parviainen et al. 2019).

Results
Due to their low S/N and to reduce the computational cost, we decided not to include light curves observed with LCOGT-CTIO and Whitin in our analysis.All other observations resulted in clear transit detections and the light curve data are included in the joint modelling in Sect.4.2 of this work.

HARPS-N Spectroscopic Follow-up
We collected observations of TOI-5398 with HARPS-N at TNG spanning the period between May 2022 and June 2023 and obtained a total of 86 spectra, with exposure times ranging from 900 to 1200 s.These spectra cover the wavelength range 383-693 nm with a resolving power of R ∼ 115 000.We performed the observations within the framework of the GAPS project and incorporated 10 hours coming from a time-sharing agreement (proposal A46TAC_32, PI: G. Mantovan).Additionally, we included eight off-transit spectra from the DDT proposal A46DDT4 (PI: G. Mantovan) in our analysis, which we binned from 600 to ⟨RV err ⟩ (m s −1 ) 3.5 ⟨S/N⟩ 5460Å 64 1200 s of exposure time.Six spectra taken in Spanish time (CAT22A_48, PI: E. Pallé) are also included in a comprehensive analysis of the object.We excluded from our final analysis one of these six spectra as it was collected during the in-transit phase of planet b.We decided to proceed in this way after considering the expected amplitude of the Rossiter-McLaughlin (Ohta et al. 2005) signature, which is further discussed in Sect.5.4.We report the details about the observations and typical S/N in Table 1.
We reduced the data collected using the HARPS-N Data Reduction Software (DRS 3.7.0),and we computed the RV through the Cross-Correlation Function (CCF) method (Pepe et al. 2002 and references therein).With this method, the scientific spectra are cross-correlated with a binary mask describing the typical features of a star with a chosen spectral type.We used a G2 mask for TOI-5398.The resulting CCFs provide us with a representation of the mean line profile of each spectrum.However, the high levels of stellar activity might distort the core of the average line profile, and the stellar rotation broadens the line.Therefore, we needed to proceed with care in selecting the half-window for the evaluation of the CCF and use a width large enough to include the continuum when fitting the CCF profile (see Damasso et al. 2020).We decided to use the G2 mask with a half-window of 40 km s −1 (instead of the default value of 20 km s −1 ) and reprocessed our data using the DRS version implemented through the YABI workflow interface (Hunter et al. 2012) at the Italian Center for Astronomical Archives5 .We obtained RVs with a dispersion of a few tens of m s −1 (∼29 m s −1 ), while their internal errors are approximately a few m s −1 (∼3.5 m s −1 ).
To assess the jitter in the RV series caused by the stellar activity, we additionally extracted a set of activity indices.The value of the CCF bisector span (BIS), as well as the CCF's full width at half-maximum (FWHM) depth and the CCF's equivalent width (W CCF , see Collier Cameron et al. 2019 for further details), are provided by the HARPS-N DRS, while the log R ′ HK index from the Ca II H&K lines was obtained using a method available on YABI (based on the prescriptions of Lovis et al. 2011 and references therein) and using the (B − V) 0 colour index quoted in Sect.3. Finally, we extracted the Hα index using the ACTIN 2 Code6 (Gomes da Silva et al. 2018Silva et al. , 2021)).Figure 4   We also obtained medium precise RVs with the Tull spectrograph at the Harlan J. Smith 2.7m telescope at McDonald Observatory.They can be found in the Appendix.

High angular resolution data
TOI-5398 was also observed on 2 December 2022 with the speckle polarimeter on the 2.5 m telescope at the Caucasian Observatory of Sternberg Astronomical Institute (SAI) of Lomonosov Moscow State University.Speckle polarimeter uses high-speed low-noise CMOS detector Hamamatsu ORCAquest (Strakhov et al. 2023).The atmospheric dispersion compensator was active, which allowed using the I c band.The respective angular resolution is 0.083 ′′ , while the long-exposure atmospheric seeing was 0.6 ′′ .We did not detect any stellar companions brighter than ∆I c = 4.0 and 6.8 mag at ρ = 0.25 ′′ and 1.0 ′′ , respectively, where ρ is the separation between the source and the potential companion (Fig. 5).

Atmospheric parameters and metallicity
Given the relatively young age of TOI-5398, we analysed the HARPS-N combined spectrum following the same methodology as, for example, in Nardiello et al. (2022) and Damasso et al. (2023).Specifically, we adopted an innovative approach to derive the stellar parameters with the equivalent width (EW) method using a combination of iron (Fe) and titanium (Ti) lines.In this way, we avoid issues related to the effect of the stellar activity that can shape the stellar spectrum at young ages (for a detailed explanation we refer the reader to Baratella et al. 2020a,b).Our initial guesses of the atmospheric parameters were estimated with Gaia DR3 (Gaia Collaboration et al. 2016) and 2MASS photometry (Cutri et al. 2003) using the tool colte (Casagrande et al. 2021) and adopting E(B − V) = 0.008 ± 0.016 (Lallement et al. 2014;Capitanio et al. 2017).The photometric estimates vary from 5978 ± 65 K in (G BP − J) to 6039 ± 34 K in (G RP − H).From these initial values and taking the Gaia parallax into account, we derived 4.44 ± 0.09 dex as input surface gravity, while an initial microturbulence ξ = 1.07 ± 0.05 km s −1 was estimated from the relation by Dutra-Ferreira et al. (2016).
The final stellar parameters (T eff , log g, ξ, [Fe/H]) were then derived through the MOOG code (Sneden 1973) and adopting the ATLAS9 grid of model atmospheres with new opacities (Castelli & Kurucz 2003).Our spectroscopic analysis provides as final atmospheric parameters T eff = 6000 ± 75 K, log g = 4.44 ± 0.10 dex and ξ = 1.12 ± 0.12 km s −1 , in excellent agreement with the initial guesses.The derived iron abundance (computed with respect to the solar Fe abundance as in Baratella et al. 2020a) is [Fe/H]= +0.09 ± 0.06, while the titanium abundance is [Ti/H]= +0.08 ± 0.05, where the errors include the scatter due to the EWs measurements and the uncertainties in the stellar parameters.Despite the error bars, there is an indication of a slightly super-solar metallicity.

Projected rotational velocity
With the same code and model atmospheres described in Sect.3.1, we adopted the atmospheric parameters (T eff , log g, ξ, [Fe/H]) previously derived to measure the stellar projected rotational velocity (v sin i ⋆ ).In particular, fixing the macroturbulence velocity to the value of 3.8 km s −1 from the relationship by Doyle et al. (2014), we applied the spectral synthesis method within MOOG for three spectral regions around 5400, 6200, and 6700 Å.Our final value of v sin i ⋆ is 7.5 ± 0.6 km s −1 .We refer the readers to Biazzo et al. (2022) for the procedure based on spectral synthesis and the description of the uncertainty measurement.

Lithium Abundance
Since the lithium line at λ 6707.8Å in the co-added spectrum of TOI-5398 resulted to be blended with the iron line at λ 6707.4Å, we applied the spectral synthesis technique as done in Sect.3.2.Using MOOG and the Castelli & Kurucz (2003) model atmospheres, we fixed the stellar parameters and v sin i ⋆ to the values derived in the previous steps.Then, adopting two different line lists by Carlberg et al. (2012) and Chris Sneden (priv. comm.) and applying the non-LTE calculations of Lind et al. (2009), we obtained a lithium abundance of log n(Li) NLTE = 2.82 ± 0.11, where the error bar considers uncertainties in the line list, in stellar parameters and in the definition of the continuum position around the Li line.With its effective temperature and lithium abundance, the position of our target in the log n(Li)−T eff appears to be in between that of the M35 cluster (∼ 200 Myr) and that of the Hyades cluster (∼ 650 Myr, see, e.g., Sestito & Randich 2005, Cummings et al. 2017).

Chromospheric activity
The lower chromosphere Ca II H&K emission was measured on HARPS-N spectra using YABI (see Sect. 2.5).The average value of the S-index, calibrated in the Mt.Wilson scale (Baliunas et al. 1995), is 0.325, corresponding to log R ′ HK = -4.43 ± 0.02 (arithmetic average and standard deviation).
For the log R ′ HK determination, we adopted (B − V) 0 = 0.58 mag, which was derived from T eff through Casagrande et al. (2006) calibration, as the observed values from Tycho2 and APASS suffer from large uncertainties.Our determination of log R ′ HK implies a stellar age of 370 Myr using the activity-age calibration and an expected rotation period of 4.25 d (corresponding to a gyrochronology age of 260 Myr), through Mamajek & Hillenbrand (2008) relationships.The difference with the observed photometric period (see Sect. 3.5) is not significant considering observational errors and intrinsic scatter of magnetic activity (see e.g.Fig. 7 of Mamajek & Hillenbrand 2008).It is possible that the star is caught by chance at a level of chromospheric activity higher than the average one or that the activity is somewhat enhanced by the presence of planetary companions.
The tidal evolution of the stellar rotation has a timescale much longer than the age of the system, even considering an extremely strong tidal interaction with a stellar modified tidal quality factor Q ′ ⋆ = 10 5 (e.g., Mardling & Lin 2002).Therefore, we do not expect that tides can affect the estimate of the stellar age based on gyrochronology.

Rotation period
The extensive datasets we collected allowed us to obtain several estimates of the rotation period of the star, a key parameter for age determination.
As mentioned in Sect.2.1 and 2.2, from the analysis of TESS and ASAS-SN light curves we obtained photometric periodicities of 7.18±0.21days and 7.34±0.01days, respectively.The larger error of the determination based on TESS data is due to the short baseline of the monitoring (single sector).
We also exploited the HARPS-N time series of RVs and several spectroscopic indicators (Sect.2.5). Figure 4 displays the GLS periodograms, computed for the frequency range 0.0001 -0.5 days −1 , or 2 -10000 days.We removed a linear trend from log R ′ HK and Hα time series before running the GLS.The RVs, BIS, and FWHM periodograms show a clear peak (orange line, normalised GLS power > 0.2) close to the first harmonic of the stellar rotation period detected using the ASAS-SN photometry (7.34 days, dotted green lines).As the period is compatible with the modulation observed in the photometry and consistently recovered in both the radial velocity and the activity indexes, we associated this signal to the stellar activity.This result could be a sign that RV variations are dominated by dark spots rather than the quenching of convective motions in the magnetised regions of the photosphere, as explained in the study by Lanza et al. (2010).Moreover, the periodograms of Hα and W CCF show peaks consistent with the photometric period.On the other hand, the second-highest peak in the RVs periodogram reveals the signal of the validated gas giant planet.We also add that the detailed modelling of the RV time series with Gaussian processing regression (Sect.4.3) yields a value of 7.37±0.03days for the rotation period.
In short, the available photometric and spectroscopic time series consistently indicate a periodicity of 7.34±0.05days, which we interpret as the rotation period of the star.The shorter periodicity seen in some spectroscopic indicators is fully compatible with being the first harmonic of the true periodicity, as observed in several other cases (see, e.g., TOI-1807 Nardiello et al. 2022).Conversely, we think it is unlikely that the true period is the shorter one (3.66 d), with this periodicity arising from the presence of active regions of comparable areas on opposite stellar hemispheres, especially considering the long-time baseline of the ASAS-SN time series.

Kinematics
Kinematic space velocities were derived adopting Gaia kinematic parameters and using the formalism by Johnson & Soderblom (1987).The resulting U, V, and W (Table 2) are at the boundary of the kinematic space of young stars (age younger than about 500 Myr, Montes et al. 2001;Maldonado et al. 2010) in agreement with the other age diagnostics.TOI-5398 is not a member of known young moving groups and a dedicated search of comoving objects within a few degrees in the Gaia DR3 catalogue does not yield convincing candidates.

Age, mass and radius
The indirect methods discussed above point towards an age of a few hundred Myr, similar to the Hyades and Praesepe open clusters, which have also super-solar metallicity (Hyades with +0.15 dex, Cummings et al. 2017, andPraesepe with +0.21 dex, D'Orazi et al. 2020).In the case of a G dwarf star with a few hundred million years, the most robust age indicator is the rotation period.An age of 680 Myr is obtained with the Mamajek & Hillenbrand (2008) calibration for our adopted rotation period.Using different photometric colours and calibrations (e.g., G-K Messina et al. 2022) yields similar results.The Lithium abundance is intermediate between members of Hyades (∼ 650 Myr) and M35 (∼ 200 Myr) of similar colours, pointing to a somewhat younger age.The level of chromospheric activity also suggests a younger age than the gyrochronology, but the discrepancy of both lithium and log R ′ HK with the expectations for an age similar to that of the Hyades cluster is marginal.Finally, kinematic parameters are fully consistent with an age of a few hundred Myr and the lack of comoving objects prevents a more precise age estimate.We further add that isochrone fitting, performed using the param7 tool (da Silva et al. 2006) does not add further relevant information (nominal age 1.6±1.6Gyr).We then adopt a system age of 650±150 Myr from the indirect indicators.
Through the param tool and by imposing the age range allowed by indirect methods to avoid the inclusion of solutions not compatible with the above results (Desidera et al. 2015), we derived the stellar mass and radius.We obtained in this way a stellar mass of 1.146±0.013M ⊙ and a stellar radius R = 1.051±0.013R ⊙ , where the uncertainties are those provided by the param interface and do not include possible systematics of the adopted stellar models.
From the combination of R ⋆ , P rot and v sin i ⋆ , we infer a system orientation fully compatible with edge-on.Indeed, for the nominal parameters, sin i ⋆ is just below unity, and taking error bars into account we estimate i ⋆ ≥ 69 deg.
The stellar parameters outlined in this and the preceding subsections serve as the reference for this study.We present them in Table 2 for reference.Notes. (a) 84th percentile.

Spectral Energy Distribution
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 (see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016); Stassun et al. (2017Stassun et al. ( , 2018)).We pulled the JHK S magnitudes from 2MASS, the W1-W4 magnitudes from WISE (Wright et al. 2010), and the G BP G RP magnitudes from Gaia (Gaia Collaboration et al. 2016), and the FUV and NUV magnitudes from GALEX (Martin et al. 2005).Together, the available photometry spans the full stellar SED over the wavelength range 0.2-22 µm (see Fig. 6).We performed a fit using PHOENIX stellar atmosphere models (Husser et al. 2013), with the free parameters being the effective temperature (T eff ) and metallicity ([Fe/H]), as well as the extinction A V , which we limited to maximum line-of-sight value from the Galactic dust maps of Schlegel et al. (1998).The resulting fit (Fig. 6) has a best-fit A V = 0.05 ± 0.02 mag, T eff = 6025 ± 100 K, [Fe/H] = 0.0 ± 0.5, with a reduced χ 2 of 1.3 (excluding the FUV and W4 fluxes, which suggest modest excesses in both the UV and mid-IR).Integrating the (unreddened) model SED gives the bolometric flux at Earth, F bol = 2.470 ± 0.029 × 10 −9 erg s −1 cm −2 .Taking the F bol and T eff together with the Gaia parallax, gives the stellar radius, R ⋆ = 1.059 ± 0.036 R ⊙ .In addition, we can estimate the stellar mass from the empirical relations of Torres et al. (2010), obtaining M ⋆ = 1.12 ± 0.07 M ⊙ .Finally, we may estimate the stellar rotation period from the above radius together with the spectroscopically measured v sin i ⋆ , giving P rot / sin i ⋆ = 7.1 ± 0.6 d.
These results are fully compatible with those presented in Sects.3.1 and 3.7.

Planets detection and vetting tests
Two candidate exoplanets orbiting TOI-5398 were identified in Sector 48 light curves in both SPOC (Jenkins et al. 2016) and QLP (Huang et al. 2020) pipeline: one giant and one sub-Neptune (P ∼ 10.59 d and P ∼ 4.77 d, respectively).The SPOC performed a transit search with an adaptive, noise-compensating matched filter (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2020)), producing Threshold Crossing Events (TCEs) for which an initial limb-darkened transit model was fitted (Li et al. 2019) and diagnostic tests were conducted to help assess the planetary nature of the signal (Twicken et al. 2018).The QLP performed its transit search with the Box Least Squares Algorithm (Kovács et al. 2002).The two transit signatures passed all TESS data validation diagnostic tests, and the TESS Science Office issued alerts for TOI-5398.01(10.59 days) and TOI-5398.02(4.77 days) on 24 March 2022.The SPOC difference image centroid offsets (Twicken et al. 2018) localised the transit source for TOI 5398.01 within 0.41 ± 2.55 arcsec and for TOI 5398.02within 2.69 ± 2.67 arcsec; all TIC v8 (Stassun et al. 2019) ob-jects other than TOI-5398 were excluded as the source of each transit signature.Mantovan et al. (2022) thoroughly validated the planetary nature of the giant, which they label TOI-5398 b, by ruling out any false positive (FP) scenarios capable of mimicking the observed transit signal.In fact, mainly due to the low spatial resolution of TESS cameras (≈ 21 arcsec/pixel), some objects initially identified as sub-stellar candidates might be FPs.As a result, vetting and validation tests are critical.Therefore, in order to better understand also the nature of the sub-Neptune candidate, TOI-5398.02,and to exclude FP scenarios, we followed the procedure adopted in Mantovan et al. (2022).The latter approach, further described in Appendix A, considers the major concerns reported in Morton et al. (2023) and ensures reliable results when using VESPA (Morton 2012(Morton , 2015)).This procedure allowed us to statistically validate the sub-Neptune exoplanet and label it as TOI-5398 c.

Photometry time-series analysis
To characterise the properties of TOI-5398 b and c, we investigated all the ground-based photometry simultaneously with the two transits of TOI-5398 b and the four transits of TOI-5398 c observed by TESS in a Bayesian framework using PyORBIT8 (Malavolta et al. 2016, Malavolta et al. 2018), a Python package for modelling planetary transits and radial velocities while simultaneously accounting for stellar activity effects.The groundbased observations rule out false positive scenarios caused by BEBs and further confirm the transits of two planetary companions orbiting TOI-5398.
In the present case, we took our TESS corrected light curve (see Sect. 2.1) and carefully considered the influence of stellar contamination from neighbour stars.We verified the stellar dilution by measuring a dilution factor, which defines the total flux from contaminants that fall into the photometric aperture divided by the flux contribution of the target star.Following Sect.2.2.2 of Mantovan et al. (2022), we determined the dilution factor (and its associated error) by calculating the contribution of the flux that falls into the TESS aperture for each star.We determined its value to be 0.00735 ± 0.00005, and we imposed it as a Gaussian prior in the modelling.Then, we selected each space-based transit event from the corrected light curve -and an out-of-transit part as long as the corresponding transit duration (both before the ingress and after the egress) -and created a mask that flags each transit and cuts the corresponding portions of the light curve.The transits of the two planets do not overlap.
We simultaneously modelled each transit (ground-and space-based) using the code BATMAN (Kreidberg 2015), fitting the following parameters: the central time of transit (T 0 ), the planetary-to-star radius ratio (R p /R ⋆ ), the impact parameter, b, the stellar density (ρ ⋆ , in solar units), the quadratic limb-darkening (LD) coefficients u 1 and u 2 adopting the LD parametrization (q 1 and q 2 ) introduced by Kipping (2013), a second-order polynomial trend to take into account the local stellar variability (with c 0 as the intercept, c 1 as the linear coefficient, and c 2 as the quadratic coefficient), and a jitter term to be added in quadrature to the errors of the photometry to account for any effects that were not included in our model (e.g., shortterm stellar activity) or any underestimation of the error bars.We applied an airmass detrending technique to each ground-based light curve, and we estimated u 1 and u 2 using PyLDTk 9 (Husser ) and applying the specific filters used for the observations.We imposed a Gaussian prior on the stellar density, whereas we imposed uniform priors on the period and T 0 .We then imposed a Gaussian prior on the eccentricities following Van Eylen et al. (2019).
We performed a global optimisation of the parameters by executing a differential evolution algorithm (Storn & Price 1997,  PyDE 10 ) and performing a Bayesian analysis of each selected light curve around each transit.The latter was achieved using the affine-invariant ensemble sampler (Goodman & Weare 2010) for Markov Chain Monte Carlo (MCMC) implemented within the package emcee (Foreman-Mackey et al. 2013).We used 4n dim walkers (with n dim being the dimensionality of the model) for 50 000 generations with PyDE and then with 100 000 steps with emcee -where we applied a thinning factor of 200 to reduce the effect of the chain auto-correlation.We discarded the first 25 000 steps (burn-in) after checking the convergence of the chains with the Gelman-Rubin (GR) statistics (Gelman & Rubin 1992, threshold value R = 1.01).Unless specifically mentioned otherwise, the same sampling configuration and process have been used throughout all occurrences of PyDE and emcee.Figures 7, 8, 9, and Table 3 present the results of the modelling.

RV time-series analysis
Using PyORBIT, we investigated the RV time series data in a Bayesian framework.We tried various approaches to model the stellar activity through the use of Gaussian processes (GP, Rasmussen et al. 2006, Haywood et al. 2014).We experimented with a number of data set combinations to constrain the GP hyperparameters and consider various planetary system architectures (one, two, or more planets).Here we outline the three most notable cases.
Due to the high computational cost of the Gaussian Processes, we only modelled the spectroscopic time series in each test.However, we included the inclination from the photometry model to determine the true masses.We modelled the stellar activity in the RV, BIS, and log R ′ HK series simultaneously, through a Gaussian process (GP) regression.We used a quasi-periodic kernel as defined by Grunblatt et al. (2015).As part of this modelling, we set the stellar rotation period P rot (Gaussian prior, as defined in Sec.2.2), the characteristics decay time scale P dec , and the coherence scale w.
In accordance with Eastman et al. (2013), we fitted the periods and semi-amplitudes of the RV signals in the linear space, and we determined the eccentricity e and the argument of periastron ω by fitting √ e cos ω and √ e sin ω.We performed a global optimisation of the parameters by running PyDE and performing a Bayesian analysis of the planetary signals and activity in the RV time series using emcee.The results of this analysis will be discussed after presenting each respective Case and in Table 4  In the present case, we modelled the stellar activity of TOI-5398 using the multidimensional GP framework developed by Rajpaul et al. (2015) and re-implemented in PyORBIT in accordance with the prescriptions in the paper (see also Barragán et al. 2022).Planetary radius (R c ) R ⊕ ...
Again, we relied on the quasi-periodic kernel and its derivatives.We modelled RV, BIS, and log R ′ HK spectroscopic time series.The three-dimensional GP model is the following: where G(t) is the GP and Ġ(t) its time derivative.The constants denoted by subscripts r and c represent free parameters linking the individual time series to G(t) and Ġ(t) (Barragán et al. 2023).
To perform the modelling with PyORBIT, we assigned the same priors described in Case 1 to planets and stellar parameters.Then, we run a global optimisation of the parameters with PyDE and a Bayesian analysis of the planetary signals and activity in the RV time series with emcee.The results of this analysis will be discussed after presenting each respective Case and in Table 4.

Case 3: Extended list of activity indexes and search for a third planet
To better disentangle planetary and stellar signals, we extended the list of activity indexes that PyORBIT simultaneously models with the two keplerian signals (described in Case 2).In particular, we included the chromospheric activity indicator Hα (Gomes da Silva et al. 2011) as well as the two CCF asymmetry diagnostics FWHM and equivalent width W CCF .
On the one hand, the Hα line complements the lower chromosphere indicator log R ′ HK by providing information about the conditions in the upper chromosphere of a star (Gomes da Silva et al. 2011).In fact, Hα and Ca H & K are emitted from different depths -and formed at different temperatures -in the chromosphere (Robertson et al. 2013;Gomes da Silva et al. 2014).Therefore, it is helpful to study the Hα and log R ′ HK indices simultaneously (Gomes da Silva et al. 2011Silva et al. , 2014) ) to learn more about the presence of distinct activity-related features and disentangle their signals from keplerian ones in the RV time series.On the other hand, according to three years of RV monitoring of the Sun (Collier Cameron et al. 2019), the line-shape parameters of the CCF appear to respond to different components of the active regions.Moreover, they help to track global temperature changes in the photosphere (see also Malavolta et al. 2017).
For the reasons listed, we decided to include in the modelling both the chromospheric indicators Hα and log R ′ HK , as well as the three CCF asymmetry diagnostics BIS, FWHM, and W CCF .We followed the multidimensional GP formalism introduced by Rajpaul et al. (2015) to examine the RV and BIS time series, using the first derivative of the GP.Conversely, we did not use the first derivative for the remaining four time series, as suggested also in Barragán et al. (2023) extension of Eq. 1, with the addition of the following supplementary terms: We performed the modelling with PyORBIT in the same way as described in the previous cases, and we show the outcomes in Figs. 10, 11, and Table 5.
We emphasise that the inclusion of the additional activity indicators reduces the uncertainties in the keplerian signals and the RV jitter term with respect to Case 1 and Case 2. Notably, the orbital parameters remain consistent across all cases.We used the Bayesian Information Criterion (BIC, Schwarz 1978) to compare the first two cases.Our analysis revealed a strong preference for Case 2 over Case 1, with a substantial ∆BIC 12 value of 148 (Kass & Raftery 1995).In general, the BIC may not be the optimal estimator for the Bayesian evidence; in this specific case, however, we believe that the extreme difference between the BIC values (∆BIC = 148) largely overcomes any possible bias, thus favouring Case 2 over Case 1.As for Case 3, a direct BIC comparison is not possible due to the different datasets employed in the analysis.Instead, based on logical grounds, we selected Case 3 as our reference model.Each additional activity indicator provides valuable information on specific aspects of stellar activity, as evidenced in previous paragraphs, and their inclusion is justified by the amplitude parameter of the covariance matrix being significantly different from zero for every extra activity indicator, i.e. those activity indicators further constrain the activity model independently from the observed radial velocities.The adopted masses for planet b and c are 58.7 +5.7 −5.6 and 11.8 +4.8 −4.7 M ⊕ , respectively.We include the final parameters of planets b and c in Table 6, while Table 4 lists the differences between this model and the other two.The significantly reduced jitter observed in Case 2 and 3, compared to Case 1, may be attributed to the limited capability of the uni-dimensional GP framework to model the different periodicities present in the spectroscopic time series.In particular, the RV and chromospheric indexes exhibit different periodicities (see also Sec. 3.5 and Fig. 4) when the RV variations are dominated by dark spots.In contrast, Case 2 and 3 use the formalism outlined in Rajpaul et al. (2015), which effectively handles variations in periodicity.
In addition to the modelling mentioned above, we looked for the presence of a third planet in our RV dataset, by applying wide uniform priors on its period P and RV semi-amplitude K.While we did not impose the orbit of the third planet to be circular, we applied a Gaussian prior on the eccentricity following Van Eylen et al. (2019).Initially, the global optimization algorithm and the first ten thousand steps of the MCMC exploration suggested a significant detection of a third planet.However, the chains of the planet's orbital period diverged quite soon, preventing us from claiming a third planet detection.We underline that the solutions for planets b and c in the system show little variation, which further strengthens the validity of their detections.We computed the observed (O) -calculated (C) diagrams for both planets removing the linear ephemeris (in Table 6) to each transit times.See the O-C diagrams in Fig. 12 and 13 for planets b and c, respectively.The possible TTV amplitude (A TTV ), computed as the semi-amplitude of the O-C, is of 2.9 +1.2 −1.0 minutes for planet b and 4 +7 −2 minutes for planet c.The associated error is derived from the subtraction of A TTV from the high-density interval at 95% of a Monte-Carlo sampling of 10 000 repetitions.
Although our data do not cover a sufficient portion of the super-period to offer direct evidence, they do suggest a possible TTV due to the gravitational interaction of planets b and c (see Fig. 12).The sparse sampling of the TTV signals did not allow us to run a dynamical fit, so we decided to run a forward dynamical model with TRADES11 (Borsato et al. 2014(Borsato et al. , 2019b(Borsato et al. , 2021)), similar  to what has been done in Tuson et al. (2023).We took the masses and the orbital parameters from Tables 5 and 6 and we integrated the orbits for about 500 days and computed the transit times of each planet.We computed the O-C diagrams (see Fig. 14) and we found that the simulated A TTV is of the order of ∼ 40 seconds and of ∼ 2 minutes for planet b and c, respectively.We used the orbital parameters from Table 6 as the starting point of the dynamical simulation, i.e., we assumed that the value we obtained from our global fit spanning about 500 days represents a specific configuration in time, which may not be true.The amplitude of the simulated TTVs may also be dependent on   the eccentricity of the two planets.Nevertheless, the outcome of the simulation is compatible (at 2σ for b and at 1σ for c) with the measured TTVs.The scope of this analysis is to show that TTVs can indeed be present in this system, while a full dynamical analysis is outside the scope of this paper and it is left to future works.5. T eq (K) 947±28 1242±37 Notes. (a) 84th percentile. (b) 16th percentile.

Peculiar architecture
TOI-5398 is a compact multi-planet system composed of a warm giant (TOI-5398 b, P orb ∼ 10.59 d) and a hot sub-Neptune planet (TOI-5398 c, P orb ∼ 4.77 d) orbiting a moderately young solaranalogue star.The peculiarity of this system comes from its com-pactness and planetary architecture, which is uncommon among known multi-planet systems with short-period giant planets.In fact, the sub-Neptune is the closest planet to the host star.There are currently very few notable examples of compact systems consisting of an inner orbit small-size planet and an outer shortperiod giant companion (e.g., Hord et al. 2022).The most famous multi-planet system with similar characteristics is WASP-47 (Hellier et al. 2012;Becker et al. 2015;Nascimbeni et al. 2023).Within this family of peculiar planetary systems, we also have Kepler-730 (Zhu et al. 2018;Cañas et al. 2019), TOI-1130 (Huang et al. 2020;Korth et al. 2023), TOI-2000(Sha et al. 2022), and WASP-132 (Hord et al. 2022).We list them in Table 7, together with a list of available datasets and a few notable parameters that we want to highlight.Among these compact systems, TOI-5398 stands out, along with WASP-47 and TOI-2000, due to its precise transit photometry and radial velocity measurements.These observations enable the precise measurement of planetary bulk densities and make it an extremely appealing target for continued monitoring with follow-up observations and surveys, including PLATO and Ariel, along with telescopes like ELTs.
Moreover, TOI-5398 is the youngest compact system with a gas giant ever confirmed with a quite good age estimation.

Variety among compact multi-planet systems
In Fig. 15, we show the architecture of compact multi-planet systems with small-size planets orbiting interior to short-period giant planets (P ≲ 10 d).For comparison, we show the next three systems with giants whose orbital periods are P < 25 d (Butler et al. 1997;Bourrier et al. 2018;Weiss et al. 2013;Nesvorný et al. 2013).The semicircular dots represent the host stars, colour-coded by their age, while the size represents their radii.Planets, on the other hand, are colour-coded by their equilibrium temperature T eq , and their size reflects their planetary masses.The inner planet of WASP-132 and the Kepler-730 planets do not yet have mass measurements, so we extracted them following Wolfgang et al. (2016) or we showed upper limits (Hord et al. 2022).
TOI-5398 hosts the youngest gas giant planet with P < 25 d and M p > 1/2 M Saturn that is known to have an inner companion.The gas giant TOI-5398 b has a radius similar to Jupiter and a mass close to 2/3 of Saturn, which is the smallest mass among giants in compact systems.As a result, its bulk density is around half that of Saturn and roughly equal to that of TOI-1130 c (0.38 g cm −3 , Huang et al. 2020).These properties, as well as its orbital period, make TOI-5398 b quite similar to the hot-Saturn planet TOI-2000 c, while they are in contrast with the properties of the larger WASP-47 b, Kepler-730 b, and TOI-1130 c, which all have radii ∼ 1.1 R J and masses ∼ 1 M J (apart from Kepler-730 b that does not have a mass measurement yet).However, the larger radius and smaller mass of TOI-5398 b compared to those of TOI-2000 c -which results in having a relatively low density compared to giant planets of similar mass (see Fig. 16 and cf., for example, Yee et al. 2022) -make the giant planet under study more like a puffy Saturn (Naponiello et al. 2022).Moreover, we calculated the equilibrium temperature T eq (cf.Eq. 4 from Cowan & Agol 2011) of TOI-5398 b, assuming zero albedo and full day-night heat redistribution following Fig. 15.Architecture of compact multi-planet systems hosting smallsize planets orbiting inner to short-period gas giants (P ≲ 10 d).Each row represents one planetary system (y-axis) and the planetary orbital periods (x-axis).The sizes of the dots correspond to the planet masses, and the colours of the points to the equilibrium temperatures (see colorbar to the right).From top to bottom, the systems are sorted in ascending order of the period of the giant.Shaded dots represent the next three systems with gas giants on P < 25 d orbits.Dots with a vertical line represent planets whose mass is multiplied by sin(i).The semi-circular dots filled with a star shape stand for the host stars, colour-coded by their age (see colorbar to the top), while the size represents their radii.
where a is the orbital semi-major axis given in the same units as R ⋆ .We obtained a value of 947±28 K, which indicates that TOI-5398 b is unlikely to be affected by the hot-Jupiter anomalous radius inflation mechanism (Thorngren et al. 2016).By contrast, the T eq of the hot-Saturn TOI-2000 c is slightly above 1000 K (Sha et al. 2022); therefore, considering the T eq = 1000 K threshold and its orbital period exceeding 10 days, we describe TOI-5398 b as a warm-Saturn planet.
The sub-Neptune TOI-5398 c is one of the very few inner companions to short-period gas giants with precise values of both mass and radius.In fact, only four small planets with R p < 4 R ⊕ in compact systems (P ≲ 15 d) have mass measurements: WASP-47 d, WASP-47 e (Vanderburg et al. 2017), TOI-1130 c (Korth et al. 2023), and TOI-2000 b.It is worth noting that they have quite different bulk densities (see Fig. 17), ranging from being composed of rocky cores to having masses and radii similar to those of Neptune and Uranus.The latter composition is true for TOI-5398 c, which shares a bulk density similar to that of the inner companion planet WASP-47 d.As a result, we provide additional evidence that inner companions to transiting giant planets tend to have the same density diversity as other small planets (Sha et al. 2022).

Ephemeris improvements
An important step of our analysis is the derivation of new and updated mean ephemeris for TOI-5398 b and c.Our best-fit relation for the warm Saturn and the sub-Neptune are:   T 0,c = 2459628.6178± 0.0009 BJD TDB + N × (4.77271 ± 0.00016), where the variable N is an integer number commonly referred to as the "epoch" and arbitrarily set to zero at our reference transit time T ref .
We emphasise that if we propagate the new ephemeris at 2030-01-01 (see Fig. 18, and 19), the level of uncertainty is significantly reduced to ∼ 5 minutes compared to the previous ∼ 197 minutes for TOI-5398 b when only TESS photometry was available.This means that when also the ground-based photometry is taken into account, the error bar for TOI-5398 b is 98% smaller compared to using TESS data alone.For TOI-5398 c, the error bar is 60% smaller than when we rely solely on TESS data.Accurately identifying the transit windows is crucial for upcoming space-based observations, given the significant investment in observing time and the time-critical nature of such observations.It is crucial to note that no further observations of TOI-5398 are planned in the current TESS Extended Mission12 .2019): solid brown indicates an Earth-like rocky core (32.5% Fe and 67.5% MgSiO3), beige a 100% water world at 1000K, and dotted beige a 1% hydrogen envelope and 99% Earth-like rocky core at the same temperature.Grey dashed curves represent densities ρ = 0.5, 1, and 5 g/cm 3 , respectively.

Planetary system formation and evolution
The obliquity between the planetary orbital plane and the stellar rotation axis is a key diagnostic for the mechanisms of formation and orbital migration of exoplanets (e.g., Naoz et al. 2011).This can be detected with in-transit RVs through the Rossiter-McLaughlin effect (RM, Ohta et al. 2005;Rossiter 1924;McLaughlin 1924;Queloz et al. 2000).Short-period giant planets are thought to form in situ close to the final orbit, or in the outer regions and migrate inward (Dawson & Johnson 2018).Different mechanisms, such as dynamical interactions (high-eccentricity migration) through planet-planet scattering (Marzari et al. 2006) or the Kozai mechanism (Wu & Murray 2003), and disc-planet interaction (Lin et al. 1996), can shrink their orbits.These mechanisms are expected to imprint different signatures in the planets' obliquity.Scattering encounters should randomise the alignments of the orbital planes, while migration through disc-planet interactions should keep the planetary orbits roughly co-planar throughout the entire process.
TOI-5398's uncommon architecture and moderately young age make it particularly promising for measuring the obliquity between the orbital plane of the giant and the spin axis of the star.First, we can access the original configuration when observing systems young enough to have avoided tidal alterations of the obliquity.Then, unlike ordinary short-period giants, we can rule out the high-eccentricity migration scenario (Mustill et al. 2015) for compact systems such as TOI-5398 and test the other formation models through detailed atmospheric characterisation.
Following Eq. 40 from Winn (2010), we determined the expected amplitude of the radial velocity variation produced by the RM effect when planet b transits (58 m s −1 ) or when planet c transits (7.0 m s −1 ).Given the activity level of our target (typical RV dispersion: 27 m s −1 ) and considering a planetary transit time scale, that is, when the activity level is significantly lower than its typical amplitude, we predict a perfectly suitable detection of the RM effect caused by TOI-5398 b and a possible detection for TOI-5398 c.

Circularisation time scale and α parameter
We calculated the circularisation time scale, denoted as τ circ , for TOI-5398 b, using Eq.6 from Matsumura et al. (2008).Assuming a circular orbit and a modified tidal quality factor Q of 10 5 (Ogilvie 2014), the calculated value for τ circ is 2.57 ± 0.88 Gyr, which significantly exceeds the age of the system.This suggests that the near-circular orbit of TOI-5398 b may well be primordial, indicating favourable conditions for preserving its close planetary companion.
Following the methodology described in Bonomo et al. (2017), we computed the parameter α, which represents the ratio of the planet's semi-major axis to its Roche limit.Bonomo et al. (2017) conclude that planets with α > 5 and circular orbits are unlikely to undergo high eccentricity migration.With an α value of 5.6, TOI-5398 b falls in the middle α range (4.55 -7.44) for short-period giants with close companions.This finding favours the notion that short-period giants with close companions should not be distinguished between hot (P < 10 days) and warm (P > 10 days) planets, that is, they belong to the same population of exoplanets, distinct in turn from the typical giants that experience the high-eccentricity migration scenario.

Atmospheric mass-loss
We investigated how the planetary masses, radii, and atmospheric mass-loss rates change with time due to photoevaporation and internal heating.
We evaluated the mass-loss rate of the two planets' atmospheres using the hydro-based approximation developed by Kubyshkina et al. (2018a,b), coupled with the planetary coreenvelope model by Lopez & Fortney (2014) and the MESA Stellar Tracks (MIST; Choi et al. 2016).For the stellar X-ray emission at different ages, we adopted the analytic description by Penz et al. (2008), with the current X-ray luminosity, L x = 10 29 erg s −1 in the 5-100 Å band.This value was derived from the rotation-activity relationships by Pizzolato et al. (2003), and closely resembles the median X-ray luminosity of Hyades stars.The stellar EUV luminosity (100-920 Å) was computed at any given time using the scaling law by Sanz-Forcada et al. (2022).The subsequent paragraphs present the outcomes of our forwardin-time simulations, which include the evolution of the XUV irradiation and the planetary structure in response to stellar behaviour.More details on our modelling of atmospheric evaporation are provided in Maggio et al. (2022) and Damasso et al. (2023).
We performed several simulations assuming different possible values for the planetary core radius at the current age.The test cases were selected by comparing them with the grids of planetary internal structures by Fortney et al. (2007), assuming cores composed of 50% rocks and 50% ices.For planet b, we explored core radii ranging from 2.4 to 7 R ⊕ , corresponding to core masses of 10 to 25 M ⊕ .As for planet c, due to its smaller size, we limited our analysis to core radii between 1 and 2.5 R ⊕ , with core masses of ∼ 10M ⊕ .In Table 8, we show the results of the simulations.
In our reference model for planet b, the core has a radius of R core = 5R ⊕ and a core mass of M core ∼ 29M ⊕ , resulting in an atmospheric mass fraction f env of ∼ 50%.The current photoevaporation rate is ∼ 5.3 × 10 11 g s −1 , and the planet will maintain a large envelope mass fraction throughout its main-sequence lifetime, with f env ∼ 48%M p at time t ∼ 5 Gyr.Its radius will only be reduced by ∼ 10 %.In the range explored, these results depend little on the assumed characteristics of the core.
Conversely, the evolution of the inner planet is very different due to the smaller distance from the host star, higher equilibrium temperature, and higher high-energy irradiation.Our reference model has a core radius R core = 1.8R ⊕ and a core mass M core ∼ 10M ⊕ , resulting in an atmospheric mass fraction f env of Notes. (a) The evaporation time scale is equal to the time required for losing half of the atmospheric mass.
∼ 3.5%M p .The current photo-evaporation rate is ∼ 3.5 × 10 12 g s −1 , and the planet will lose its entire envelope in ≲ 200 Myr from now.The planetary size will decrease to match that of the core.However, a larger core radius implies a smaller atmospheric mass fraction and shorter evaporation time scales.For example, for a core radius R core = 1R ⊕ , the planet would keep a residual atmospheric envelope even at t = 5 Gyr.

TOI-5398's global formation history
The different masses of the two planets could be either primordial or, as suggested by their different evaporation rates discussed in Sect.5.4.2, the result of distinct photoevaporation histories.However, these two scenarios have different implications regarding formation regions and bulk compositions.Our preliminary exploration of the formation tracks of the two planets using the methodology from Johansen et al. ( 2019) highlights the following possibilities.
For the mass of TOI-5398 c to be primordial or close to the original one, the planet should not have captured significant quantities of disk gas after reaching its pebble isolation mass.This condition is satisfied if the planet had started its formation beyond about ten au and comparatively late (∼ 1 Myr) in the life of its native disk.However, in this scenario, TOI-5398 b should still be able to capture its present gaseous envelope.For this to occur, planet b should have started its formation at an earlier time (∼ 0.1 Myr) than planet c: as a result, planet b would already be close to its current orbit while planet c is forming and migrating, and the two planets would have to cross orbits to reach their current architecture.Such an encounter would likely result in a planet-planet scattering event: this scenario would result in higher eccentricities and inclinations than the currently observed ones.Even in the case that the dynamical excitation created by the planet-planet scattering event is removed by the interactions of the two planets with the disk gas, the density of planet c appears too low for a realistic mixture of rock and ice resulting from its growth track (as a comparison, the density of the ice-rich dwarf planet Pluto is 2 g cm −3 ).
The other possible scenario results from the two planets having started forming early in the lifetime of their native disk at a few au from their host star.In this case, the simulated growth tracks favour planet c to have possessed an extended primordial atmosphere of comparable mass to that of planet b.In this scenario, the present-day mass disparity between the two plan-ets would result from their different photoevaporation histories, with planet c experiencing a much higher mass loss than its outer counterpart.Their gas accretion phases would have occurred close to their final orbits, in the innermost and hottest regions of the native disk, which suggests that their atmospheres could exhibit stellar composition.

Planetary bulk composition prediction
Accurate data on mass, eccentricity, and radius allow for measuring precise inner bulk densities and exploring the differences in planetary structure and evolution, from inflated Hot Jupiters (HJ) to "over-dense" Warm Jupiters (WJs, Fortney et al. 2021).While the prediction of hotter interiors and larger radii for HJs (Guillot et al. 1996), compared to Jupiter, has been proven true, the mechanism(s) behind some HJs having anomalously large radii remain a challenge to explain (Thorngren & Fortney 2018).For non-inflated giant planets (F ⋆ < 2 × 10 8 erg s −1 cm −2 ), Thorngren et al. (2016) found a relation between planet mass and bulk metallicity, which confirms a key prediction of the core-accretion planet formation model (Mordasini et al. 2014).A recent study (Müller & Helled 2023b) presents the current knowledge of mass-metallicity trends for warm giant exoplanets (Teske et al. 2019;Müller et al. 2020;Müller & Helled 2023a), and raises some doubts about its extent and existence.Müller & Helled (2023b) link this ambiguity to theoretical uncertainties on the assumed models and the need for accurate stellar age and atmospheric measurements.
Understanding this relationship and open questions regarding giants require characterising planets and host stars, focusing on the metal enrichment of planetary atmospheres (Miller & Fortney 2011).It is particularly true for warm giants -scarce among confirmed planets13 -unaffected by the radius inflation mechanism, as we can reasonably constrain their bulk metal enrichment and interpret atmospheric features more safely (Thorngren et al. 2016).
TOI-5398 b is a perfect case study if we consider both its low insolation F ⋆ and its good age estimation.Therefore, we estimated its planetary bulk heavy-element mass fraction (or bulk metallicity), using evolution models14 from Müller & Helled (2021).Figure 20 shows the radius evolution for vari-ous bulk metallicities (coloured lines) for three different atmospheric heavy-element mass fractions.The bulk heavy-element mass fraction of TOI-5398 b, depending on the adopted atmospheric metallicity, varies between 20 and 30 per cent of its total mass.Our result follows the mass-metallicity trend from Thorngren et al. ( 2016), but the heavy-element mass appears to be slightly lower than expected.Therefore, we may infer that the trend from Müller & Helled (2023a), which predicts a lower heavy-element mass for a given planetary mass compared to Thorngren et al. (2016), might offer a more plausible explanation of our result.

Future atmospheric characterisation
To break the degeneracy in determining the planetary bulk composition, it is crucial to perform atmospheric measurements and to get information on metal enrichment.Müller & Helled (2023b) show that atmospheric measurements by JWST and Ariel can significantly reduce this degeneracy and that this is particularly promising for warm giant planets (Müller & Helled 2023a).The precise characterisation of the TOI-5398 b atmosphere will be crucial to validate or disprove formation and evolution theories.
TOI-5398 is a fascinating compact system, as planet b has the highest Transmission Spectroscopy Metric (TSM, Kempton et al. 2018) value (∼ 300) among warm giant planets (10 < P < 100, M p > 0.1 M J ) currently known, making it ideal for atmospheric characterisation by JWST.Indeed, the TSM parameter quantifies the expected signal-to-noise in transmission spectroscopy for a given planet, and according to Kempton et al. (2018), a giant planet's atmosphere is considered amenable to JWST observations when its TSM value is greater than 90.In Fig. 21, we include all confirmed planets, with colour-coding only for those with a TSM > 90, and M p > 0.1 M J .In this plot, we show the planetary masses versus radii, where planets are colourcoded by their orbital period.We colour-code planets with periods longer than 10 days in orange to highlight the warm-giant planets, defined as giant planets with periods exceeding 10 days (e.g., Yee et al. 2021;Gan et al. 2023 and references therein).Instead, the size of the dots tracks the TSM.For comparison, WASP-47 b has only a modest TSM value of ∼ 47 (Bryant & Bayliss 2022).

Characterisation with JWST/NIRSpec
TOI-5398 b is an ideal candidate for precise transmission spectroscopy and atmospheric characterisation, with a focus on the carbon-to-oxygen (C/O) ratio.This ratio is essential for understanding planetary formation mechanisms and grasping planetary atmosphere composition, shedding light on volatile content and atmospheric chemistry.The availability of carbon and oxygen determines different chemical reactions as well as the stability of molecules of planetary atmospheres.Understanding the C/O ratio helps in predicting the composition and behaviour of atmospheric constituents like carbon dioxide (CO2), carbon monoxide (CO), methane (CH4), and water (H2O) (Keyte et al. 2023).
To test the feasibility of atmospheric characterisation using JWST, we investigated three different atmospheric scenarios for TOI-5398 b.We assumed equilibrium chemistry as a function of temperature and pressure using FastChem (Stock et al. 2018) and three different C/O ratios: 0.5, 1.0 and 1.5.We used FastChem within TauREx3 (Al-Refaie et al. 2021) using the taurex-fastchem15 plugin.TauREx is a retrieval code that uses a Bayesian approach to infer atmospheric properties from observed data, utilising a forward model to generate synthetic spectra by solving the radiative transfer equation throughout the atmosphere.We used all the possible gases contributions within FastChem and the active absorption contribution given by K, Na, HCN, H2CO, CH4, CO, CO2, C2H2, C2H4, H2O, NH3, SiO, TiO, VO and SO2 opacities.
After generating the transmission spectra using Tau-REX+FastChem, we simulated a JWST observation using Pandexo (Batalha et al. 2017), a software tool specifically developed for the JWST mission.The software allows users to model and simulate various atmospheric scenarios, incorporating factors such as atmospheric composition, temperature profiles, and molecular opacities.We simulated a NIRSpec observation in bots mode, using the s1600a1 aperture with g235h disperser, sub2048 subarray, nrsrapid read mode and f170lp filter.We simulated one single transit and an observation 1.75 T 14 long to ensure a robust baseline coverage.We fixed this instrumental configuration for all three scenarios.In Fig 22, we show the resulting spectra for the different C/O ratios and their best-fit models.
We performed three atmospheric retrievals on the NIR-Spec/JWST simulations using a Nested Sampling algorithm with the nestle16 library with 1000 live points.We used the cuda transmission model with the taurex-cuda17 TauREx plugin.We fitted three parameters: the radius of the planet R p , the equilibrium temperature of the atmosphere T eq , and the C/O ratios.
Using NIRspec with the g235h disperser wavelength range (1.66 µm -3.07 µm), we can assess the C/O ratio under the three assumptions (see Tab 9).In particular, when assuming C/O ratios of 0.5 and 1.0, we can retrieve the correct value within 1σ error bar, while under the C/O = 1.5 assumption, we retrieved a value of 1.87 ± 0.15, within 2.5σ error bar.
The three atmospheres can be explained with three distinct sets of parameters (Fig 23).The results of atmospheric retrievals confirm and quantify the feasibility of atmospheric characterisation using NIRSpec@JWST.Furthermore, they demonstrate that TOI-5398 b is an excellent candidate for comprehensive atmospheric analysis, to measure the C/O ratio and, therefore, to constrain planet formation theories for this system.

Conclusions
In this study, we presented the discovery of the youngest transiting planetary system containing a sub-Neptune planet orbiting interior to a Saturn-mass planet with P < 15 days.Using HARPS-N radial velocity measurements of the host star TOI-5398 and multidimensional Gaussian processes, we modelled the stellar activity and confirmed the planetary nature of both candidates identified in the TESS light curve measuring their masses.Furthermore, our methodology allowed us to accurately determine the stellar parameters.With a Transmission Spectroscopy Metric value of around 300, the warm Saturn TOI-5398 b is the most suitable warm giant planet for atmospheric characterisation using JWST.Such investigations are crucial to validate or disprove existing formation and evolution theories.By measuring atmospheric chemistry, we can gain information on metal enrichment and effectively break the degeneracy in determining the planetary bulk composition.
The presence of two planets in this relatively young system offers the opportunity to examine distinct evolutionary paths over the initial hundreds of millions of years following system formation, under the influence of the same host star.In this study, we provided a characterisation of the system, with a special focus on the future evolution of the planetary atmospheres.Future works will focus on investigating the past evolution of the system.We explored the evolution of the atmospheres of both planets, considering the decay of stellar activity and XUV irradiation with time.We estimated that given reasonable assumptions regarding core radius and mass, planet b probably retains a substantial atmosphere, with a mass fraction ∼ 0.5M p , making it amenable to be probed with transmission spectroscopy.Conversely, planet c is expected to possess a tiny atmospheric envelope.At the current age, the mass loss rate of planet c exceeds that of planet b by a factor of 7, implying that planet c will completely lose its residual atmosphere within a few hundred million years, while planet b will retain a thick atmosphere even at the solar age.
Notably, TESS observed TOI-5398 during sector 48 and no further observations are planned in the current Extended Mission.Consequently, our ground-based light curves play a pivotal role in refining the ephemeris of both planets.The improved ephemeris we calculated are vital for future follow-up observations and surveys, including those conducted by CHEOPS, JWST, and upcoming missions such as PLATO and Ariel, along with telescopes like ELTs.

Fig. 1 .
Fig. 1.SAP-corrected light curve of TOI-5398, observed in Sector 48.The local background value is colour-coded.Yellow points -with local background values > 4σ above its mean value along the light curvehave been excluded from all analyses.

Fig. 2 .Fig. 3 .
Fig. 2. GLS periodogram of the TESS SAP photometry corrected for time-correlated instrumental signatures.The vertical lines indicate the period at 7.18 days and its first harmonic.
displays the spectroscopic time series.

Fig. 4 .
Fig. 4. Left: HARPS-N spectroscopic time series used in this work.The time series for RV, BIS, log R ′ HK , Hα, FWHM, and W CCF are shown in the panels in order from top to bottom.Right: GLS periodogram of the RVs and the spectroscopic activity indicators under analysis.The primary peak of each periodogram is indicated by a vertical orange line.The dotted green lines represent the stellar rotation period described in Sect.2.2.The signals along the red dotted vertical lines correspond to the transit-like signals with periods 4.77 and 10.59 d.

Fig. 5 .
Fig. 5. High angular resolution speckle imaging of TOI-5398 in I c filter using SAI-2.5 m telescope.

Fig. 6 .
Fig. 6.Spectral energy distribution of TOI-5398.Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband.Blue symbols are the model fluxes from the best-fit PHOENIX atmosphere model (black).The Gaia spectrum is overlaid as a grey swathe and shown in closer detail in the inset plot.

Fig. 7 .
Fig. 7. Photometric modelling of TOI-5398 b planetary signal.In the top panel, we display the TESS phase-folded transits after normalisation along with the transit model (black line).In the panel below, we show the residuals.

4. 3
.1.Case 1: Two-planet system and activity modelling trained on spectroscopy -uni-dimensional GP In the first case, we tested a multi-planet system model formed by TOI-5398 b and TOI-5398 c.PyORBIT simultaneously modelled the stellar activity and the signals of the two planets in the RV series.In the model, we included the inclination measured from photometry and the stellar mass M ⋆ as derived in Sect. 3 to determine the true masses of the planets.We then imposed a Gaussian prior on the eccentricities followingVan Eylen et al. (2019) and included the RV semi-amplitudes K b and K c in the RV modelling.Moreover, we used Gaussian priors on both orbital periods (P b , P c ) and central time of the first transits (T 0,b , T 0,c ) by considering the parameters outlined in Sect.4.2. .

Fig. 9 .
Fig. 9. Ground-based photometric data, simultaneously modelled with TESS transits.Different datasets are represented by distinct colours.Top: Individual transits of TOI-5398 b on the left and TOI-5398 c on the right.Bottom: Simultaneous transit of TOI-5398 b and c, observed during the night of 21 April 2022.

Table 4 .
Comparison of the three different

4. 4 .
Search for TTVs To investigate the potential presence of dynamical interactions between TOI-5398 b and c we performed a search for Transit Timing Variations (TTVs; e.g., Agol et al. 2005; Holman & Murray 2005; Borsato et al. 2014, 2019a, 2021) of planet b and c.

Fig. 10 .
Fig. 10.Phase-folded RV fit of TOI-5398 b planetary signal.The reported error bars include the jitter term, added in quadrature.The shaded area represents the RV model's ±1σ uncertainties.The bottom panel displays the residuals of the fit.
Fig. 12. O-C plot representing the observed (O) and calculated (C) transit times for the linear ephemeris of TOI-5398 b (see Table6).Each dataset is shown in a distinct colour.

Fig. 13 .
Fig. 12. O-C plot representing the observed (O) and calculated (C) transit times for the linear ephemeris of TOI-5398 b (see Table6).Each dataset is shown in a distinct colour.

Fig. 14 .
Fig. 14.Synthetic O-C diagrams (top-panel planet b, bottom-panel planet c) computed from the dynamical simulation with TRADES and with parameters from Table5.
Zhu et al. (2018) Notes.(a)Precise Radial Velocity.(b)Transit Time Variations.(c) Transmission Spectroscopy Metric of the giant planet in the system.

Fig. 16 .
Fig. 16.Mass-density distribution of all confirmed planets from the NASA Exoplanet Archive with mass and radius determination better than 20%.The red dots represent TOI-5398 b and c, while the remaining planets mentioned in Table 7 are represented by dots of different colours.We also plot theoretical mass-radius curves for planets of various pure compositions from Zeng et al. (2019, online at https:// lweb.cfa.harvard.edu/~lzeng/planetmodels.html): solid red indicates a pure-Iron core, brown an Earth-like rocky core (32.5% Fe and 67.5% MgSiO3), light-blue a 100% Water world at 1000K, and blue a 100% cold-Hydrogen world.

Fig. 17 .
Fig. 17.Mass-radius distribution of all confirmed planets with R p < 4 R ⊕ from the NASA Exoplanet Archive with mass and radius determination better than 20%.The red dot represents TOI-5398 c, the beige one TOI-2000 b, the green one TOI-1130 c, and the violet dots are respectively WASP-47 d and e.We include the Solar System planets Earth, Venus, Uranus, and Neptune.We add theoretical mass-radius curves from Zeng et al. (2019): solid brown indicates an Earth-like rocky core (32.5% Fe and 67.5% MgSiO3), beige a 100% water world at 1000K, and dotted beige a 1% hydrogen envelope and 99% Earth-like rocky core at the same temperature.Grey dashed curves represent densities ρ = 0.5, 1, and 5 g/cm 3 , respectively.

Fig. 18 .Fig. 19 .
Fig. 18.Ephemeris uncertainty of TOI-5398 b propagated until the start of 2030.On the x-axis, we show the date, while on the y-axis the error bars in minutes.In grey, we show the uncertainty propagation considering TESS data alone.In red, we show the improvement considering also Ground-based photometry data.

Fig. 20 .
Fig. 20.Radius evolution in time (x-axis) for various bulk metallicities (coloured lines, in units of planet b masses).The three plots adopt three different atmospheric heavy-element mass fractions.The light-blue dot represents TOI-5398 b.

Fig. 21 .Fig. 22 .
Fig. 21.Mass radius distribution of all confirmed planets that are present in the TEPCat catalogue (Southworth 2011).The orbital period of a planet is colour-coded when it has TSM > 90, M p > 0.1 M J , and R p > 8 R ⊕ .Planets without these characteristics are coloured grey.The dot size tracks the TSM.The planetary mass is shown on a logarithmic scale.

Table 1 .
Observations from TESS and HARPS-N summarised.

Table 3 .
Priors and outcomes of the model of planet b and c from the analysis of the photometric time series.In Appendix E, we show the second part of this Table, moved for improved readability of the main text.

Table 5 .
Priors and outcomes of the model of planet b and c from analysing spectroscopic series with a multidimensional GP framework (Case 3).

Table 7 .
Confirmed compact multi-planet systems, sorted by age.

Table 8 .
Atmospheric mass-loss simulations Mass loss rate (g/s) time scale a Mass (M ⊕ ) Radius (R ⊕ ) f atm (%) Mass loss rate (g/s)

Table 9 .
Retrieval results for the three different scenarios.