© The Authors 2024

1 Introduction

The characterisation of exoplanetary atmospheres is a key frontier in astrophysics. It offers valuable insights into the physical and chemical processes occurring in distant worlds. The field has witnessed significant advancements, with various observational techniques contributing to our understanding of their diverse compositions. Among these, high-resolution spectroscopy emerges as a powerful tool, since this technique allows one to resolve individual spectral lines, providing detailed information about the chemical composition and atmospheric dynamics. The ability to break degeneracies induced by broad spectral features in low-resolution spectra makes high-resolution spectroscopy indispensable for accurate atmospheric studies, which focus on the analysis of individual atomic lines and matching with model templates. High-resolution transmission spectroscopy has predominantly been applied to hot and massive exoplanets, such as hot Jupiters (HJs), because of their high equilibrium temperatures due to the small distances from their host star (Guilluy et al. 2022). The closeness and the consequently high equilibrium temperature, T_eq, favour photo-evaporation processes, and enhance the scale height defined as $$H=\frac{k_{\text{B}}{T}_{\text{eq}}}{\mu g_{\text{p}}},$$(1)

where k_B is the Boltzmann’s constant and μ is the mean molecular weight. In many cases, this leads to a large annulus that scatters the starlight during the transit, resulting in a strong signal. However, the large HJ’s mass increases the surface gravity, g_p, and hence diminishes the scale height. Therefore, the ideal targets for transmission spectroscopy are close and hot inflated planets with lower densities. The poorly explored class of warm, Saturn-like exoplanets combines these characteristics, presenting an exciting opportunity to broaden our understanding of planetary atmospheres beyond the better-studied HJs.

The TOI-5398 system (see Table 1 for the parameters) is a fascinating target for atmospheric studies due to its compact planetary architecture consisting of two planets orbiting close to the host star (Mantovan et al. 2024a), with the outer one 0.0980 au distance from the star. The outer planet, TOI-5398 b, is a warm, Saturn-like exoplanet with distinct features that make it an ideal candidate for transmission spectroscopy. It has an equilibrium temperature of 947 K, resulting in a scale height of H_b = 1106 km (with μ = 1.3 assuming a hydrogen-dominated atmosphere, instead of a hydrogen-and-helium-dominated atmosphere as in Fossati et al. 2022), which is larger than the mean value of an HJ (~450 km) and the largest transmission spectroscopy metric (TSM, Kempton et al. 2018) value (TSM = 288) among all known warm giant planets. The target’s position at the border of the Neptunian valley and savannah (Owen & Wu 2017; Armstrong et al. 2019) makes it even more fascinating, since its predicted evolution will lead it to the centre of the savannah in the coming hundreds of millions of years, as is shown in Fig. 1 and in Fig. E.1 from Mantovan et al. (2024b). The system evolution was strongly shaped by the closeness to the host star, a young (650 Myr) G-type star responsible for the photo-evaporation of the planets, as is discussed by Mantovan et al. (2024b). Based on their findings, planet b underwent minimal mass loss, and hence the present size and composition of this gas giant are probably primordial, while planet c, on the other hand, was more profoundly influenced by photo-evaporation, as is hinted at by its modal density.

The stellar X-ray and extreme ultraviolet (XUV) flux plays a crucial role in influencing the atmospheric properties of young exoplanets. Understanding the interactions between the XUV flux and the atmospheres of distant worlds provides critical insights into the atmospheric escape processes, photochemistry, and overall atmospheric evolution. The present study aims to explore the atmospheric composition of this intriguing exoplanet and contribute to the growing body of knowledge on diverse planetary systems. We aim to explore the composition of the upper layers of the atmospheres using different techniques involving both the visible and near-infrared (nIR) parts of our spectra, deepening our knowledge of the role of the XUV flux in the TOI-5398 system. For this purpose, we focus on lines arising from atomic species that are sensitive to photo-evaporation such as He I and Hα, and also search for other atomic species.

The paper is organised as follows. In Sect. 2, we describe the observational data. In Sect. 3, we illustrate the different methodologies employed to analyse the data in both the visible and nIR parts of the spectrum. In Sect. 4, we report the results we obtained via the single-line analysis and cross-correlation technique, while in Sect. 5 we evaluate the match between our results and the predicted photo-evaporation models. Eventually, we give a short summary that highlights the main conclusions in Sect. 6.

Fig. 1 Known exoplanets as a function of their radius and period from the NASA Exoplanet Archive (Akeson et al. 2013) within the Neptunian desert and savannah. The color represents the density of planets. TOI-5398 b is highlighted with a purple diamond, while planet c is highlighted with a cyan one. The arrow indicates the expected position of planet b within the savannah at the end of its evolution according to Mantovan et al. (2024b). The marker size is proportional to the planetary mass.

2 Observations

We observed the system TOI-5398 as part of the DDT proposal A46DDT4 (PI: G. Mantovan) on March 25 2023 UT with the GIARPS observing mode (Claudi et al. 2017) of the Telescopio Nazionale Galileo (TNG). GIARPS allows for simultaneous monitoring at high resolution in the visible band (0.39–0.69 μm) with HARPS-N (R ≈ 115000) and in the nIR (0.95–2.45μm) with GIANO-B (R ≈ 50000). The S/N as a function of the orbital phase for both the HARPS-N and GIANO-B datasets is shown in Fig. 2.

For the GIANO-B observations, we employed an ABAB nodding pattern (Claudi et al. 2017) to optimise thermal background noise subtraction. GIANO-B is an echelle spectrograph that covers four spectral bands in the nIR (Y, J, H, and K) divided into 50 orders. For our analysis, we focussed on order #39 in the Y band, where the He I triplet is located. The HARPS-N spectrum, on the other hand, is made up of 69 orders along its wavelength range. We focussed on orders #56 and #64 for the single-line analysis of the Na I doublet and Hα, respectively, while we employed the whole spectrum for the cross-correlation with template analysis. TOI-5398 b was scheduled for observations before, during, and after the planetary transit. Coincidentally, during the selected transit window, TOI-5398 c was also transiting (see Fig. 3). The impact of planet c’s transit on the analysis of planet b is further discussed in the next sections.

We observed the target with an exposure time of 600 s for HARPS-N and 300 s for GIANO-B, collecting 47 and 82 spectra, with an average signal-to-noise ratio (S/N) of 43 and 33 for HARPS-N and GIANO-B, respectively. A summary log of the observations is provided in Table 2. The same dataset was used in Mantovan et al. (2024b) to analyse the geometry of the system and evaluate the impact of photo-evaporation on the two planets of the system.

Fig. 2 S/N as a function of the orbital phase for both the HARPS-N and GIANO-B datasets. The air mass is colour-coded. The dashed lines represent the ingress and egress of TOI-5398 b.

3 Methods

3.1 Extraction of single-line transmission spectra

The GIARPS observing mode of the TNG is crucial for simultaneously observing both the visible and nIR parts of the planetary spectrum. Recent studies by Guilluy et al. (2024) have demonstrated that combining data from HARPS-N and GIANO-B enables the retrieval of valuable information about the atmospheres of exoplanetary targets. In particular, the joint analysis of the He I triplet and Hα plays an important role in understanding if a detection may be affected by stellar activity (e.g., Guilluy et al. 2020, 2023). The two different datasets require different frameworks to extract the planetary signal around the studied lines. We discuss the different methodologies to extract the transmission spectra in the next subsections. Eventually, we apply a uniform analysis on the extracted spectra.

We have not included in our analyses any correction to take into account the presence of planet c. Despite planet c being closer, and hence hotter than planet b, the inner planet is almost five times denser than its outer companion, resulting in a smaller scale height (~880 km vs. ~1100 km). Therefore, we think that planet c does not play a major role in our analysis, as we deeply discuss in Sect. 5.

The main difference between the nIR and visible datasets concerns the modelling and removal of the Rossiter–McLaughlin effect (RME) and the centre-to-limb variation (CLV). In fact, we neglect the impact of the RME and CLV on the Helium triplet, as several studies have shown that these effects have little (or no) impact on the helium lines (e.g., Nortmann et al. 2018; Salz et al. 2018; Allart et al. 2019; Fossati et al. 2022; Allart et al. 2023; Guilluy et al. 2023) As was stated before, we do not correct for planet c RME+CLV in the visible dataset as the RME+CLV impact depends on the size of the obscured stellar disc and planet c has a radius three times smaller than planet b. Therefore, the surface covered by planet c is almost ten times smaller than the one obscured by planet b. This means that the intensity of the RME effect should be ten times smaller than planet b’s one. Moreover, Mantovan et al. (2024b) discussed the same problem regarding their Rossiter analysis to retrieve the orbital obliquity and showed that the radial velocity RME of planet c is comparable with the stellar activity.

Fig. 3 Radial velocity profile of the two planets observed as part of the DDT proposal A46DDT4. The transit of planet c completely overlaps with planet b one.

3.1.1 GIANO-B nIR region

We employed the G0FI0 pipeline (Rainer et al. 2018) to extract spectra from the raw GIANO-B images. This process involved dark subtraction, flat-field correction, and the removal of bad pixels. Additionally, we performed a preliminary wavelength calibration using a U-Ne lamp spectrum as a template in the vacuum wavelength frame. We focussed on the msld spectra, with the echelle orders separated. To enhance the initial wavelength solution, we employed the same approach described in our previous works (e.g., Giacobbe et al. 2021; Guilluy et al. 2022). We aligned all the spectra to the Earth’s atmospheric rest frame, assuming it as the frame of the observer (disregarding any ~10 m s⁻¹ differences due to winds), by measuring any shifts relative to an average spectrum taken as a template for the night. Subsequently, we refined the wavelength solution by utilizing an Earth’s atmospheric transmission spectrum generated using the Sky Model Calculator¹. For the remaining analysis, our focus was solely on order #39, which includes the He I triplet. The magnitude of the wavelength calibration refinements is ~0.56 km s⁻¹, approximately one fifth of a pixel.

To account for the Earth’s atmosphere contamination, we utilised the ESO software Molecfit². (Smette et al. 2015; Kausch et al. 2015) to correct for the transmission telluric lines, with a particular focus on the H₂O lines at approximately 1083.51 nm. The GIANO-B nodding acquisition mode automatically corrects for OH emission lines. However, in our previous works (e.g., Guilluy et al. 2023, 2024), we observed that the A-B subtraction could leave residuals at the position of the deepest OH doublet (around 1083.43 nm) due to variations in atmospheric seeing during observations. To address this, we masked any potential residuals in the final transmission spectra, as is described in Guilluy et al. (2023).

To separate the stellar contribution from the potential planetary signal, we performed transmission spectroscopy on the GIANO-B spectra according to the following approach. To start with, we moved the spectra to the star’s rest frame and then we normalised the spectra to unity by median division, excluding the spectral region around the He I triplet. Following that, we generated a master stellar spectrum, S_out(λ), by averaging the out-of-transit spectra (i.e. with an orbital phase smaller than t₁ or greater than t₄, the two ingress and egress contact points). We performed a visual comparison between S_out(λ) and a master-in spectrum, which was derived by averaging the in-transit between t₂ and t₃ (that is, the contact points of total transit). An absorption feature is readily discernible in the stellar spectrum, precisely coinciding with the position of the He I triplet (Fig. 4). To better investigate this feature, we derived individual transmission spectra, T(λ,i), by dividing each spectrum by S_out(λ). Finally, we linearly interpolated the transmission spectra in the planet’s rest frame. The left panel of Fig. 5 displays the 2D transmission spectroscopy map, while the right panel shows the He I spectroscopic light curve computed within a passband of 0.075 nm centred on the peak of excess absorption in the planet rest frame (Allart et al. 2019). To quantify the contrast, c, of the extra absorption observed at the position of the He I triplet, we employed the differential evolution (DE) Markov chain Monte Carlo (MCMC) method (Ter Braak 2006). This involved fitting a Gaussian profile to the mean in-transit transmission spectrum, computed by averaging the 2D maps between the transit contact points, t₂ and t₃. We fitted the peak position, the full width at half maximum (FWHM), the contrast value (c), and an offset for the continuum. Correlated noise in the transmission spectrum was accounted for using Gaussian process (GP) regression within the same DEMCMC framework, employing a covariance matrix described by a squared exponential kernel (Guilluy et al. 2024; Sicilia et al. 2024). Additionally, uncorrelated noise was taken into account through a jitter term, σ_j. The correction with the GP regression model is illustrated in Fig. 6, with detailed posterior distributions provided in Fig. A.3. The best-fit parameters obtained from the DE-MCMC Gaussian analysis are summarised in Table 3. Parameter values and their 1σ uncertainties were determined from the medians and the 16%-84% quantiles of their posterior distributions.

Table 3 also reports the effective He I radius (e.g., Chen et al. 2018) that would produce the observed absorption contrast, c. We then normalised it to the atmospheric H_eq to compute the quantity, ${\delta }_{{\text{R}}_{\text{p}}}/H_{\text{eq}}$ (Nortmann et al. 2018), which represents the number of scale heights probed by the atmosphere in the spectral range under consideration. Here, H_eq was computed by assuming a mean molecular weight of 1.3. (Fossati et al. 2022).

Fig. 4 Master-out (black) and master-in (red) spectra in the star frame. The dashed grey line indicates the spectrum before the Molecfit correction, while vertical dashed black and blue lines denote the position of the stellar He I triplet and OH emission line, respectively. In-transit absorption is visible by eye, especially in the zoomed-in panel.

Fig. 5 Analysis results for the He I triplet. Left panel: 2D maps of transmission spectra in the planet rest frame in the region of the He I triplet, as a function of wavelength and planetary orbital phase. The contact points, t1, t2, t3, and t4, are marked with horizontal blue lines. The regions affected by OH contamination are masked. Some residuals are left at the position of the Si line (~1083 nm). This is due to the depth of the line, which can cause difficulties in spectral extraction (see, e.g., Krishnamurthy et al. 2023). Dashed red lines mark the position of the He I lines. Right panel: light curve in the planet rest frame computed in a bandpass of 0.075 nm (equivalent to ~20 km s−1). Black points were computed with a phase bin of 0.002. Vertical dashed black lines indicate the position of the transit contact points, t1, t2, t3, and t4.

3.1.2 HARPS-N visible region

We obtained the transmission spectra of the individual lines using a method similar to the one outlined in Wyttenbach et al. (2015), and thus compared spectra during and outside of transit. The raw data had already been processed using version 3.7 of the HARPS-N Data Reduction Software (DRS; Pepe et al. 2002). At this point, the ID spectra we employed included the stellar signal, the planetary signal, and telluric contamination, all presented in the solar system barycentric reference frame with the wavelength information provided in the air reference frame.

Correcting for telluric absorption is crucial, especially in the red region of the HARPS-N spectra, where the Earth’s atmosphere shows significant H₂O and O₂ absorption features. Figure A.1 illustrates an example of removing telluric lines near the Hα and Na I doublet regions. The software Molecfit provides telluric-corrected spectra in the solar barycentric reference frame. Consequently, we adjusted the spectra to the stellar reference frame by accounting for the star’s systemic radial velocity, v_sys. The telluric-corrected spectra were then normalised to a common continuum level within a narrow wavelength range around each absorption line of interest. Before transitioning to the planetary reference frame, we must consider the combined effects of the RME (Rossiter 1924; McLaughlin 1924) and CLV, which result from stellar rotation and limb darkening, respectively. We addressed the RME+CLV effects by employing a similar approach to the one outlined by Yan et al. (2017), utilising stellar models. Specifically, we used ATLAS9 stellar models (Kurucz 1992, 2005, 2014, 2017) to generate the disc-integrated stellar model based on the system parameters listed in Table 1. We computed spectra for both a non-rotating star and a rotating one. To determine the portion of the stellar disc obscured by the planet at each phase (referred to as the obscured region), we used PyLightCurve (Tsiaras et al. 2016) to calculate the planet’s transit path. We then created a grid with a pixel size of 0.01 R_* × 0.01 R_* to approximate the stellar disc and calculated the stellar flux for the obscured region at each phase. The intensity of each pixel was determined based on the limb darkening angle and the radial velocity shift relative to the stellar rotation axis. Using the ATLAS9 models, we computed the spectrum of the obscured region by summing the spectra of each pixel, each appropriately shifted by its own radial velocity. We derived the limb darkening coefficients for our stellar model using the ExoTethys (Morello et al. 2020) function Sail with a quadratic limb darkening law, consistent with ATLAS9. The integrated flux, F_λ, corresponds to the spectrum of the non-rotating model. Rotational broadening was obtained by summing the contributions of each pixel, each shifted by its own radial velocity. The observed master-out spectrum is broadened by stellar rotation, so we used the rotationally broadened model to normalise the continuum of the non-rotating model. Finally, we subtracted the modelled obscured region from the master-out spectra and divided all spectra by this result. Once the spectra in the stellar reference frame had been corrected, we shifted all the in-transit spectra to the planet reference frame by correcting for its radial velocity profile in a circular orbit scenario using the parameters listed in Table 1. The final transmission spectrum of each line was calculated as the error-weighted average of all the full in-transit spectra (t₂–t₃) in the planet reference frame. We then performed a Gaussian fit, returning the absorption depth, the FWHM, and the velocity shift with respect to the line wavelength corresponding to a zero radial velocity. In the beginning, we used the planetary radius to represent the size of the obscured region of the stellar disc. We then translated the line absorption depths, returned by the Gaussian fit, into planetary radii and we used these values to repeat the RME + CLV removal with an increased size that accounts for both the planetary radius and the atmospheric extension. We repeated the loop, retrieving a new depth expressed in planetary radii, until the threshold convergence (0.001 Rp) was reached. Once we found the optimal value for the planet+atmosphere region of the stellar disc obscured, we used the same DE-MCMC+GP tool employed in the GIANO-B analysis to fit our final corrected transmission spectra. Figure 6 shows the correction with the GP regression model, while the posterior distributions are provided in Fig. A.3. The best-fit parameters obtained from the DE-MCMC Gaussian analysis are reported in Table 3. As for the He I contrast, we derived the corresponding effective radius for both Hα, and the Na I doublet. As the correlated noise is weaker in HARPS-N spectra compared to the GIANO-B ones, we assessed the necessity of the GP regression model correction using the Bayesian information criterion (BIC, Kass & Raftery 1995). We considered the Gaussian+GP model to be the favourite compared to the simple Gaussian model (which also considered uncorrelated jitter) only if accompanied by a lower BIC and a difference in BIC (ΔBIC) > 10, along with a Bayesian evidence >150. For all the three investigated lines (i.e. Hα, NaD1, and NaD2), the Gaussian+GP model emerged as the preferred choice. We provided the ΔBIC and the corresponding Bayesian evidence in Table A.1. We computed the spectroscopic light curves similarly to the He I ones, considering a width of 20 km s⁻¹.

Fig. 6 Transmission spectra of He I triplet and Hα and Na I doublet after GP correction. For each investigated line, the GP correction is shown. Left panel: transmission spectrum centred on the line (in the planet rest frame) with the GP regression model overplotted, along with the lσ uncertainty intervals (in blue), and the Gaussian+GP model (in red). Right panel: final transmission spectrum after removing the GP model. Vertical dotted black lines indicate the position of the investigated line. The error intervals for the Gaussian fit were computed by displaying 1000 Gaussian fits within the lσ uncertainties of the derived parameters, spanning the 16%–84% quantiles.

3.2 Cross-correlation with templates

To optimise the detection of certain elements with hundreds of weak lines, the cross-correlation function (CCF) is more suitable than transmission spectroscopy analysis (Brogi et al. 2012; Hoeijmakers et al. 2019; Borsa et al. 2021), especially in the case of individual nights, since it allows one to boost the S/N. We used the CCF method in the visible band, searching for any exoplanetary signal. Before applying the CCF technique to our HARPS-N data, we employed the same telluric profile retrieved from the HARPS-N s1d spectra using Molecfit, and we used it to correct the HARPS-N e2ds orders from telluric lines for each exposure (Hoeijmakers et al. 2020). To improve the accuracy of the CCF analysis, we chose to omit the initial five orders and the last one from all exposures due to their low S/N. We constructed the mono-dimensional spectra from normalised telluric-corrected e2ds orders. These 1D spectra were then cross-correlated with the templates at 2000K, 2500K, and 3000K provided by Kitzmann et al. (2023) after convolving them with the HARPS-N resolution (R ~ 115000). Since these templates were derived using the Sun’s radius as a reference, before any CCF computation, we multiplied them by the ratio (R_⊙/R_*)², as is stated in Kitzmann et al. (2023), using R_*, from Table 1. The templates were shifted in radial velocities between [−200 km s⁻¹, +200 km s⁻¹ ] with a step of 1 km s⁻¹. The resulting CCF was then shifted in the stellar reference frame, interpolated in radial velocity between [−150 km ^s−1, +150 km s⁻¹] maintaining the 1 km s⁻¹ step, and divided by the mean of the out-of-transit (master-out). To achieve a more precise normalisation and address any imperfections, we divided the CCF by the median values across the exposures. To remove any residual fluctuations, we performed a Fourier transform, discarding all the frequencies beyond 100 km s⁻¹.

In this way, we investigated the presence of all the species given by Kitzmann et al. (2023), including both neutral and ionised elements, with at least one absorption line in the visible band covered by our order selection. The resulting CCFs for some elements exhibit the characteristic Doppler shadow. We removed it following the Rainer et al. (2021) approach: we shifted the data in the Doppler shadow reference frame and fitted each column in a region where the Doppler shadow fell with a fifth-degree polynomial. However, the specific geometry of the TOI-5398 system is such that the Doppler shadow and the planetary trace overlap (see Casasayas-Barris et al. 2022; Sicilia et al. 2024, for a similar scenario). Therefore, in this scenario, the Doppler shadow signal may completely mask the potential planetary signal, and even if the Doppler shadow is modelled and removed, the potential planetary signal may also be eliminated. Figure A.4 shows the CCF obtained with the Na I, Ca I, and Fe I templates at 2000K. The Doppler shadow for these elements was modelled out for CaI and FeI.

4 Results

4.1 Single line and cross-correlation analyses

By comparing the out-of-transit spectra with the in-transit ones, we were able to discern clear He I, Hα, NaD1, and NaD2 absorption signals. Applying the framework described in Sect. 3, we estimated contrasts of the excess absorption of ${3.57}_{-0.50}^{+0.51}\%$ (6.9σ), ${1.39}_{-0.28}^{+0.28}\%$ (5.0σ), ${3.79}_{-0.47}^{+0.50}\%$ (7.9σ), and ${3.36}_{-044}^{+0.46}\%$ (7.4σ) for He I , Hα, NaD1, and NaD2, respectively. These values correspond to effective planetary radii of ~2.33 Rp, ~1.65 Rp, ~3.81 Rp, and ~2.28 Rp, the extensions of which are well below the planet’s Roche lobe radius (~5.8Rp), computed using Eq. (2) from Eggleton (1983). This indicates that the planetary atmosphere is not escaping due to the gravitational pull of the star, as is discussed in the next section. Assuming a scale height (H) of ~1106 km, and considering 1σ uncertainties, we found that the He I atmosphere probes a number of scale heights varying between ~70 and ~88 H, while the Hα atmosphere varies between 32 and 46 H, NaD1 between 70 and 96 H, and NaD2 between 64 and 86 H. We found the He I absorption signal to show a clear blueshift, ~8 km s⁻¹, while the optical species do not display either a significant blueshift or a redshift. He I and Ha have similar FWHMs, larger than the Na I doublet ones. For the He I and Hα lines, our detections seem to be supported by the spectroscopic light curves, which hint at the presence of an He I tail that we discuss later. The Na doublet light curves do not show a clear transit signal, despite showing a decrease in flux in the second half of the transit. This is possibly related to the overlapping RME, as is discussed later. We also investigated the presence of other metal species in the atmosphere, analysing a wide forest of lines belonging to Mg I, Ca I, Fe I, and Fe II, but unfortunately we did not detect any significant line. The cross-correlation of the extracted spectrum with Kitzmann et al. (2023) templates did not return any detection of the investigated species. As was already stated, this framework is particularly challenging in our case due to the overlap of the possible signals with the Doppler shadow. This also affected the extraction of single lines in the visible spectrum of TOI-5398 b. An under-correction of RME may cause one to miss the signal, while an over-correction may lead to a false detection. We compared the results obtained applying the method described in Sect. 3.1.2 (already employed in Guilluy et al. 2024 and D’Arpa et al. 2024) with the results obtained using SLOPpy (Sicilia et al. 2022), which employs a different method to model the RME+CLV, as is described in B.

4.2 Comparison with ATES predicted spectra

To support the interpretation of our He I observations, we conducted simulations of the atmospheres of the two planets and their respective He I absorptions using the 1D hydrophotoionisation self-consistent code ATES (Caldiroli et al. 2021) and the new transmission probability module (TPM, Biassoni et al. 2024). The ATES code derives values pertaining to the 1D atmospheric radial profile such as temperature, pressure, density, and velocity, along with the mass loss rate, M.

For these simulations, we used the planetary parameters listed in Table 1 for planets b and c, while the stellar spectrum was calculated as follows. For the UV component, we adopted a stellar spectrum from the Phoenix library (Husser et al. 2013), to which we added the Lyman alpha emission following the prescriptions of Linsky et al. (2020) and in the manner described in Locci et al. (2022). For the EUV and X-ray components, we produced a synthetic spectrum from a composite chromosphere-corona plasma emission measure distribution versus temperature, EMD(T), constructed as follows: the coronal part (log T > 5.5 K) was interpolated from the grid of EMDs provided by Wood et al. (2018) for a star with a surface X-ray flux of Fx = 1.5 × 10⁶ erg s⁻¹ cm⁻²(corresponding to Lx = 1 × 10²⁹ erg s⁻¹), while the chromospheric section (log T from 4.0 to 5.5 K) was approximated as a power law following Sanz-Forcada et al. (2011), with the constraint that the Feuv/Fx ratio follows the scaling law indicated by Johnstone et al. (2021). Then, the optically thin XUV spectrum (1.24–1700 A) was computed adopting the plasma emissivities from the CHIANTI v7.13 atomic database as in Maggio et al. (2023).

We started our simulations for both planets assuming a constant He/H number abundance of 1/12 alongside the overall atmosphere. This provided us M of 10^11,70 g s⁻¹and 10^11,44 g s⁻¹for planets b and c, respectively, indicating that nowadays the innermost planet, c, is affected less than planet b by the atmospheric escape effect. However, this does not exclude the possibility that in the past planet c may have experienced significantly greater atmospheric loss. We investigated the He I absorption profiles for both planets using the TPM module. For this analysis, we firstly simulated the atmospheric profiles changing the He/H number abundance, using values of 0.01, 0.02, 0.04, and 1/12. We then repeated the simulations, varying the XUV flux to include the enhanced stellar activity of a young star. TOI-5398 is a G-type star with an age comparable to the one of the Hyades cluster stars. According to Fig. 4 of Preibisch & Feigelson (2005), in such stars, the typical X ray luminosity is centred at Lx = 1 × 10²⁹ erg s⁻¹, which is the value we adopted, with a variation from the pre-main-sequence through the main sequence that is almost a factor of two. Therefore, we repeated the same simulations described above, varying the XUV luminosity by a factor of two. We used the results as a function of the He/H number abundance to estimate the uncertainties of the values retrieved with the nominal luminosity. The resulting He I absorption profiles obtained with these different He/H number abundances for planet c are always below the observed signal, while for planet b, the simulations are consistent with the observed signal only with a He/H number abundance between 0.01 and 0.02 (see Fig. 7).

Fig. 7 He I simulated absorption profiles for both planets changing the He/H number abundance. The green line corresponds to the He I observed signal with lσ error bars given in Table 3. The blue and red bands correspond to the simulations obtained for the two planets, varying the X-ray luminosity by a factor of two.

5 Discussion

Within the single-line analysis context, the He I triplet in the nIR, and the Hα and Na I doublet in the visible, are among the pivotal lines studied in exoplanetary atmospheres, thanks to their strength and to their role as markers of features such as photo-evaporation and stellar activity. The Na I doublet has been detected in several hot exoplanets with temperatures ranging from 963 K in WASP-69 b (e.g., Casasayas-Barris et al. 2017; Khalafinejad et al. 2021) to 3921 K in KELT-9 b (e.g., Langeveld et al. 2022; D’Arpa et al. 2024). Therefore, according to our detection, TOI-5398 b will be among the coldest (947 K) planets showing the Na I doublet, to our knowledge. The He I and Ha distribution has recently been discussed by Orell-Miquel et al. (2024), who explored the high-resolution spectroscopy observations from CARMENES and GIARPS checking for He I and Hα signals in 20 exoplanetary atmospheres. From their work, we can see that the joint detection of He I and Hα has been observed only in two targets: WASP-52 b (Chen et al. 2020; Kirk et al. 2022; Vissapragada et al. 2020) and HAT-P-32 b (Czesla et al. 2022; Zhang et al. 2023), with WASP-52 b being the only one with a detection of all three species detected in our study.

The height distribution of the species in our atmosphere suggests that He I and Na I lie in an external atmospheric layer with respect to Hα. Hα forms through absorption from the n=2 level of neutral hydrogen, which is likely excited by the strong stellar Lyα line and/or by the photoionizing radiation of our young star followed by subsequent recombination and population of the n = 2 level as in solar-like stellar atmospheres (Cram & Mullan 1985). To justify the observed Hα absorption, it is necessary for hydrogen to have a sufficient number of atoms in the excited level n=2. If we go too high in the planetary atmosphere, hydrogen finds itself in a less dense medium and has enough time to decay to the fundamental level n=1 before being able to absorb a photon in the Hα line. In the case of the He I triplet, the situation is different because the lines form through absorption from a metastable level. Even if we are high in the atmosphere, the level decays over much longer timescales than those required for photon absorption in the lines. Moreover, in a less dense medium, the probability that the atom decays to the fundamental level due to collision is low, facilitating the process that produces the He I triplet.

Indeed, absorption of He I at 1083 nm is seen in certain cases, even at a great distance from the planet, along almost its entire orbit (see Zhang et al. 2023, for example). The sodium retrieved height is more challenging to explain, since the sodium atom has a low ionisation potential, and therefore, far from the planet, it risks being ionised. However, the absorption of the two doublet lines starts from the fundamental level, so the atom does not need a dense surrounding medium to bring it to an excited level through thermal collisions to produce the line. Therefore, if it does not ionise too much, some absorption in the Na I doublet at a certain distance from the planet might be observable. In this case, we stress the mismatch between the results obtained with SLOPpy, where the Na I doublet is found at 1.12 Rp. Since the species is not detected via cross-correlation with templates, the signal we detect may be due to a spurious signal arising from a residual of the RME+CLV contamination overlapping with the planetary signal. This hypothesis could explain the fact that the Na doublet lines do not show a clear transit signal in the light curve despite showing a flux decrease in the second half of the transit. Future observations and improvements in the RME+CLV modelling will be necessary to clarify this point.

The distribution of Hα and Na I resemble the same behaviour observed in WASP-52 b, where Chen et al. (2020) find both sodium lines lying at an effective radius larger than the Hα one. The extended atmosphere observed for our target may be explained taking into account the photo-evaporation induced by the host star. In fact, planets close to their host star experience a strong irradiation that causes the exoplanetary atmosphere to heat up, and hence to expand. The larger thermal energy induces an increase in velocity, making the particles belonging to the atmospheric species reach the Roche lobe and escape from the planetary gravitational well. This is particularly relevant for younger stars, which generally have a larger XUV flux such as TOI-5398. Mantovan et al. (2024b) evaluated photo-evaporation in the system during a period spanning from 3 Myr to 650 Myr. The two planets experienced a different evolutionary path, with planet b being less influenced by atmospheric escape, while planet c, on the other hand, experienced the loss of almost 60% of its mass, resulting in a less extended atmosphere nowadays. At this point of its evolutionary path, planet c should show a less extended atmosphere, while planet b remains dominated by an extended envelope despite experiencing a currently larger mass loss. The simulations we performed with ATES retrieve M of 10^11,70 g s⁻¹ and 10^11,44 g s⁻¹ for planets b and c, respectively, confirming that part of the outer layers of the atmosphere are currently undergoing atmospheric escape. This scenario is ideal for transmission spectroscopy due to the larger extension of the scattering annulus surrounding the planet.

From our results, we observe that TOI-5398 b shows an extended atmosphere that reaches ~2.4 Rp, still well under the Roche lobe at ~5.8 Rp, computed using Eq. (2) from Eggleton (1983). This does not question the photo-evaporation ongoing for the target, since we obtained the height from a geometrical point of view from the absorption depth of the lines detected. Since the latter is related to the point where the optical depth is one, the height does not represent an individual layer where the species lie, since it is distributed according to the density profile. Therefore, we are able to find particles overcoming the Roche lobe despite observing an absorption that peaks well below.

The past and current status of the two planets, along with the ATES models, should help us to address the possible contamination of planet c during the transit that we observed. In Fig. 3, we show that the two transits almost completely overlap and only three points of planet b are unaffected by planet c. While it is impossible to completely exclude the uncertain contribution of planet c, its role in the observed signal seems to be secondary, as is jointly supported by the two diagnostics: planet c has a less extended atmosphere due to its past evaporation and the spectra retrieved by ATES (Fig. 7) confirm our suggestions, since the He I arising from the planet c is not able to justify the absorption we observe, even considering two times the nominal XUV flux. In particular, in order to match the He I triplet retrieved absorption, we should increase the XUV flux by more than ten times its value. On the other hand, planet b is expected to show hints of photo-evaporation and the ATES spectra manage to mimic the 3.57% absorption with a He/H number fraction ranging between 0.01 and 0.02. While retrieving the He/H number fraction is not the main aim of our comparison with the TSM module, we are able not only to evaluate the contribution of the two planets but also to discern a range of He/H number fraction as a function of the XUV flux adopted, which in a young star may vary. The use of models including the Hα line may be able to break this degeneracy.

Observing the spectroscopic He I light curve in Fig. 5, we note that the signal seems to continue beyond T4. This may be an indication of the presence of an extended atmosphere or some sort of cometary tail trailing the planet. The tail may be further proof of planet b being photo-evaporating. As further evidence, we have absorption of He I at 8 km s⁻¹, which could originate from the cometary tail driven by stellar winds. However, the tail does not seem to be present in the other species' light curve, as it is in the He I one. This may arise from the lower number of points in the visible band due to a larger exposure time.

6 Conclusions

In this paper, we have performed a high-resolution nIR+optical transmission spectroscopy study of the Saturn-like planet TOI-5398 b, analysing one transit collected on March 25 2023 with the GIARPS observing mode at the TNG. The same dataset has been used before by Mantovan et al. (2024b) to constrain the orbital parameters and evaluate the evolutionary path of the system.

We detected the presence of metastable He I with the GIANO-B spectrograph and of a Hα and Na I doublet with HARPS-N employing single-line analysis. According to the ExoAtmospheres database³, this represents one of only a few simultaneous detections of He I and an Hα and Na I doublet in the atmosphere of an exoplanet.

The cross-correlation with templates did not return any detection, probably due to the overlapping RME, the low temperature of the planet, and the exiguous dataset. In fact, the difficulty of extracting the planetary signal in the atmosphere of smaller planets is enhanced by the fact that we observed only one transit. As is discussed in Guilluy et al. (2024), the study of individual lines across different nights is crucial to identify possible contamination due to stellar activity, and hence more transits will be necessary to further corroborate our detections.

Despite us attributing the detected signal to planet b according to the evolutionary history of the system and to the ATES simulations, future individual observations of both planets will also be needed to definitively clarify the origin of the signals detected. Future observations with longer baselines would also help in investigating the possible presence of a cometary tail, along with 3D hydrodynamical simulations such as the ones performed in Matsakos et al. (2015) and Colombo et al. (2024).

Acknowledgements

We sincerely thank the Referee for their comments stressing some important factors and analyses that we did not consider initially, and that we think improved the quality of the manuscript. The authors acknowledge financial contribution from PRIN INAF 2019 and from the European Union – Next Generation EU RRF M4C2 1.1 PRIN MUR 2022 project 2022CERJ49 (ESPLORA), PRIN MUR 2022 project PM4JLH ("Know your little neighbours: characterizing low-mass stars and planets in the Solar neighbourhood") and PRIN MUR 2022 (project No. 2022J7ZFRA, EXO-CASH). The authors acknowledge the support of the ASI-INAF agreement 2021-5-HH.0 and the project “INAF-Astrofisica Fondamentale GAPS2” and INAF GO Large Grant 2023 GAPS-2. GMa acknowledges support from CHEOPS ASI-INAF agreement no. 2019-29-HH.0. Part of the research activities described in this paper were carried out with contribution of the Next Generation EU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures, Project IR0000034 – “STILES – Strengthening the Italian Leadership in ELT and SKA”. The authors acknowledge the Italian center for Astronomical Archives (IA2, https://www.ia2.inaf.it), part of the Italian National Institute for Astrophysics (INAF), for providing technical assistance, services and supporting activities of the GAPS collaboration.

Appendix A Additional figures and tables

Fig. A.1 Example of telluric correction for the Hα and Na I doublet regions. After normalisation, the DRS processed data (in orange) was telluric corrected with Molecf it using the telluric model in blue, here shifted upwards for clarity.

Fig. A.2 Same as Fig. 5 but for the lines detected in the visible range. The dashed black lines represent the wavelength position of the line.

Fig. A.3 Posterior distribution of the investigated parameters in the DE-MCMC analysis for the four investigated lines. The excess of absorption c [%], offset [%], peak position, and FWHM correspond to the parameters we used in the Gaussian fit, while the jitter term σj, the semi-amplitude of the correlated noise h, the correlation length λ were used to parametrise the SE kernel within the GP.

Fig. A.4 Cross-correlation and Kp-Vsys maps of Na I, Ca I, and Fe I obtained with the Kitzmann et al. (2023) templates at 2000K. In the first column the CCF are depicted in terms of ppm. The dotted horizontal lines indicate the T1 and T4 contact points, while the red dashed line is the planetary trace. The Doppler shadow follows the white slanted line. Second columns are the Kp-Vsys maps. The planetary signal is expected to be limited between the green pointer in the middle of the figures, while the red lines point at the minimum value of the Kp-Vsys between Kp 75 km s−1, 125 km s−1 and RV -50 km s−1, +50 km s−1. The last column is the Kp-Vsys 1D map evaluated at the Kp minimum value. The grey region is the 3σ edge, while the vertical navy dashed line refers to the RV 0 km s−1.

Appendix B Comparison with SLOPpy

We compared the results obtained in the visible range with the method described in this work with SLOPpy to evaluate the impact of different correction of the RME+CLV. The two methods differ by many factors but rely on the same idea of correcting for the part of the stellar disc obscured by the planet. Our method sets the extension of the atmosphere increasing iteratively the radius retrieved from the Gaussian fit until we reach the convergence of 0.001 Rp, as described in Sect. 3.1.2 and hence we retrieve the R_eff from the absorption depth. On the other hand, SLOPpy fits the size of the atmosphere by comparing the 2D map containing the modelled RME+CLV and the planetary signal employing an MCMC approach. We stress that the extension retrieved by SLOPpy is not the physical extension of the atmosphere but a corrective factor. Therefore, the major difference in the two frameworks is the extension adopted to perform the RME + CLV correction. We performed the comparison in the following way: we employed SLOPpy which returns the extracted transmission spectrum and fit it using an MCMC method to find the best Gaussian fit values and the best corrective radius factor to employ in the RME+CLV correction. Then, we applied the same GP+DEMCMC analysis described in Sect. 3 on the spectrum extracted with SLOPpy. The results are shown in Fig. B and Table B.1. The FWHM and the velocity shift are within 1σ for both the Ha and Na I doublet, while, for the latter, the contrasts we retrieve are larger than SLOPpy’s ones. A deep comparison between the two frameworks is out of the scope of this work and hence we limit to say that, despite the difference in depth, in both cases we detect the presence of the two Na lines anf the Hα in the atmosphere of planet b.

Fig. B.1 Same as Fig. A.3 but for SLOPpy.

References

All Tables

All Figures

Fig. 1 Known exoplanets as a function of their radius and period from the NASA Exoplanet Archive (Akeson et al. 2013) within the Neptunian desert and savannah. The color represents the density of planets. TOI-5398 b is highlighted with a purple diamond, while planet c is highlighted with a cyan one. The arrow indicates the expected position of planet b within the savannah at the end of its evolution according to Mantovan et al. (2024b). The marker size is proportional to the planetary mass.

In the text

	Fig. 2 S/N as a function of the orbital phase for both the HARPS-N and GIANO-B datasets. The air mass is colour-coded. The dashed lines represent the ingress and egress of TOI-5398 b.
In the text

	Fig. 3 Radial velocity profile of the two planets observed as part of the DDT proposal A46DDT4. The transit of planet c completely overlaps with planet b one.
In the text

	Fig. 4 Master-out (black) and master-in (red) spectra in the star frame. The dashed grey line indicates the spectrum before the Molecfit correction, while vertical dashed black and blue lines denote the position of the stellar He I triplet and OH emission line, respectively. In-transit absorption is visible by eye, especially in the zoomed-in panel.
In the text

Fig. 5 Analysis results for the He I triplet. Left panel: 2D maps of transmission spectra in the planet rest frame in the region of the He I triplet, as a function of wavelength and planetary orbital phase. The contact points, t1, t2, t3, and t4, are marked with horizontal blue lines. The regions affected by OH contamination are masked. Some residuals are left at the position of the Si line (~1083 nm). This is due to the depth of the line, which can cause difficulties in spectral extraction (see, e.g., Krishnamurthy et al. 2023). Dashed red lines mark the position of the He I lines. Right panel: light curve in the planet rest frame computed in a bandpass of 0.075 nm (equivalent to ~20 km s−1). Black points were computed with a phase bin of 0.002. Vertical dashed black lines indicate the position of the transit contact points, t1, t2, t3, and t4.

In the text

Fig. 6 Transmission spectra of He I triplet and Hα and Na I doublet after GP correction. For each investigated line, the GP correction is shown. Left panel: transmission spectrum centred on the line (in the planet rest frame) with the GP regression model overplotted, along with the lσ uncertainty intervals (in blue), and the Gaussian+GP model (in red). Right panel: final transmission spectrum after removing the GP model. Vertical dotted black lines indicate the position of the investigated line. The error intervals for the Gaussian fit were computed by displaying 1000 Gaussian fits within the lσ uncertainties of the derived parameters, spanning the 16%–84% quantiles.

In the text

	Fig. 7 He I simulated absorption profiles for both planets changing the He/H number abundance. The green line corresponds to the He I observed signal with lσ error bars given in Table 3. The blue and red bands correspond to the simulations obtained for the two planets, varying the X-ray luminosity by a factor of two.
In the text

	Fig. A.1 Example of telluric correction for the Hα and Na I doublet regions. After normalisation, the DRS processed data (in orange) was telluric corrected with Molecf it using the telluric model in blue, here shifted upwards for clarity.
In the text

	Fig. A.2 Same as Fig. 5 but for the lines detected in the visible range. The dashed black lines represent the wavelength position of the line.
In the text

	Fig. A.3 Posterior distribution of the investigated parameters in the DE-MCMC analysis for the four investigated lines. The excess of absorption c [%], offset [%], peak position, and FWHM correspond to the parameters we used in the Gaussian fit, while the jitter term σj, the semi-amplitude of the correlated noise h, the correlation length λ were used to parametrise the SE kernel within the GP.
In the text

Fig. A.4 Cross-correlation and Kp-Vsys maps of Na I, Ca I, and Fe I obtained with the Kitzmann et al. (2023) templates at 2000K. In the first column the CCF are depicted in terms of ppm. The dotted horizontal lines indicate the T1 and T4 contact points, while the red dashed line is the planetary trace. The Doppler shadow follows the white slanted line. Second columns are the Kp-Vsys maps. The planetary signal is expected to be limited between the green pointer in the middle of the figures, while the red lines point at the minimum value of the Kp-Vsys between Kp 75 km s−1, 125 km s−1 and RV -50 km s−1, +50 km s−1. The last column is the Kp-Vsys 1D map evaluated at the Kp minimum value. The grey region is the 3σ edge, while the vertical navy dashed line refers to the RV 0 km s−1.

In the text

	Fig. B.1 Same as Fig. A.3 but for SLOPpy.
In the text