The GAPS Programme at TNG LXXXI. The puzzle of the HeI absorption signal in the atmosphere of V1298 Tau b with GIARPS@TNG

Introduction

Understanding exoplanet evolution is crucial to explaining system architectures and atmospheric stability. A key feature is the “radius gap” (or “evaporation valley”) around $\sim$1.8 $R_{\oplus}$, separating dense, rocky super-Earths from larger, volatile-rich mini-Neptunes (e.g., Fulton et al. [33]; Fulton & Petigura [32]). Its origin likely involves both formation and evolutionary processes: migration may shape the sub-Neptune peak, while photoevaporation can produce the super-Earth population (Burn et al. [8]). Alternatively, mechanisms such as core-powered mass loss could produce the observed bimodal distribution on their own (e.g., Gupta & Schlichting [39]; Loyd et al. [55]). Studying atmospheric escape in planets near the radius gap therefore lets us test whether evaporation drives these demographics to understand the transition from gas-rich planets to stripped cores, especially in young systems where escape is most efficient (Sanz-Forcada et al. 2011 and references therein).

The primary driver of atmospheric erosion is stellar X-ray and extreme ultraviolet (XUV) irradiation (e.g., Lammer et al. [47]; Lecavelier des Etangs et al. [48]). As a star ages and spins down, its magnetic activity and XUV output decrease significantly (e.g., Owen [72]). Consequently, for most exoplanets orbiting mature stars, the atmospheric evaporation responsible for the radius gap has largely concluded. Planets around young stars ($\lesssim100$ Myr) therefore offer a unique opportunity to observe atmospheric escape in action during its most transformative epoch. In this context, the metastable helium triplet at 10832.06 Å, 10833.22 Å, and 10833.31 Å (vacuum wavelength He i, e.g., Seager & Sasselov [86]; Oklopčić & Hirata 2018; Spake et al. [94]) has emerged as a premier diagnostic of escaping atmospheres. The He i $\lambda$10833 Å triplet is an important astrophysical diagnostic, tracing neutral helium in its triplet state. It forms mainly through photoionisation by ultraviolet radiation ($\lambda< 504$ Å), followed by recombination into the metastable orthohelium $2\,{}^{3}\mathrm{S}$ state. Once in this state, He atoms can absorb incoming photospheric near-infrared photons, producing the He i 10833 Å absorption feature. Unlike traditional tracers such as Ly$\alpha$, the He i triplet is less sensitive to interstellar medium absorption and geocoronal contamination, and its accessibility from the ground enables high-resolution instruments to probe the dynamics of planetary outflows in detail (e.g., Nortmann et al. [67]; Allart et al. [2]).

The V1298 Tau planetary system is a key laboratory for studying atmospheric photoevaporation in the early stages of planetary evolution. The host star is a young ($\sim$10–30 Myr, David et al. [21]; Suárez Mascareño et al. [97]), rapidly rotating ($P_{\rm rot} \approx$2.91 days, Finociety et al. [30]) K-type dwarf. Its strong magnetic activity complicates radial velocity (RV) measurements (e.g., Suárez Mascareño et al. [97]; Di Maio et al. [22]; Damasso et al. [17]), but provides an ideal environment for studying extreme irradiation effects on young planetary atmospheres. Among its four known transiting planets (David et al. [21]), V1298 Tau b ($P \sim$ 24.14 days, $R \sim$ 9.9 $R_{\oplus}$; Table 2) is a prime target for transmission spectroscopy. Hubble Space Telescope (Barat et al. [5]) and James Webb Space Telescope (Barat et al. [4]) observations revealed an H/He-dominated atmosphere, detecting CO$_2$, H$_2$O, CO, CH$_4$, SO$_2$, and OCS, alongside a high “internal temperature” $T_{\rm int}$ of $\sim$500 K ($T_{\rm int}$ is a parameter of the temperature–pressure profile model from Parmentier & Guillot [75], as adopted by Barat et al. [4], that accounts for the intrinsic internal heat flux). These studies imply a low planetary mass of 12–15 $M_{\oplus}$, independently confirmed by Livingston et al. [51], whose dynamical estimate of $\sim$13.1 $M_{\oplus}$ was derived by modelling the transit timing variations (TTVs). These values are drastically lower than previous RV estimates (e.g., $M_{\rm P}$ = 203 $\pm$ 60 $M_{\oplus}$, Suárez Mascareño et al. [97]), and indicate that V1298 Tau b may be an inflated “gas-dwarf” progenitor (Barat et al. [4]) of mature sub-Neptunes. Despite predictions of ongoing atmospheric escape driven by the host star’s intense XUV irradiation (Poppenhaeger et al. [78]; Maggio et al. [56]), He i observations of V1298 Tau b remain inconclusive. Previous campaigns yielded either non-detections from partial transits (Vissapragada et al. [103]) or ambiguous absorption signals blended with stellar variability or the transit of planet d (Gaidos et al. [34]). This discrepancy between predictions and observations may stem from four main challenges: TTV-induced ephemeris uncertainties, the difficulty of capturing the long transit ($\sim$7 hours) from a single ground-based site, V1298 Tau d atmospheric contamination, and the strong, highly variable stellar He i background originating in the star’s active chromosphere (e.g., Guilluy et al. [35]). In young, active stars such as V1298 Tau, stellar signals can mimic or mask planetary absorption, creating a puzzle that is difficult to resolve without simultaneous monitoring of activity indicators.

Parameter	Prior	Value	Reference
Stellar parameters
Effective temperature, $T_{\rm eff}$ [K]	...	$5050^{+100}_{-100}$	SM21
Stellar radius, $R_{\star}$ [$R_{\odot}$]	...	$1.36 \pm0.097$	This work$^{(a)}$
Stellar mass, $M_{\star}$ [$M_{\odot}$]	...	$1.124 \pm0.015$	This work$^{(a)}$
Systemic velocity, $v_{\rm sys}$ [km s$^{-1}$]	...	$16.15 \pm0.38$	David et al. [21]
Projected rotational velocity, $v \sin i_{\star}$ [km s$^{-1}$]	$\mathcal{N}(24.0, 2.0)$	$24.5^{+1.3}_{-1.2}$	This work
Limb darkening coefficient, $\epsilon$	$\mathcal{N}(\mathrm{LDTk}, 0.2)$	$0.66^{+0.16}_{-0.17}$	This work
Stellar RV semi-amplitude, $K_{\star}$ [m s$^{-1}$]:
...	...	$41 \pm12$	SM21
...	...	$2.7 \pm1.1$	This work using $(b)$ and $M_{\rm p}/M_{\star}$ from L26
Planetary and Orbital parameters
Orbital period, $P$ [days]	$\mathcal{N}(\mathrm{Ephem}, 0.01)$	$24.1394^{+0.0047}_{-0.0047}$	This work
Transit epoch, $T_{0,1}$ [BJD$_{\rm TDB}$]	$\mathcal{N}(\mathrm{Ephem}, 0.1)$	$2458877.639^{+0.012}_{-0.012}$	This work (Night 1)
Transit epoch, $T_{0,2}$ [BJD$_{\rm TDB}$]	$\mathcal{N}(\mathrm{Ephem}, 0.1)$	$2459191.4276^{+0.0092}_{-0.0080}$	This work (Night 2)
Transit epoch, $T_{0,3}$ [BJD$_{\rm TDB}$]	...	$2460591.50759 \pm0.00059$	Leonardi et al. (in prep.) (Night 3)
Transit epoch, $T_{0,4}$ [BJD$_{\rm TDB}$]	...	$2460615.6503 \pm0.0050$	Leonardi et al. (in prep.) (Night 4)
Scaled semi-major axis, $a/R_{\star}$	$\mathcal{N}(27.0, 2.0)$	$27.35^{+0.96}_{-0.94}$	This work
Radius ratio, $R_{\rm p}/R_{\star}$	$\mathcal{N}(0.070, 0.005)$	$0.068^{+0.004}_{-0.004}$	This work
Inclination, $i$ [deg]	$\mathcal{N}(89.0, 1.0)$	$89.46^{+0.31}_{-0.32}$	This work
Spin-orbit angle, $\lambda$ [deg]	$\mathcal{U}(-180.0, 180.0)$	$31 \pm14$	This work
Planetary mass, $M_{\rm p}$ [$M_{\rm Jup}$]:
...	...	$0.63 \pm0.19$	Derived from $K_{\star}$ by SM21
...	...	$0.042 \pm0.017$	Derived from $M_{\rm p}/M_{\star}$ by L26
Planet RV semi-amplitude, $K_{\rm p}$ [km s$^{-1}$]:
...	...	$76.6 \pm1.2$	This work$^{(c)}$ using $K_{\star}$ from SM21
...	...	$76.6 \pm1.2$	This work$^{(c)}$ using $K_{\star}$ derived from L26
Eccentricity, $e$	Fixed	$0$	Consistent with a circular orbit (e.g., L26)
Transit duration, $T_{14}$ [h]	...	$6.95^{+0.32}_{-0.38}$	This work$^{(d)}$
Full transit duration, $T_{23}^{(4)}$ [h]	...	$6.00^{+0.31}_{-0.41}$	This work$^{(d)}$
Impact parameter, $b$	...	$0.26^{+0.15}_{-0.15}$	This work$^{(e)}$
Planetary radius, $R_{\rm p}$ [$R_{\rm Jup}$]	...	$0.902 \pm0.060$	This work
Noise model and GP parameters
RV jitter, $\sigma_{\rm RV}$ [m s$^{-1}$]	$\mathcal{U}(0.0, 50.0)$	$23.9^{+5.0}_{-5.1}$	This work
GP amplitude 1, $A_{\rm GP}^{1}$ [m s$^{-1}$]	$\mathcal{U}(0.0, 50.0)$	$23.0^{+2.1}_{-2.1}$	This work
GP timescale 1, $\tau_{\rm GP}^{1}$ [days]	$\mathcal{U}(0.0, 1.0)$	$0.016^{+0.011}_{-0.010}$	This work
GP amplitude 2, $A_{\rm GP}^{2}$ [m s$^{-1}$]	$\mathcal{U}(0.0, 50.0)$	$30.3^{+4.4}_{-4.3}$	This work
GP timescale 2, $\tau_{\rm GP}^{2}$ [days]	$\mathcal{U}(0.0, 1.0)$	$0.105^{+0.020}_{-0.016}$	This work

In this work, we used high-resolution GIANO-B observations of the metastable helium triplet to investigate the atmospheric outflow of V1298 Tau b while carefully accounting for its complex stellar environment. Our observational setup avoided simultaneous planet d transiting (Sect. 2), preventing spectral contamination, and guaranteeing an unambiguous analysis of V1298 Tau b’s atmospheric signal. The paper is organised as follows: Sect. 2 describes the observations, while Sect. 3 includes an analysis of the Rossiter-McLaughlin effect, the data reduction, and our findings, which are presented and contextualised within the framework of stellar activity. In Sect. 4, we discuss the role of star–planet interactions in shaping the He i signature, contextualize our results within previous literature, and model planetary evolution. We summarise our conclusions in Sect. 5.

Observations

We observed four transits of V1298 Tau b, namely UT 28 January 2020 (Night 1), UT 07 December 2020 (Night 2), UT 07 October 2024 (Night 3), and UT 31 October 2024 (Night 4). A log of the observations is provided in Table 1. Observations were acquired exploiting the GIARPS observing mode (Claudi et al. [16]) at the Telescopio Nazionale Galileo (TNG). GIARPS enables simultaneous high-resolution observations in both the visible (0.39–0.69 $\mu$m) with HARPS-N ($R \approx$ 115,000) and the near-infrared (0.95–2.45 $\mu$m) with GIANO-B ($R \approx$ 50,000). During Night 3 and 4, some instrumental problems occurred that prevented us from using the GIARPS mode for the entire observing night. Since He i was the main focus of our research, it was decided to prioritise observations with GIANO-B. Therefore, we only have a few HARPS-N pre-transit and ingress spectra for Night 3, and a few in-transit spectra for Night 4.

	Date [UT]	Instrument	$N_{\mathrm{obs}}$	Exp. time [s]	Airmass min→mean→max	S/N min→mean→max	PI Program	PI Name
Night 1	2020-01-28	GIANO-B	56[55]	300	1.01→1.21→2.05	26.1→40.1→52.0	GAPS2	G. Micela
		HARPS-N	47[47]	400		23.1→39.1→50.9
Night 2	2020-12-07	GIANO-B	80[72]	300	1.01→1.21→2.07	26.0→44.2→53.82	GAPS2	G. Micela
		HARPS-N	71[71]	400		11.1→28.7→56.5
Night 3	2024-10-07	GIANO-B	58[53]	300	1.01→1.12→1.57	29.8→41.4→48.8	CAT24B_81	J. Orell Miquel
Night 4	2024-10-31	GIANO-B	76[72]	300	1.01→1.25→1.83	30.7→39.9→48.0	CAT24B_81	J. Orell Miquel

For the GIANO-B observations, we applied an ABAB nodding pattern (Claudi et al. [16]) to improve thermal background noise subtraction. GIANO-B is an echelle spectrograph covering four near-infrared bands ($Y$, $J$, $H$, $K$), divided into 50 orders. The HARPS-N observations were performed with the B-fibre positioned on the sky. This instrument, also an echelle spectrograph, provides spectral coverage across 69 distinct orders for each fibre.

Based on the TTV dynamical model from Leonardi et al. (in prep.), we calculated the orbital phases of all known planets in the system for each observing run. None of the other planets were transiting during our observations, ensuring that our GIANO-B and HARPS-N spectra are free from atmospheric contamination by secondary planets (Fig. A.1). Night 1 and Night 4 comprise observations taken before and during the early phases of the transit, while Night 2 and Night 3 cover spectra acquired during and after the transit. During Night 4, the target was not fully in focus for some observations. GIANO-B experienced guiding problems around 01:30 UT (orbital phase $\sim$−0.0022), but the tracking stabilised by 01:52 UT (orbital phase $\sim$0.0029), allowing observations to continue.

Furthermore, we observed that on Night 3 and 4, towards the end of the night, specifically from UT 04:45:06 on Night 3 and UT 05:12:46 on Night 4 until the end of the night, an emission line appears in the GIANO-B raw spectra at the position of He i in the Earth’s rest frame (Fig. A.2 in the Appendix). Unfortunately, on these two nights the difference between the systemic velocity ($v_{\rm sys}$) and the barycentric Earth radial velocity ($v_{\rm berv}$) is small (on average $\sim$4.2 km s$^{-1}$ and $\sim$−5.8 km s$^{-1}$ for Night 3 and 4, respectively), so this emission feature falls directly over the stellar helium line. The line corresponds to the metastable He($2^{3}$S) airglow. First detected in 1959 (Mironov et al. [61]), it probes the Earth’s upper thermosphere via solar-induced fluorescence of He($2^{3}$S), peaking at altitudes above 250 km with densities of a few atoms cm$^{-3}$ (Waldrop et al. [105]). Crucially, because this emission arises just before nautical twilight and grows progressively stronger over time, the standard A-B (B-A) nodding subtraction cannot remove it effectively. For our analysis, we therefore masked the pixels affected by this feature in the individual images. We emphasise that for future observations, especially those taken near twilight, the difference between $v_{\rm sys}$ and $v_{\rm berv}$ should be kept large enough to prevent this airglow from contaminating the stellar helium line.

Notes. Columns denote the observing night, instrument, total number of acquired spectra (with the analysed subset in brackets), exposure time, airmass, signal-to-noise ratio (S/N), and Principal Investigator (PI) program and name. S/N values are computed at order #39 (He i) for GIANO-B, and order #64 (H$\alpha$) for HARPS-N. For GIANO-B, spectra with S/N below 25 were discarded. Airmass and S/N for HARPS-N observations during Nights 3 and 4 are omitted, as the low number of spectra precluded the H$\alpha$ analysis (Sect. 2).

Analysis and results

Rossiter–McLaughlin effect analysis

We used the HARPS-N spectra gathered in Night 1 and Night 2 to obtain an independent constraint on the stellar obliquity of V1298 Tau relative to the normal of planet b’s orbital plane, and to measure two new, accurate values of the mid-transit time ($T_0$) at these two epochs. The latter is essential given that V1298 Tau b is known to exhibit large TTVs due to mutual perturbations from the other planets in the system. Previous work has already detected the RV anomaly produced by the Rossiter-McLaughlin (RML) effect for V1298 Tau b. Gaidos et al. [34] obtained in-transit spectra during a single night with the IRD spectrograph on the 8.2-m Subaru telescope, measuring a projected spin-orbit angle of $\lambda= 15^{+15}_{-16}$ deg and $\lambda= 2^{+12}_{-4}$ deg from the apparent RV shift and from changes in the line profile, respectively. Johnson et al. [41] observed a partial transit with Keck/HIRES (missing egress) and LBT/PEPSI (missing the second half and egress), obtaining $\lambda= 4^{+7}_{-10}$ deg from the extracted RV sequences for both instruments. The consensus is that V1298 Tau b lies on an aligned orbit.

For our independent assessment of the spin-orbit alignment of V1298 Tau b, we extracted RVs from the HARPS-N spectra using the template-matching approach implemented in the TERRA pipeline (Anglada-Escudé & Butler [3]). For a young, active star such as V1298 Tau, this approach is preferable to the method based on a Gaussian fit to the cross-correlation function (CCF) of the observed spectra with a binary mask, as implemented in the standard HARPS-N Data Reduction Software pipeline (Dumusque et al. [24]). To model the observed RML signal, we adopted the formalism of Ohta et al. [68], which is well suited for fitting the RML RV signal from template matching. This model was implemented within a differential evolution Markov chain Monte Carlo (DE-MCMC) algorithm (Ter Braak [99]; Eastman et al. [25]). On both Night 1 and Night 2, the RV sequences show significant variability attributable to stellar activity. We included a Gaussian process (GP) framework in the RML modelling to capture this correlated RV noise, adopting a Matérn-3/2 covariance kernel (Stangret et al. [96]).

In the Ohta et al. [68] formalism, the RML model depends on the following parameters: $\lambda$, $T_0$, $P$, the projected stellar rotational velocity $v \sin i_\star$, a linear limb-darkening coefficient $\epsilon$, the inclination angle of the planetary orbit $i$, the radius ratio $R_{\mathrm{p}}/R_\star$, and the orbital semi-major axis scaled in units of the stellar radius $a/R_\star$. The common approach is to adopt $\lambda$, $T_0$, $v \sin i_\star$, and $\epsilon$ as free parameters, while fixing the others to high-precision determinations from the literature (e.g., Palle et al. [74]; Stangret et al. [96]). In the case of V1298 Tau b, this approach risks introducing biases and inconsistencies, primarily because of (a) the heterogeneity of the methods used to model spectroscopic and photometric datasets of variable quality, (b) the strong dynamical perturbations in the system, which make the use of instantaneous (osculating) orbital elements (e.g., $P$) a potentially poor choice, and (c) the intrinsic difficulty of obtaining consistent estimates of the stellar parameters for the young primary.

Table A.1 compares published stellar and planetary parameters drawn from a range of literature sources, including David et al. [21] (D19a), David et al. [21] (D19b), Gaidos et al. [34] (G21), Vissapragada et al. [103] (V21), Feinstein et al. [29] (F21), Feinstein et al. [28] (F22), Johnson et al. [41] (J22), Suárez Mascareño et al. [97] (SM21), Finociety et al. [30] (F23), Sikora et al. [88] (S23), and Livingston et al. [51] (L26). As the table shows, the complete set of parameters is not always reported across studies, and discrepancies between reported values can reach or exceed the $3\sigma$ level, the most significant differences being, as expected, those in the $P$ values. We thus chose to treat the full set of RML model parameters as free, adopting wide uninformative priors or broad Gaussian priors that comfortably bracket the ranges of literature values, in order to optimise the exploration of the parameter space (Table 2). The global RML+GP model contains a total of 14 free parameters, including separate values of the epoch of transit centre ($T_{0,1}$, $T_{0,2}$), the GP amplitude ($A_{\mathrm{GP}}^{1}$, $A_{\mathrm{GP}}^{2}$), and the GP timescale ($\tau_{\mathrm{GP}}^{1}$, $\tau_{\mathrm{GP}}^{2}$) for Night 1 and Night 2, respectively. The priors on $T_{0,1}$ and $T_{0,2}$ were centred on the predicted ephemerides (Ephem) derived from the dynamical TTV model of L26, while the prior on the linear limb-darkening coefficient was centred on the value obtained with LDTk (Parviainen & Aigrain [76]) for a star with the same effective temperature as V1298 Tau and a spectral range corresponding to that of HARPS-N.

The global DEMCMC fit converged to a set of posterior distributions whose medians and $1\sigma$ uncertainties (obtained from the $\pm34.13\%$ intervals of the posteriors) are reported in Table 2. Figure 2 shows the phase-folded HARPS-N RV time series, with the stellar activity contribution removed in the top panel and the RML effect removed in the bottom panel. The RML-only and GP-only model components are also overplotted. Overall, the central value of $\lambda$ ($31 \pm14$ deg) differs from zero at the $2.2\sigma$ level, indicating a possible departure from exact spin-orbit alignment. Nevertheless, the large uncertainty on the spin-orbit angle makes the retrieved value consistent with previous literature estimates.

Night 3 and Night 4 $T_0$ derivation

Because the HARPS-N spectra obtained during Night 3 and Night 4 were insufficient to carry out an RML analysis, we relied on results from the ongoing TTV follow-up of the V1298 Tau system (Leonardi et al. in prep.) to constrain the time of inferior conjunction. For Night 3, V1298 Tau b was observed simultaneously with CHEOPS ([6]) as part of the GO programme PR250033 (PI: Leonardi), allowing a precise determination of the transit mid-time ($T_{0,3} = 2460591.50759 \pm0.00059$ BJDTDB). However, no simultaneous photometric observations were available for Night 4, preventing a direct measurement of the mid-transit time. In this case, we adopted the predicted time of inferior conjunction from the dynamical TTV model of Leonardi et al. (in prep.) ($T_{0,4} = 2460615.6503 \pm0.0050$ BJDTDB). As a consistency check, we compared the mid-transit times derived from the RML analysis for Nights 1 and 2 with the predictions of the TTV dynamical model. The two independent determinations agree within $1.3\sigma$, providing additional confidence in the adopted ephemeris.

GIANO-B He I transmission spectroscopy

To probe the He I absorption in the extended atmosphere of V1298 Tau b, we focused on GIANO-B echelle order #39. Telluric features, particularly water vapour lines near 10835.1 Å, were corrected using the ESO Molecfit software ([89, 43]). This spectral region also contains two OH emission doublets close to the He I triplet: the first at ∼10832.1 Å and ∼10832.4 Å, and the second at ∼10834.2 and ∼10834.3 Å (vacuum wavelengths). Although the A-B nodding acquisition automatically removes these emission features, variations in atmospheric seeing during the night can leave residual signals. We therefore visually identified the positions of these OH lines in the raw spectra and masked them to exclude any potential contamination from our final transmission spectra (e.g., [36, 19] (Guilluy et al. 2024)).

The next step was to isolate the planetary signal via transmission spectroscopy. We shifted all spectra to the stellar rest frame and normalised them to unity. We then computed a master out-of-transit spectrum, $S_{\mathrm{out}}(\lambda)$, by averaging the out-of-transit spectra with S/N $> 30$. To avoid potential early absorption signals before first contact ($t_1$), particularly evident in Night 1 (Figure 3), we restricted the pre-transit baselines to orbital phases $< -0.009$ for Night 1 and $< -0.0065$ for Night 4.

Individual transmission spectra, $T(\lambda,i)$, were obtained by dividing each spectrum by $S_{\mathrm{out}}(\lambda)$ and shifting it to the planetary rest frame (using the parameters from Table 2). The resulting 2D transmission maps are shown in Figure 3. For each night, we also computed a weighted-average in-transit transmission spectrum between contact points $t_2$ and $t_3$ (Figure 4). We modelled these average profiles with a DEMCMC algorithm (Ter Braak [99]), fitting a Gaussian profile (peak position, full width at half maximum, absorption depth $c$, and continuum offset). Correlated noise was mitigated through GP regression with a squared-exponential kernel (e.g., [19] (Guilluy et al. 2024)), while a jitter term $\sigma_j$ accounted for uncorrelated noise. Using the publicly available PARVATI tool (Rainer et al. in prep.), we extracted the equivalent widths (EWs) of the GP-corrected weighted-average in-transit transmission spectra within a $\pm3\sigma_{\mathrm{gauss}}$ range around the line peak, resulting in 211.1 $\pm$ 4.9, 106.3 $\pm$ 2.1, 19.2 $\pm$ 1.7, and 123.9 $\pm$ 2.1 mÅ for

Notes. $^{(a)}$ Weighted average of the values reported in the literature. $^{(b)}$

$^{(c)}$

$^{(d)}$

$^{(e)}$

Night 1, Night 2, Night 3, and Night 4, respectively. The main parameters of the Gaussian fit, as well as the EWs for each night, are reported in Table 3.

Night	$\lambda_0$ [Å]	$v$ [km/s]	$c$ [%]	$\sigma_{\mathrm{Gauss}}$ [nm]	EW(He i) [mÅ]
1	$10831.84^{+0.12}_{-0.11}$	$-39.5^{+3.2}_{-3.1}$	$8.47^{+0.72}_{-0.71}$	$0.1127^{+0.013}_{-0.011}$	$211.1 \pm4.9$
2	$10831.72^{+0.12}_{-0.12}$	$-42.7^{+3.3}_{-3.4}$	$3.98^{+0.40}_{-0.39}$	$0.1114^{+0.016}_{-0.013}$	$106.3 \pm2.1$
3	$10832.78^{+0.17}_{-0.33}$	$-13.2^{+4.8}_{-9.2}$	$1.58^{+0.73}_{-0.81}$	$0.0453^{+0.048}_{-0.018}$	$19.2 \pm1.7$
4	$10832.43^{+0.12}_{-0.12}$	$-23.0^{+3.2}_{-3.2}$	$4.49^{+0.43}_{-0.42}$	$0.1122^{+0.015}_{-0.014}$	$123.9 \pm2.1$

Finally, we calculated the EWs of the individual transmission spectra using a uniform integration passband across all four nights (Figure 5). Specifically, this window was defined by computing the weighted mean of the Gaussian centroid ($l_0$) and $\sigma_{\rm Gauss}$ derived from the DEMCMC fits to the nightly mean transmission spectra, setting the integration bounds to $l_0 \pm1\sigma_{\rm Gauss}$, that is, 10830.96–10833.14 Å. We selected a restricted $1\sigma_{\rm Gauss}$ range (rather than a broader window) to limit contamination from the telluric line at $\sim1083.5$ nm. Since uneven phase coverage across the observing nights can make full-transit averages misleading, we divided the $t_2$–$t_3$ transit chord into three phase intervals: $R_1$ ($-t_2 < \phi< -0.0015$), $R_2$ ($-0.0015 < \phi< 0.0015$), and $R_3$ ($0.0015 < \phi< t_3$). The temporally resolved EWs for each night and phase bin are reported in the bottom panel of Figure 5. The results reported here were obtained using the SM21 mass value. However, they do not change if the mass determination of L26 is adopted, since the derived $K_{\rm p}$ for a 24-day orbital period depends only weakly on the planetary mass. Given that our measured $T_{14}$ values for both Night 1 and Night 2 (Table 2) are higher than literature estimates (David et al. (2019a,b; Feinstein et al. 2021; Suárez Mascareño et al. 2021; Livingston et al. 2026)), though still compatible within $1.4\sigma$, we tested the robustness of our analysis by shifting the transit center to $T_0 \pm1\sigma$. The resulting contrast ($c\%$) and EWs of the final averaged transmission spectra remain compatible with our findings.

The puzzle of the He I signal

Figure 5 illustrates the time series of the excess-absorption EWs derived from the GIANO-B transmission spectra around the He I

triplet. The EWs exhibit an asymmetric shape, characterised by stronger absorption before orbital phase $\sim0$ than after. This asymmetry, corroborated by the 2D maps in Figure 3, suggests the presence of a helium structure preceding the planet along its orbit. Furthermore, the 2D maps and Figure 4 reveal a noticeable blueshift. However, this shift is not uniform across all observations: Night 1 appears to be in good agreement with Night 2, while Night 3 and Night 4 exhibit a less pronounced displacement. Notably, Night 3 also displays a smaller FWHM and depth $c$ than the other nights (Table 3). Several caveats must be considered in this analysis. First, the averaged transmission spectra were obtained over different time intervals owing to the uneven phase coverage of the observing nights. Second, the data are affected by telluric contamination from OH emission lines and H$_2$O absorption. Although the OH lines were masked and the water vapour line was corrected using molecfit, some residuals may persist in the final spectra. A further complication is that the three He i triplet lines are blended with one another and with the Si i photospheric line (Fig. A.3).

In-transit EW variability among the different nights is also visible. The maximum variability is observed in Range 1 between Night 4 and Night 1 at $4.1\sigma$. While in Range 2 the largest variability is observed between Night 2 and Night 4, and in Range 3, between Night 3 and Night 4 but the nights remain consistent within $3\sigma$. The variability in Range 1 between Night 1 and Night 4 is also visible during ingress.

Given the high activity level of V1298 Tau, we investigated whether this could drive the observed variability. To this end, we examined the stellar He i line (Sect. 3.5) alongside the log $R'_{\mathrm{HK}}$, H$\alpha$, and Na i activity indicators (e.g., Guilluy et al. [35]) extracted from the HARPS-N spectra (Sect. 3.6).

Investigation of the stellar helium line

The master-out spectra exhibit significant night-to-night variability at the He i triplet (Appendix Fig. A.3, top panel). Using PARVATI, we calculated the stellar EWs(He i) by integrating the flux over the 10831.5–10835.0 Å range. This range differs from the 10830.96–10833.14 Å interval used for the individual planetary transmission EWs (Sect. 3.3). The stellar integration range was chosen to capture the He i core while avoiding nearby H$_2$O and blended Si i lines, whereas the planetary range was shifted to capture the absorption feature blueshifted relative to the helium triplet in the 2D maps (Figure 3). The resulting EWs are 265.0 ± 1.5, 219.2 ± 1.2, 172.9 ± 1.4, and 269.2 ± 2.7 mÅ for Night 1, Night 2, Night 3, and Night 4, respectively, yielding a maximum variation of ∼96 mÅ. This high variability is also evident in the EWs of the individual spectra (Fig. A.5, upper panel). To provide context, we monitored the nearby photospheric Si i line at ∼10830.1 Å (Fig. A.5, bottom panel). While the Si i line remains remarkably stable across all four nights, the He i triplet fluctuates significantly. The difference between the He i and Si i EWs (Fig. A.5, lower panel) reveals an in-transit excess, raising the question of whether this signal stems from an evaporating planetary atmosphere.

However, any in-transit signal must be interpreted with caution, given the well-documented intrinsic variability of V1298 Tau ([103, 34, 46, 1]). Previous studies have reported severe stellar He i fluctuations: ∼0.2 Å over a two-month baseline ([103]), variations of up to 50% during a single transit (likely due to the extended gas tail of planet d, [34]), and short-term changes of ∼80 mÅ over just three weeks ([46]). Due to this variability, Krolikowski et al. [46] suggested a highly conservative 95{sup}th percentile detection limit of ∼6.7%.

Ultimately, stellar activity remains a critical bottleneck that can easily mimic or mask planetary He i absorption (e.g., [82, 11, 35]). Our findings reinforce that, for active young stars like V1298 Tau, acquiring contiguous out-of-transit spectra immediately before and after the transit is essential to robustly isolate planetary signals ([102, 46]). This stands in stark contrast to older, inactive systems, where multi-epoch baselines are usually sufficient (e.g., [71] (Gully-Santiago et al. 2024)).

Stellar activity monitoring in the visible

log $R'_{\mathrm{HK}}$

We first examined the log $R'_{\mathrm{HK}}$ chromospheric activity index using HARPS-N data. The index was extracted from the spectra via the YABI platform (Hunter et al. [40]), hosted at the IA2 Data Center³, assuming a $B-V$ colour index of 0.99 mag. The resulting time series are shown in Figure 6. As detailed in Sect. 2, instrumental issues during Night 3 and Night 4 limited the datasets to only a few pre-transit/ingress spectra and a small number of in-transit spectra, respectively. Conversely, Night 2 features a detectable rise in the chromospheric index, which we attribute to surface inhomogeneities. This interpretation is corroborated by simultaneous photometric observations that we performed with the Nordic Optical Telescope (NOT), which revealed a strong anomaly on the light curve that is attributable to a large starspot (complex) on the surface of the parent star (Appendix B). This perfectly aligns with the cross-correlation function (CCF) analysis by Di Maio et al. [22], who mapped a persistent circumpolar spot visible throughout the rotation, alongside rotationally modulated, low-latitude active regions. Such a scenario is further supported by TESS light-curve analyses (private communication). Recently, Murphy et al. [64] found that V1298 Tau has starspots with distinct umbral and penumbral structures with thermal and geometrical properties similar to those of the Sun.

H-alpha

To assess whether the variability observed in the He i absorption is driven by stellar activity, we analyzed the Hα chromospheric line at 6562.8 Å collected simultaneously to our near-infrared observations as already done in Guilluy et al. (2024) and in D’Arpa et al. [19]. While the He i triplet is primarily populated in the upper chromosphere and transition region via XUV irradiation, the Hα line originates in the lower-to-middle chromosphere through the $n = 3 \to2$ transition. As a result, Hα serves as a robust tracer of chromospheric features, such as active regions (plages), filaments, and flares. Unlike other young stars of similar activity, which typically display Hα in emission, Feinstein et al. [29] noted that V1298 Tau shows a purely absorption profile without a significant core component. Following the strategy adopted for He i, we analysed the master out-of-transit spectra, which revealed a large variability (bottom panel of Fig. A.3 in the appendix), with EWs(Hα) values of 279.2 ± 1.4 mÅ, 270.6 ± 1.4 mÅ, and 329.4 ± 2.4 mÅ, for Night 1, Night 2, and Night 4, respectively. Left and middle panels of Fig. A.4 show how the EWs of the He i (left panel) and Hα (right panel) stellar lines (derived from the master-out spectra) vary as a function of $\log R'_{\mathrm{HK}}$. While the Hα line exhibits solar-like behaviour characterised by weaker absorption in more active regimes due to chromospheric filling-in ([13, 12]), the He i triplet show a departure from this pattern. In the Sun, He i absorption typically intensifies in active regions, with the EW increasing alongside stellar activity. In the solar case, Hα and He i become typically correlated during relatively strong flares, often preceded by a marked dimming (Kerr et al. [44]). In the solar chromosphere, the He i triplet forms via absorption from the $2^{3}\mathrm{S}$ atomic level, which is populated through the recombination of photoionised helium atoms (the photoionization-recombination PR mechanism). Ionising helium requires photons with wavelengths shorter than 503 Å, originating from strong emission in the transition region and corona of active regions. As the activity level increases, this high-energy emission grows, increasing the population of metastable helium and thus resulting in a deeper absorption line and a larger EW. For active dwarf stars, where the chromosphere is hot and dense enough, there is an alternative mechanism that contribute to populate the He i state directly from the ground state via collision (the collisional excitation, CE, mechanism, Sanz-Forcada & Dupree [83]). In dwarf stars, the increasing levels of activity do not cause more substantial He i absorption. On the contrary, the absorption appears rather shallow (if not partly in emission) for very active stars. In V1298 Tau, we observe a decrease in EW(He i) with increasing activity, suggesting that CE mechanism is also contributing to the population of the He i $2^{3}\mathrm{S}$ level. However, since the averaged $\log R'_{\mathrm{HK}}$ values in Fig. A.4 span different temporal windows (Figure 6), this trend warrants caution.

Following a similar methodology applied to the He i triplet, we then derived Hα transmission spectra. We extracted them following the framework described in detail in D’Arpa et al. [20]. The resulting transmission spectra 2D maps are displayed in the left panel of Figure 7, while the right panel illustrates the time-series EWs of these transmission spectra. Notably, Night 2 exhibits a reduced EW during the transit, accompanied by a distinct emission-like feature in the 2D tomography maps. This signal is most likely attributable to stellar activity. This pseudo-emission feature, which manifests as a relative decrease in the stellar Hα absorption depth, could stem from either a stellar flare occurring during the transit (inducing chromospheric filling-in) or the occultation of a chromospheric active region. Simultaneous NOT photometry strongly supports the latter scenario (Appendix B), indicating that the planet crossed a large starspot (complex) during Night 2. In a standard framework, crossing a chromospheric active region (where the Hα line is typically shallower due to localized emission filling-in) is expected to remove a lower-absorption component from the integrated disk, thereby resulting in a net absorption anomaly in the transmission spectrum (an increase in the observed absorption depth). However, our observations reveal the opposite behavior (i.e., a net pseudo-emission signature). To reconcile this, we hypothesize a geometric configuration where the chromospheric active region (plage) is significantly more extended than both the planetary disk and the underlying photospheric spot. This is consistent with solar observations, where chromospheric faculae generally encompass much larger areas than their photospheric counterparts. In this framework, the specific chord crossed by the planet may be locally less active (characterized by a lower emission contrast) than the surrounding, unocculted sections of the plage. Consequently, as the planet blocks this relatively quieter zone, the emission from the dominant, unocculted active region decreases the net integrated Hα absorption depth of the star during the transit. Regarding the He i signal during Night 2, we observe the opposite trend: an increased absorption (a deeper line profile) at mid-transit (Figure 3 and Figure 5). As discussed in Sect. 3.6, the stellar He i absorption depth in V1298 Tau decreases with increasing stellar activity due to chromospheric filling-in. Therefore, a spot-crossing event involving an active region with a shallower He i profile results in a net absorption peak in the in- and out-of-transit ratio, as the “less deep” component is temporarily blocked by the planet. The simultaneous decrease in Hα absorption depth and the increase in He i absorption can be physically reconciled by assuming that the formation regions of these two indicators are spatially decoupled (e.g., forming at different chromospheric heights). This atmospheric stratification allows the planet to simultaneously occult a region that appears less active in Hα relative to its immediate surroundings, while remaining sufficiently active in the He i-forming layers to generate the observed absorption enhancement. We note that a remarkably similar behavior, an increase in He i absorption coupled with an Hα emission anomaly, has been reported by Orell-Miquel et al. [70] for the 650-Myr-old sub-Neptune HD 235088 b. Yet, these considerations are qualitative; further observations providing full transit coverage, including robust pre- and post-transit baselines, are essential to disentangle the planetary signal from the complex stellar activity of this system. In contrast, the Hα transmission spectroscopy for Night 1 reveals a marginal absorp tion signal. Due to the low $K_{\mathrm{p}}$ value (Table 2), which is a consequence of the planet’s relatively long orbital period, it is challenging to discern whether this feature is better aligned with the planetary or the stellar rest frame. However, the fact that Night 2 displays a pronounced emission signal suggests that the features observed in both nights are dominated by stellar activity, and reflect the complex stellar surface.

Nai

In addition to Hα and $\log R^{\prime}_{\mathrm{HK}}$, we also analysed the Na i doublet. Our analysis reveals a significant night-to-night variability. We observe a clear emission signal in Night 2. Conversely, Night 1 is affected by significantly larger uncertainties that preclude a definitive conclusion. A comprehensive discussion about the analysis is reported in Appendix C.

Discussion

Different stellar activity levels can lead to different He i planetary absorption signals

An additional complication in separating the stellar contamination from the planetary signal arises from the fact that the observed planetary absorption may be altered by variations in the stellar flux caused by activity. In order to test how different activity levels affect the helium absorption profile, we simulated the evaporation of V1298 Tau b with the hydrodynamic code ATES (Caldiroli et al. 2021), using different XUV irradiation levels. We scaled the reference XUV spectral energy distribution (SED, [57]) by factors of 0.5, 1, and 2 relative to its baseline integrated flux $F_{0} = 3.73 \times10^{4}\,\mathrm{erg\,s^{-1}\,cm^{-2}}$. The simulations were performed for both the tabulated mass values (Table 2) fixing the H/He ratio at 100/10; the results are presented in Table A.2. Although core-powered mass loss is not accounted for in ATES and could complicate the atmospheric escape picture, we expect it to be subdominant over the range of planetary masses and equilibrium temperatures considered for V1298 Tau b ([73, 62]).

In the case of the higher planetary mass value from SM21, the ATES simulations indicate that increasing the XUV flux leads to a decrease in the predicted transit depth. Since the outflow is highly ionized, a higher number of ionizing photons lowers the density of singly ionized helium – which would otherwise supply the metastable state via recombination – in the diffuse, optically thin region of the outflow, which dominates the absorption signal (Fig. A.6, bottom panel). Accordingly, the predicted depth and equivalent width of the metastable helium line decrease as the XUV flux increases. The opposite trend is found for the lower-mass L26 case: increasing the XUV flux produces a stronger outflow that remains dominated by neutral and singly ionized species. As a result, the population of metastable helium grows with irradiation (Fig. A.6, top panel), leading to a larger predicted line depth and equivalent width.

It is important to note that reproducing the exact absolute values of the observed He i depth and EW remains a general challenge for 1D hydrodynamic frameworks. As recently discussed by Taylor et al. [98], 1D models frequently underestimate the total observed EW because they omit additional non-thermal broadening mechanisms. Processes such as 3D atmospheric dynamics, plasma turbulence, or magnetic interactions with an active host star can significantly broaden the absorption line profile. This extra velocity dispersion yields observed EWs that are systematically larger than those predicted by purely radial outflow models. The observed depths for Night 2, Night 3, and Night 4, plotted as a function of $\log R^{\prime}_{\mathrm{HK}}$ (last panel of Fig. A.4), appear to follow the trend predicted by the simulations when assuming SM21 mass: higher stellar activity leads to a greater XUV flux and, consequently, a reduced He i absorption depth. The only clear exception is Night 1. The discrepancy in Night 1 can be attributed to two factors: an enhanced star-planet interaction, suggested by the leading tail observed exclusively on this night (Sect. 4.3), and potential residual telluric contamination, given that the 10835.1 Å telluric line completely overlaps with the stellar He i line during this particular run (Sect. 3.4; Figure 1). It is worth noting, however, that because the measured absorption depths are derived from averages over different temporal windows (e.g., Night 3 covers only the final part of the transit), this trend should be interpreted with caution.

Comparison with previous works

Previous analyses of the He i triplet in V1298 Tau b resulted in non-detections from partial transits ([103]) or ambiguous results ([34]). Specifically, Gaidos et al. [34] analyzed multiple observations of V1298 Tau b across several nights in 2020, covering the transit along with pre- and post-transit baselines. To minimize baseline systematics and stellar variability, we construct transmission spectra exclusively using out-of-transit references obtained on the same night as the transit, whereas Gaidos et al. [34] used out-of-transit spectra obtained on the preceding night. In their study, they found an elevated absorption during the transit relative to the pre-transit baseline; however, the absorption decreased throughout the transit, with an equivalent width variation of $\Delta\mathrm{EW}(\mathrm{He\,I}) \approx0.19$ nm. Gaidos et al. [34] noted that this behavior could be explained by a He-bearing cloud preceding V1298 Tau b in its orbit. Their observation sequence began only 1 h after the transit of planet d, introducing potential contamination from this interior planet. We avoided this issue by selecting observation windows free from orbital phase overlap (Sect. 2). While a direct quantitative comparison between our EW(He I) values and those of Gaidos et al. [34] is not feasible due to differences in calculation methodology, we can compare the shape of the absorption. Notably, the feature shown in their Fig. 8 appears fully compatible with both the wavelength amplitude and the blueshift relative to the He I rest wavelengths observed in our dataset (Figure 3).

The role of star-planet interactions in shaping the He I transit signature

In the previous sections, we highlighted the difficulty in disentangling activity-induced pseudo-signals from real planetary absorption for a highly active star such as V1298 Tau. We also demonstrated how the expected metastable He I absorption is strongly modulated by the incident stellar flux, which governs the population of the ${}^{2}{}^{3}\mathrm{S}$ state. These combined factors could potentially contribute to the diverse absorption levels observed across the different investigated epochs. Furthermore, Night 1 exhibits a distinct blueshifted leading tail, a feature that is counterintuitive based on simple orbital kinematics and the conservation of angular momentum, which would typically predict a red-shifted signal in a trailing wake. Crucially, this blueshifted signature is only visible during Night 1 and not during Night 4. According to the log $R'_{\mathrm{HK}}$, Night 1 shows a higher activity level than Night 4. This may suggest that the anomalous velocity component responsible for the blueshift is not an intrinsic feature of the planetary outflow under quiescent conditions, but is instead directly modulated by variations in the stellar environment (XUV flux, wind properties, and magnetic configuration) and the local plasma regime governing the star-planet interaction. Additional insight is provided by the averaged in-transit transmission spectra (Figure 4), where Night 1, displaying the clearest pre-ingress absorption, is also associated with the strongest and broadest He I. This combination disfavours a scenario in which the leading feature is produced by a purely static, magnetically trapped reservoir alone. Instead, it suggests that the absorbing material has a significant dynamical component (contributing to both the line broadening and the net Doppler shift, Table 3), while still being shaped and, at times, concentrated by confinement within the planet’s environment.

The nature of the star-planet interaction depends on the local Alfvénic Mach number $M_{\mathrm{A}} = u/u_{\mathrm{A}}$, defined as the ratio of the stellar wind velocity to the local Alfvén speed, which is defined as $v_{\mathrm{A}} = B/\sqrt{4\pi\rho}$, where $B$ and $\rho$ are the magnetic field strength and plasma density of the wind at the planetary orbit. When this ratio exceeds unity (super-Alfvénic regime), a bow shock can form; when it is below unity (sub-Alfvénic regime), no shock develops. The Alfvénic Mach number at the planetary orbit depends on stellar wind parameters, such as the stellar mass-loss rate and surface magnetic field strength, that remain poorly constrained, making it difficult to determine unambiguously which regime applies. Moreover, variations in stellar activity level may cause the system to alternate between the two regimes because such variations are accompanied by changes in the intensity and in the configuration of the stellar field close to the planet.

We conjecture that V 1298 Tau b is located in a region where the Alfvénic Mach number of the stellar wind is of the order of unity ($M_{\mathrm{A}} \sim1$) for most of the time. In such a regime, assuming a planetary dipole field of $B_{\mathrm{p}} \sim25$ G and a stellar dipole field of $B_{\mathrm{s}} \sim200$ G (Finociety et al. [30]), the extension of the planetary magnetosphere in the direction of the impinging stellar wind is of $\sim12 R_{\mathrm{p}}$, where $R_{\mathrm{p}}$ is the radius of the planet (Appendix § 10). For a planet mass of $\sim0.042\,M_{\mathrm{Jup}}$ (L26), the low surface gravity of the planet and the heating of the planetary atmosphere by the strong stellar XUV flux make its evaporation remarkable, thus filling all the magnetosphere with plasma. The large pressure scale height of the plasma confined inside the magnetosphere, $H_{\mathrm{p}} \sim1.6 R_{\mathrm{p}}$, makes the density difference between the surface of the planet and the boundary of the magnetosphere less than a factor of two for a magnetospheric plasma temperature of $T = 10^{4}$ K (Appendix § 11).

To explain why the leading tail of helium is observed only during Night 1, we can conjecture that, only during that night, the magnetosphere of the planet was remarkably compressed by the interaction with a faster and denser stellar wind streamer with an Alfvénic Mach number $M_{\mathrm{A}} \sim10$. The compression of the magnetosphere by a factor of about two in the direction of the streamer made its internal density increase by approximately the same factor up to a distance of $\sim6 R_{\mathrm{p}}$ in the direction of the streamer. During Night 1, the master-out spectrum was taken before the transit and before the compressed magnetosphere was transiting over the stellar disc. Later, the compressed magnetosphere started to transit over the stellar disc, thus providing an excess absorption with respect to the master-out spectrum. This corresponds to the detection of a leading helium tail before the first contact during Night 1. The angle $\theta_{0}$ formed by the impinging wind streamer with the direction tangent to the planetary orbit depend on the properties of the stellar wind, specifically, its velocity components, and the planet orbital velocity (Llama et al. [52]). It controls the duration of the helium absorption in the leading tail by projecting the radius of the compressed magnetosphere by a $\cos\theta_{0}$ factor along the orbit of the planet.

During Night 4, the leading tail was not observed. We interpret such a result as the stellar wind returning to its quiescent regime ($M_{\mathrm{A}} \sim1$) leading to a magnetosphere with a lower density and a larger extension in the direction of the wind with only a weak shock, if any, at the boundary. On that night, the master-out corresponded to such a lower density in the extended magnetosphere. Its inner part was denser by a factor of about two than its boundary (cf. Appendix § 11), thus explaining the absorption detected in transit during Night 4.

Finally, during Night 2 and Night 3, only the second part of the transit was observed. This corresponded to the transit of the trailing planetary magnetosphere that was not remarkably impacted by the stellar wind, assumed to be in its quiescent regime, and the density of which was higher close to the planet. On those nights, we saw no significant density differences out of transit when using a master-out obtained after the fourth contact because the density in the extended magnetosphere behind the planet was varying slowly due to the large pressure scale height. If we adopt a planet mass of $\sim0.63\,M_{\mathrm{Jup}}$ (SM21), the density gradient across the magnetosphere becomes remarkable, and our interpretation is no longer valid (cf. App. § 11). Adopting the planet mass by L26, our model can also account for the broadening of the neutral helium triplet and its variations from one night to the other (cf. Figure 4), at least qualitatively. We discuss this in App. § 10.

Planetary evolution

The hydro-D simulations with ATES provide us with an estimate of the atmospheric mass loss rate at the present age of the system. In this respect, we have recomputed the stellar age by matching the position of V1298 Tau in the temperature-luminosity diagram with the MESA evolutionary tracks ([14]). To this end, we assumed the stellar mass and effective temperature reported in Table 2, and $L_{\mathrm{bol}} = 0.943 \pm0.018\,L_{\odot}$. These values lead to an age of $12 \pm1$ Myr, which is in agreement with the results of Finociety et al. [30], but significantly younger than the age of $23 \pm4$ Myr obtained by David et al. [21] using magnetic-inhibited convection models of stellar interiors ([27]). In the following, we adopt this younger age. We recommend the reading of Gaidos et al. [34] for a detailed analysis of the uncertainties on the age of this system.

In order to investigate whether photoevaporation may produce relevant variations of the planetary radius and mass, an educated guess of the long-term evolution of the stellar XUV emission is required. Following the approach first developed in Locci et al. [53] and subsequently updated to study the evolution of single systems (e.g., [56, 18, 65]), we have simulated the potential photoevaporation history of V1298 Tau b, from an age of 3 Myr to about 5 Gyr. In our modelling, we adopt the ATES code ([95] (Caldiroli et al. 2021; Caldiroli et al. 2022)), the planetary core-envelope models by Fortney et al. [31] and Lopez & Fortney [54], and the MESA Stellar Tracks ([14]). Initially, we assumed a planetary mass of $\sim0.042\,M_{\mathrm{Jup}}$ (L26) and computed the internal structure of the planet, obtaining a core mass of $8.8\,M_{\oplus}$ and a core radius of $1.7\,R_{\oplus}$. Moreover at the current system age, we obtained an atmospheric mass fraction of 33.8%.

We then evolve the system forward and backward in time, updating the planetary radius to account for natural gravitational contraction and photoevaporation-driven atmospheric mass loss. For the stellar XUV luminosity history, we perform two simulations: the first is based on Penz et al. [77] coupled with the X-ray-to-EUV scaling law by Sanz-Forcada et al. [84] (hereafter P08+SF25), while the second adopts the semi-empirical modeling of Johnstone et al. [42] (hereafter Jo21). In both cases, we estimate a current mass loss rate of $\sim1.45 \times10^{12}\,\mathrm{g\,s^{-1}}$.

The predictions of the two simulations at 3 Myr are similar. In both cases, the planet begins its evolution with a mass of approximately $13.4\,M_{\oplus}$, an atmospheric mass fraction of about 34.2%, and a radius of $11.6\,R_{\oplus}$. Since the P08+SF25 tracks predict a higher XUV luminosity, the corresponding evaporation rates are larger in this scenario, leading to a more substantial mass loss. At 5 Gyr, the P08+SF25 model predicts a planetary radius of approximately $4.5\,R_{\oplus}$ and an atmospheric mass fraction of 13.5%, whereas the Jo21 model provides a radius of $5.2\,R_{\oplus}$ and an atmospheric mass fraction of 20.2%. In Figure 8 we show the temporal evolution of the planetary radius, atmospheric mass fraction, and mass-loss rate. Both evolutionary scenarios indicate that the planet loses more than half of its primordial atmosphere and its radius contracts evolving into a mature sub-Neptune, in agreement with Barat et al. [4] and L26. Alternatively, when we performed the simulation by assuming a mass of $\sim0.63\,M_{\mathrm{Jup}}$ (SM21), we obtained that the deeper gravitational potential ensures that V1298 Tau b remains scarcely affected by photoevaporation induced by high-energy radiation, thus confirming the conclusions of Maggio et al. [56]. Conversely, a direct comparison with the findings of Maggio et al. [56] is not feasible when considering the mass from L26, due to the substantial discrepancy in the assumed mass values.

Conclusions

In this work, we investigate the metastable He I triplet of V1298 Tau b by analysing four transit events observed with the high-resolution GIANO-B spectrograph. Due to the long transit duration ($\sim7$ hours), the dataset exhibits uneven phase coverage and lacks comprehensive out-of-transit baselines, missing contiguous pre- and post-transit observations. This poses significant challenges for characterising a planet orbiting such an active host star. Further complexity arises from the intrinsic blending of the He I triplet components with the nearby photospheric Si I line, due to the fast rotation of the star, as well as significant telluric contamination from OH emission and H$_2$O absorption. Finally, the low planetary $K_{\mathrm{p}}$ complicates the disentanglement of the signals, making it difficult to definitively align the observed features with either the stellar or the planetary rest frames.

We first modelled the RML effect, refining the projected spin-orbit angle to $\lambda= 31 \pm14$ deg. While this value suggests a potential departure from exact alignment, the large uncertainties make it fully consistent with previous literature estimates. We then performed He I transmission spectroscopy, reporting the detection of He I absorption in all four analyzed nights; however, we observe significant night-to-night variability. The measured absorption depths (c) are $8.47 \pm0.71\%$, $3.98 \pm0.40\%$, $1.58 \pm0.77\%$, and $4.49 \pm0.43\%$, with corresponding EWs of $154.0 \pm4.5$ mÅ, $81.4 \pm1.9$ mÅ, $18.8 \pm1.7$ mÅ, and $113.4 \pm2.0$ mÅ for Night 1, Night 2, Night 3, and Night 4, respectively. We need to emphasise that these values are derived from time-averaged in-transit spectra, which were observed over different phase intervals due to the uneven coverage of the dataset. Specifically, for Night 3, the observations only cover the final part of the transit.

In agreement with previous studies (e.g., [103, 34, 46, 1]), the master out-of-transit spectra also exhibit variability, which is driven by the strong activity of the young star.

To disentangle stellar activity from planetary signals, we also analysed the log $R'_{\mathrm{HK}}$ and H$\alpha$ indicators, despite the limited data coverage for Night 3 and Night 4. The variability in the master-out EWs of both He I and H$\alpha$ against log $R'_{\mathrm{HK}}$ underscores the significant impact of the active regions on these lines. In the H$\alpha$ transmission spectra, we observed contrasting behaviours: Night 2 exhibits a pronounced in-transit pseudo-emission “bump” (likely attributable to the occultation of a large stellar starspot or spot complex, as indicated by NOT simultaneous observations), whereas Night 1 shows only a marginal absorption signal. Although these H$\alpha$ features are clearly dominated by stellar activity, which likely contaminates the He I measurements, a key distinction emerges. While the H$\alpha$ indicator displays opposite trends across different epochs, a behaviour further supported by the Na I analysis, which reveals an emission signal during Night 2 but is affected by significantly larger uncertainties in Night 1; the He I triplet shows an absorption signal across four independent measurements spanning a four-year baseline, albeit with varying depths. This may support a planetary origin for the He I signal, although the underlying stellar activity inevitably modulates its observed depth.

This interpretation is supported by our ATES simulations, which demonstrate that stellar XUV flux variations significantly alter the absorption depth. We show that the response of the He i triplet to XUV radiation is highly sensitive to planetary mass, shifting from an over-ionization regime, and thus reduced absorption, for the higher-mass scenario (SM21), to an optically thick outflow dominated by neutral and singly ionized species, resulting in enhanced absorption, for the lower-mass case (L26).

We attribute the blueshifted leading tail, seen only during Night 1, to a compression of the planetary magnetosphere likely triggered by a fast and dense stellar wind streamer. This explains why the tail vanished during the less active Night 4, despite the pre-transit coverage that would have allowed its detection.

Finally, we investigated the evolution of V1298 Tau b to assess whether photoevaporation induces significant variations in its mass and radius. Assuming the lower mass by L26, the planet loses more than half of its primordial atmosphere and its radius contracts, depending on the model, up to around 4.5 $R_{\oplus}$, evolving into a mature sub-Neptune. Alternatively, adopting the higher mass from SM21, the deeper gravitational potential ensures that V1298 Tau b remains largely unaffected by photoevaporation.

Our findings highlight that without continuous activity monitoring and full transit coverage for all epochs, establishing a robust baseline to disentangle stellar contamination from the true planetary signal remains a formidable challenge. Ultimately, for a star as variable as V1298 Tau, co-adding multiple nights with uneven coverage remains problematic. A complete transit with simultaneous multi-wavelength monitoring represents the only definitive path to fully resolve this complex “puzzle", requiring dedicated proxies for the stellar wind and XUV emission. In this context, future UV-to-near-infrared space concepts like the Pollux instrument (Neiner et al. [66]) for the Habitable Worlds Observatory (HWO, Dressing et al. [23]) will enable a simultaneous evaluation of stellar high-energy behaviour and the planetary He i response, ultimately unveiling the true nature of signals detected around such active host stars.

Acknowledgements. We thank the referee for helpful feedback that improved our work. This work is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF at the Spanish Observatorio Roque de los Muchachos of the IAC in the frame of the program Global Architecture of the Planetary Systems (GAPS). We thank A. Vidotto, and A. Presa for the interesting discussion on star-planet magnetic interactions in photoevaporating exoplanets. Based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofísica de Canarias. The NOT data were obtained under program ID 62-009; P.I.: Luigi Mancini. The ALFOSC instrument is provided by the Instituto de Astrofísica de Andalucía (IAA) under a joint agreement with the University of Copenhagen. We 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. We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through projects PID2021-125627OB-C32, CEX2021-001131-S, PID2022-137241NB-C42, PID202241216NB-I00 and PID2024-158486OB-C32. This work is also supported by the ERC Grant (ERC Advanced Grant SPEAR, GA 101200674), funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency (ERCEA). Neither the European Union nor the granting authority can be held responsible for them. JS acknowledges support from STFC under grant number ST/Y002563/1. FB acknowledges support from INAF Go Large Grant 2024 “Strengthening pylons for BRIDGES”. The authors acknowledge the support from ASI-INAF agreement 2021-5-HH.0. CDM acknowledges the support 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”.

References

1 INAF – Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025, Pino Torinese, Italy

2 Department of Astronomy, University of Texas at Austin, 2515 Speedway, Austin, TX 78712, USA

3 INAF – Osservatorio Astronomico di Palermo, P.zza Parlamento 1, 90134 Palermo, Italy

4 INAF – Osservatorio Astrofisico di Catania, Via S. Sofia 78, 95123, Catania, Italy

5 Dipartimento di Fisica e Astronomia “Ettore Majorana”, Università di Catania, Via S. Sofia 64, 95123 Catania, Italy

6 Scuola Superiore di Catania, Università di Catania, Via Valdisavoia 9, 95123 Catania, Italy

7 Como Lake Center for Astrophysics (CLAP), DiSAT, Università degli Studi dell’Insubria, Via Valleggio 11, 22100 Como, Italy

8 Dipartimento di Fisica e Astronomia, Università degli Studi di Padova, Vicolo dell’Osservatorio 3, 35122 Padova, Italy 2

9 INAF– Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy

10 Department of Physics, University of Rome “Tor Vergata”, Via della Ricerca Scientifica 1, 00133, Roma, Italy

11 INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807, Merate (LC), Italy

12 INAF – Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, 80131, Naples, Italy

13 INAF – Osservatorio Astronomico di Trieste, via Tiepolo 11, 34143, Trieste, Italy

14 INAF – Fundación Galileo Galilei, Rambla José Ana Fernandez Pérez 7, 38712 Breña Baja, Tenerife, Spain

15 Department of Physics, University of Turin, Via Pietro Giuria 1, 10125 Torino, Italy

16 Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042, Graz, Austria

17 Department of Astronomy, University of Michigan, 1085 S University, Ann Arbor, MI 48109, USA

18 Center for Astrophysics and Space Science, NYUAD Research Institute, New York University Abu Dhabi, PO Box 129188, Abu Dhabi, UAE

19 Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain

20 Instituto de Astrofísica de Canarias (IAC), 38205 La Laguna, Tenerife, Spain

21 Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain

22 Centro de Astrobiología, CSIC-INTA, Camino bajo del Castillo s/n, 28692 Villanueva de la Cañada, Madrid, Spain

23 Astrophysics Group, Keele University, Staffordshire, ST5 5BG, United Kingdom

24 Centro de Astronomía de la Universidad de Antofagasta, Angamos 601, Antofagasta, 1270300, Chile

Additional figures and tables

This section provides additional figures and tables that complement the analysis presented in the main text.

Reference	$M_{\star}$ ($M_{\odot}$)	$R_{\star}$ ($R_{\odot}$)	$v\sin i_{\star}$ (km s$^{-1}$)	$P$ (d)	$R_{\mathrm{p}}/R_{\star}$	$a/R_{\star}$	$b$	$i$ (deg)
D19a	$1.10 \pm0.05$	$1.305 \pm0.070$	$23 \pm2$	$24.13889^{+0.00043}_{-0.00044}$	$0.07111^{+0.00117}_{-0.00061}$	$28.7^{+1.5}_{-2.3}$	$0.23^{+0.16}_{-0.15}$	$89.51^{+0.33}_{-0.38}$
D19a$^{(a)}$	$1.099^{+0.049}_{-0.049}$	$1.314^{+0.052}_{-0.064}$	$\ldots$	$24.13861^{+0.00102}_{-0.00090}$	$0.0713^{+0.0021}_{-0.0020}$	$27.6^{+1.5}_{-1.0}$	$0.29^{+0.16}_{-0.19}$	$89.35^{+0.41}_{-0.33}$
D19b	$1.101^{+0.049}_{-0.051}$	$1.345^{+0.056}_{-0.051}$	$\ldots$	$24.1396 \pm0.0018$	$0.0700 \pm0.0023$	$27.0 \pm1.1$	$0.46^{+0.13}_{-0.24}$	$89.00^{+0.46}_{-0.24}$
G21	$1.145 \pm0.075$	$1.317 \pm0.027$	$23 \pm2^{(b)}$	$\ldots$	$\ldots$	$27.0 \pm1.1^{(b)}$	$\ldots$	$\ldots$
V21	$1.104^{+0.046}_{-0.045}$	$1.333^{+0.030}_{-0.024}$	$\ldots$	$\ldots^{(c)}$	0.0128^{+0.0137}_{-0.0089}^{(d)}	$\ldots$	$0.439^{+0.041}_{-0.016}$	$\ldots$
F21$^{(e)}$	1.101^{+0.049}_{-0.051}^{(f)}	1.345^{+0.056}_{-0.051}^{(f)}	$23.65^{+0.87}_{-1.17}$	$\ldots$	$\ldots$	$\ldots$	$\ldots$	$\ldots$
F22	$1.095^{+0.049}_{-0.047}$	$1.33^{+0.04}_{-0.03}$	$\ldots$	$24.1315^{+0.0033}_{-0.0034}$	$0.0636 \pm0.0018$	$\ldots$	$0.45^{+0.05}_{-0.04}$	$\ldots$
J22	$1.135 \pm0.047$	$1.407 \pm0.018$	$24.77 \pm0.19$	$24.141341 \pm0.000023$	$0.071998 \pm0.00062$	$26.06 \pm0.46$	$0.564 \pm0.026$	$88.759 \pm0.077$
SM21	$1.17 \pm0.06$	$1.28 \pm0.07$	$23.8 \pm0.5$	$24.1399 \pm0.0015$	$0.0698 \pm0.0024$	$\ldots$	$\ldots$	$>88.7$
F23	$1.26 \pm0.06$	$1.43 \pm0.03$	$24.9 \pm0.5$	$24.1315^{(g)}$	$\ldots$	$\ldots$	$\ldots$	$\ldots$
S23	$1.157^{+0.057}_{-0.058}$	$1.355^{+0.032}_{-0.030}$	$\ldots$	$24.140410 \pm0.000022$	$0.0673 \pm0.0017$	$27.25^{+0.40}_{-0.49}$	$0.451^{+0.033}_{-0.030}$	$89.052^{+0.072}_{-0.087}$
L26$^{(h)}$	$1.10 \pm0.05^{(i)}$	$1.32 \pm0.05^{(i)}$	$\ldots$	$24.140006 \pm0.000017$	$\ldots$	$\ldots$	$\ldots$	$\ldots$
L26$^{(j)}$	$\ldots$	$\ldots$	$\ldots$	$24.140420 \pm0.000012$	$0.0661^{+0.0019}_{-0.0017}$	$\ldots$	$0.26 \pm0.17$	$\ldots$
L26$^{(k)}$	$\ldots$	$\ldots$	$\ldots$	$24.13934 \pm0.00030$	$\ldots$	$\ldots$	$\ldots$	$\ldots$

Notes. $^{(a)}$ Recommended values based on simultaneous variability and transit fits; $^{(b)}$ Adopted from D19a and D19b; $^{(c)}$ No period value reported, V21 rather mention that ‘ We fixed the mean ephemerides to the best-fit solution from Livingston et al. (2021, in preparation)...’, the paper being eventually published as L26; $^{(d)}$ Typo in Table 2 of V21, in which the radius ratio value of V1298 Tau c is used. The likely correct value for V1298 Tau b is $R_{\mathrm{p}}/R_{\star} = 0.0642^{+0.0048}_{-0.0047}$, which is erroneously reported for V1298 Tau d; $^{(e)}$ The F21 work focuses on V1298 Tau c, but it is relevant for the independent estimate of $v\sin i_{\star}$ for the primary; $^{(f)}$ Adopted from D19b; $^{(g)}$ Adopted from F22; $^{(h)}$ Based on analytic TTV modeling reported in Table 1 and Supplementary Table 3 of L26; $^{(i)}$ Adopted from D19a; $^{(j)}$ Based on joint modeling of the K2 and TESS photometry reported in Extended Data Table 1 of L26; $^{(k)}$ Based on results of N-body modeling reported in Extended Data Table 2 of L26.

	SM21: $M_{\mathrm{p}} = 0.63M_{\mathrm{Jup}}$			L26: $M_{\mathrm{p}} = 0.042M_{\mathrm{Jup}}$
$F_{\mathrm{XUV}}/F_{0}$	$\dot{M}$ (g s$^{-1}$)	$c$ (%)	EW(He I) (mÅ)	$\dot{M}$ (g s$^{-1}$)	$c$ (%)	EW(He I) (mÅ)
0.5	$8.13 \times10^{9}$	5.90	29.43	$2.24 \times10^{12}$	15.60	154.9
1.0	$1.62 \times10^{10}$	4.78	25.90	$3.98 \times10^{12}$	20.58	231.6
2.0	$3.16 \times10^{10}$	3.72	21.98	$6.92 \times10^{12}$	26.00	328.3

NOT light curve

A transit of V1298 Tau b was photometrically monitored with the NOT telescope on December 7, 2020, simultaneously with GIARPS observations on Night 2. We used the Alhambra Faint Object Spectrograph and Camera (ALFOSC) in its imaging mode and selected the engineered diffuser #1 (von Essen et al. [104]). The observations were performed through a $I$-band (#12) filter. The 6.5′ × 6.5′ FOV was suitably chosen to include three good reference stars, namely HD 284152, HD 284153, and HD 284154. The observations started at 19:35 UT and ended at 04:39 UT. The night started with a few clouds, which cleared by 23:00 UT. The airmass varied from 2.28 to 1.01 and then to 1.98. The moon illuminated fraction at the midpoint of exposure was 50%. The exposure time was 35 s at the beginning of the observations and was gradually reduced to 6 s. A total of 918 images were collected.

The data were analysed using the defot pipeline (Southworth et al. [91, 90]). Following a standard approach, which has been used in many previous works (e.g., Southworth et al. [93, 92]), the raw science images were initially calibrated using master-bias and master-flat frames. These calibration frames were created by median-combining individual bias and sky flat-field images captured just after sunset. Subsequently, we performed aperture photometry by selecting a set of five comparison stars and cross-correlating each image in the time series with a reference frame. We tried different aperture sizes and retained those (22 px, 45 px, 65 px) that gave the light curve with the smallest scatter. Finally, differential-magnitude light curves were extracted relative to each comparison star. Instrumental magnitudes were converted into differential-magnitude light curves, which were normalized to zero out of transit. This normalization

was achieved by fitting second-order polynomials to the out-of-transit data. The differential magnitudes were calculated relative to a weighted ensemble of five comparison stars that we already mentioned. To minimize the scatter in the out-of-transit data, we simultaneously optimized the comparison star weights and the polynomial coefficients. Timestamps were converted from JD$_{\mathrm{UTC}}$ to BJD$_{\mathrm{TDB}}$ following the recommendations of Eastman et al. [26].

Due to adverse weather conditions in the early part of the night, the first data points of the light curve resulted in being very scattered. Therefore, we were forced to exclude the first 222 points from the analysis and considered only those starting from 22:32 UT, thereby losing all the ingress data points. The final light curve, composed of 696 data points, presents a clear bump that peaks close to the transit midpoint and lasts about half the entire transit. We interpreted this anomaly as the consequence of a starspot complex on the stellar photosphere that was occulted by the planet during its transit; top panel of Fig. B.1.

Following the approach of previous works in which similar cases have come across (e.g., Mohler-Fischer et al. [63]; Mancini et al. [59, 58, 60]; Biagiotti et al. [7]), we used the prism and gemc codes⁴ (Tregloan-Reed et al. [101, 100]), which allowed us to model the light curve by fitting both the transit and the starspot-crossing event simultaneously, thanks to a Bayesian approach and the use of genetic algorithms (Ter Braak [99]). Utilizing a quadratic limb-darkening (LD) law and assuming a circular orbit, prism+gemc models light curves by fitting the following parameters: the sum of the fractional stellar and planetary radii ($r_{\mathrm{p}} + r_{\star}$, $r_{\mathrm{p}}$ and $r_{\star}$ are defined as the planetary radius $R_{\mathrm{p}}$ and the stellar radius $R_{\star}$ scaled by the semi-major axis $a$, respectively), their ratio $r_{\mathrm{p}}/r_{\star}$, the orbital inclination $i$ and the LD coefficients. It also fits the mid-transit time, the longitude and latitude of the starspot centre ($\phi$, $\theta$), its radius $r_{\mathrm{s}}$, and its contrast $\rho$, the latter defined as the surface brightness ratio between the starspot and the surrounding photosphere.

We started the modelling of the light curve by fixing the orbital period at $P_{\mathrm{orb}} = 24.14141$ days and the mid-transit time at $T_0 = 2459\,191.4932383439$ (private communication), and assuming the starspot as a spherical cap. By imposing this geometry, a large ($25^{\circ} \pm4^{\circ}$), high-latitude starspot at $48^{\circ} \pm5^{\circ}$ is needed to reproduce such a long-lasting bump. Although it may be a long active region located along the transit chord, such a large starspot is very suspicious because it would produce an apparent blueshift as it moves across the half of the stellar disk that is moving away from the observer. In reality, prism+gemc scans the stellar surface within the transit chord, which, when using circular starspots, causes (in some cases) large starspots that graze the transit chord. Therefore, we tried a different approach and attempted to model an elongated series of small starspots instead of a single starspot. We imposed a starspot on the pixels within the transit cord and apply a contrast gradient to the starspot. The gradient follows a cosine function where the starspot centre has a contrast given by the best-fitting parameter and the starspot edge has a contrast of 1.0. The equation we used is

where $\xi$ is the angle between the spot centre and the pixel being evaluated and ranges between 0 (starspot centre) and 90 degrees (starspot edge). The best-fitting parameters indicate a starspot complex within the transit cord, whose sizes are approximately 29 degrees in length and 9 degrees in height, bottom panel of

Parameter	Single starspot	Starspot complex
$R_{\mathrm{p}}/R_{\star}$	$0.0823 \pm0.0020$	$0.0807 \pm0.0080$
$a/R_{\star}$	$33.8 \pm2.5$	$33.4 \pm7.0$
$i$ (deg)	$89.23 \pm0.18$	$89.58 \pm0.39$
$u_{\mathrm{A}}$	$0.93 \pm0.17$	$0.990 \pm0.047$
$v_{\mathrm{A}}$	$0.07 \pm0.16$	$0.010 \pm0.045$
$\theta$ (deg)	$1.11 \pm0.45$	$1.13 \pm0.36$
$\phi$ (deg)	$48.2 \pm4.8$	$9.3 \pm9.9$
$r_{\mathrm{spot}}$ (deg)$^{(a)}$	$24.6 \pm3.8$	$29 \pm10$; $9.24 \pm0.92$
$\rho_{\mathrm{spot}}$ (deg)	$0.338 \pm0.050$	$0.270 \pm0.070$
$T_{\mathrm{spot}}$ (K)	$3890 \pm136$	$3711 \pm204$

Fig. B.1. The parameters of both fits (single starpsot and starspot complex) are reported in Table B.1.

NaI analysis

For the Na i analysis, we employed the SLOPpy (Spectral Lines Of Planets with python) pipeline (Sicilia et al. [87]), a standard user-friendly and reliable tool optimised for spectral reduction and the extraction of planetary transmission spectra in the VIS obtained from high-resolution observations. Although the main body of this work investigates a total of four observational nights, the sodium line analysis was exclusively performed on the first two nights, which are the only ones featuring high-resolution data from the HARPS-N spectrograph. During the data reduction and selection process, one out-of-transit spectrum from Night 1 was discarded due to an extremely low signal-to-noise ratio (S/N $\approx2.6$) to prevent the introduction of severe systematics or spurious noise into the final transmission spectrum.

The spectra were corrected for stellar spectroscopic effects that typically contaminate transit observations, namely the Rossiter-McLaughlin (RM) effect and center-to-limb variations (CLV). These corrections were implemented directly through SLOPpy, utilizing synthetic spectral models generated via PySME (Wehrhahn [106]). The stellar atmosphere models were computed using the line list from the third version of the Vienna Atomic Line Database (VALD3$^{4}$) combined with MARCS stellar atmosphere models, adopting non-LTE grids.

The final transmission spectrum, corrected for stellar systematics, was modeled using a Markov Chain Monte Carlo MCMC sampling procedure within SLOPpy, applying a prior on the theoretical planetary orbital radial velocity semi-amplitude ($K_{\mathrm{p}}$). The model assumes a double Gaussian profile for the Na i atmospheric absorption or emission lines and a flat spectrum otherwise. The free parameters allowed to vary during the fitting process include $K_{\mathrm{p}}$, which is required to model the atmospheric lines in the planetary rest frame; the line contrasts ($c_{\mathrm{D1}}$ and $c_{\mathrm{D2}}$); the FWHM of the double Gaussian profile; and the residual radial velocity of the atmospheric wind ($v_{\mathrm{wind}}$), defined as the shift relative to the transition in the planetary rest frame.

Parameter	Night 1	Night 2
$c_{\mathrm{D1}}$ (%)	$+1.3^{+3.8}_{-2.2}$	$+4.4^{+0.5}_{-0.7}$
$c_{\mathrm{D2}}$ (%)	$-3.5^{+7.3}_{-3.2}$	$+5.0^{+3.6}_{-0.6}$
$\mathrm{FWHM}_{\mathrm{shared}}$ (km s$^{-1}$)	$5.6^{+3.5}_{-5.0}$	$5.8^{+0.5}_{-0.6}$
$v_{\mathrm{wind}}$ (km s$^{-1}$)	$+7.9^{+3.5}_{-9.6}$	$-22.4^{+0.35}_{-0.24}$
$K_{\mathrm{p}}$ (km s$^{-1}$)	$76.7^{+10.8}_{-9.3}$	$52.6^{+10.0}_{-9.3}$
jitter	$0.0200 \pm0.0008$	$0.0193 \pm0.0004$

We assumed a shared FWHM and $v_{\mathrm{wind}}$ for both D lines. To account for any unmodeled instrumental or environmental systematics, such as variations in seeing, a free jitter parameter ($\mathrm{jitter}_{0}$) was also included in the fitting procedure. The median values of the resulting posterior distributions were adopted as the best-fit parameters, with error bars corresponding to the $1\sigma$ statistical uncertainties at the 15th and 84th percentiles. The final transmission spectra are shown in Figure C.1.

The MCMC fit results for the two nights reveal a significant discrepancy in signal quality and parameter constraints, as summarised in Table C.1. For Night 1, the retrieved parameters exhibit exceptionally large uncertainties. Although the nominal best-fit profile suggests a disparate behaviour between the two lines—with D2 in absorption and D1 in emission—both contrasts are fully compatible with zero within their error bars. Consequently, the data from Night 1 are inconclusive and can neither confirm nor rule out an atmospheric signal. Conversely, Night 2 yields a much cleaner and tightly constrained signal, characterised by a distinct double emission profile in both components of the doublet. The positive contrast values robustly confirm an emission signal well above the noise floor, accompanied by a strong blue-shift.

We note that this pronounced discrepancy between the two observational epochs closely follows the variability observed in the H$\alpha$ line. Conversely, the He I triplet displays a consistent behaviour, with transit absorption being detected in both nights.

Modelling an absorption leading tail in V1298 Tau b

A young planet such as V1298 Tau b is highly inflated and is probably characterized by a strong hydromagnetic dynamo action because of the intense convective flows in its interior. The timescale for the synchronization of its rotation with the orbital motion, $\tau_{\mathrm{sync}}$, can be estimated by means of the tidal model of Leconte et al. [49], adopting $k_2\Delta t = 3/(2nQ^\prime)$, where $k_2$ is the planet Love number of order 2, $\Delta t$ the tidal time lag, $n$ the orbital mean motion, and $Q^\prime$ the modified tidal quality factor of the planet that we assume similar to that of Jupiter, $Q^\prime\sim10^5$. This gives a timescale of $\tau_{\mathrm{synch}} \sim0.3$ Myr, adopting a gyration radius similar to Jupiter’s. This is much shorter than the age of the star estimated to be between 10 and 30 Myr, therefore, the planet rotation period can be considered to be equal to its orbital period. In spite of that, the planetary dynamo can be remarkably strong owing to the large internal thermal energy of the planet responsible for its highly inflated radius ([81, 15, 107]). As a consequence, we believe that a planetary dipole field with an intensity of $B_{\mathrm{p0}} = 25$ G at the poles can be a reasonable assumption. On the other hand, the stellar dipolar magnetic field at the stellar poles can be assumed to be $B_{\mathrm{s0}} = 200$ G, according to the Zeeman Doppler Imaging maps of V1298 Tau (cf. Finociety et al. [30]).

The average radius of the planetary magnetosphere can be estimated by assuming pressure equilibrium at the boundary of the magnetosphere in front of the stellar wind flow. Specifically, we assume that the magnetic pressure of the planetary field at the boundary of radius $r_{\mathrm{M}}$, measured from the centre of the planet, is equal to the sum of the thermal pressure of the stellar wind, $p_{\mathrm{w}}$, the magnetic pressure of the stellar wind field, $B_{\mathrm{w}}$, and the ram pressure of the stellar wind having a density $\rho_{\mathrm{w}}$ and a relative velocity $u_{\mathrm{w}}$ at the distance $a$ of the planet

where $\mu$ is the magnetic permeability of the plasma that can be assumed to be equal to that of the vacuum. Both the stellar and the planetary fields can be assumed to scale as a dipole field with the distance from the centre of the star and of the planet, respectively. This gives $B_{\mathrm{p}}(r_{\mathrm{M}}) = B_{\mathrm{p0}}(r_{\mathrm{M}}/R_{\mathrm{p}})^{-3}$ and $B_{\mathrm{s}}(a) = B_{\mathrm{s0}}(a/R_{\mathrm{s}})^{-3}$, where $R_{\mathrm{p}}$ and $R_{\mathrm{s}}$ are the radius of the planet and of the star, respectively.

We assume that the thermal pressure of the stellar wind is dominated by the pressure of the magnetic field of the wind or by the ram pressure, so that $p_{\mathrm{w}}(a) \ll\max[B_{\mathrm{s}}^2(a)/2\mu,\rho_{\mathrm{w}}(a)u_{\mathrm{w}}^2(a)/2]$, the planet being relatively close to its host star. The kinetic ram pressure of the wind can be expressed by introducing the Alfven Mach number $M_{\mathrm{A}}$ of the wind at the distance of the planet

Substituting into Eq. (D.1), we find

We can assume that the stellar wind of V1298 Tau is highly inhomogeneous and remarkably time dependent with streamers of widely different velocity and density. Considering a typical situation where the planet is close to the transition between the sub-Alfvenic and the super-Alfvenic regimes of the stellar wind, $M_{\mathrm{A}} \sim1$, and Eq. (D.3) gives $r_{\mathrm{M}}/R_{\mathrm{p}} \sim12.2$. On the other hand, when the planet is reached by a fast and dense streamer of the stellar wind with $M_{\mathrm{A}} \sim10$, the radius of the hemisphere of the magnetosphere impacted by the streamer is compressed to $r_{\mathrm{M}}/R_{\mathrm{p}} \sim6.35$, according to Eq. (D.3). A simple sketch of the compressed planetary magnetosphere is plotted in Fig. D.1 where the star is located at a large distance, $a/R_{\mathrm{p}} \sim400$, approximately along the positive $x$-axis.

The planet V1298 Tau b is subject to an intense photoevaporation due to the strong XUV radiation coming from its young and active star. The ionized evaporation flow moves along the magnetic field lines of the planetary magnetosphere and fills it with plasma. Considering a typical temperature of $\sim10^{4}$ K, and the radius and mass of the planet as determined by L26, the pressure scale height is $H_{\mathrm{p}} \sim1.6 R_{\mathrm{p}}$, suggesting that the evaporated plasma fills the magnetosphere up to its outer boundary considered above. Specifically, for that temperature, the density variations between the planet surface and the limit of the magnetosphere are within a factor of about two, while in a planet with a stronger surface gravity remarkably larger variations will occur (Sect. E).

The slow rotation of the planet induced by tidal synchronization implies that the centrifugal force is negligible in comparison to gravity and plasma pressure in the planetary magnetosphere, contrary to the case of Jupiter that is a fast rotator in this context. Quantitatively, even at the outer boundary of the magnetosphere of V1298 Tau b, $r = r_{\mathrm{M}} = 12.2 R_{\mathrm{p}}$, the ratio of the centrifugal to the gravitational acceleration is 0.69 and scales as $r^{3}$ with the distance from the centre of the planet within the magnetosphere.

The mean density of the magnetospheric plasma in front of the planet, in the hemisphere impacted by the wind streamer, is approximately proportional to the inverse of the volume of the hemisphere itself. Assuming that the cross section of the hemisphere orthogonal to the wind streamer remains approximately constant, that volume is proportional to $r_{\mathrm{M}}$. Therefore, it decreases by a factor of $\sim2$ when a wind streamer with $M_{\mathrm{A}} \sim10$ impacts onto the magnetosphere of our planet producing an increase in the plasma mean density by a factor of 2. That increase in the density of the magnetospheric plasma in front of the planet can be invoked to interpret the overdensity of the neutral helium as observed in the tail preceding the planet before its transit.

The streamer is not moving radially away from the star, rather it forms an angle with the radial direction that depends on the conservation of the angular momentum of the streamer’s plasma, its radial velocity, and the orbital velocity of the planet. As a consequence, the compressed magnetosphere can lead the planet along its orbit, thus producing an absorption in the helium triplet that is observable before the ingress of the planet on the stellar disc.

The present magnetospheric model assumes a pressure balance at the boundary of the planet magnetosphere and a negligible penetration of the stellar wind plasma inside the magnetosphere itself. On the other hand, if the wind plasma penetrates into the magnetosphere, we have to model the magnetospheric boundary as a strong MHD perpendicular shock because of the relatively high Alfven Mach number ($M_{\mathrm{A}} \sim10$). In that case, the density inside the magnetosphere can increase up to a factor of $\sim4$ ([79]), thus further enhancing the absorption in the leading tail over the effect of the magnetospheric compression considered above, while the velocity of the plasma past the shock can lead to an additional Doppler shift and an increase in the width of the absorption of the helium triplet owing to its velocity component in the direction of the observer and the extension of the flow to the hemisphere of the magnetosphere past the planet.

The neutral helium atoms are formed by recombination of the single-ionized helium, therefore, their density mirrors that of those ions. Nevertheless, once neutral helium atoms form, their dynamics is no longer constrained by the magnetic field and they can freely flow evaporating away from the planet. The helium triplet line blueshifts and widths reflect such an evaporating motion (cf. Figure 4). The blueshift is produced by the component of the evaporation flow towards the observer that is at its maximum during Night 1 when the evaporation is likely faster because of the higher level of stellar activity. The sonic radius in the case of our very inflated planet is only $\sim0.3 R_{\mathrm{p}}$ for a temperature of the evaporation flow of $10^{4}$ K implying that the wind flow becomes remarkably supersonic within the boundary of the magnetosphere. However, even in this regime, the line broadening that can be explained by the evaporating flow is still of the order of $\pm c_{\mathrm{s}}$, where $c_{\mathrm{s}} \sim10\ \mathrm{km\ s}^{-1}$ is the sound speed in the evaporating flow. Nevertheless, the line broadening during Night 1 is larger by a factor of at least two, suggesting that there is a component of the neutral helium absorption produced by atoms originated by charge exchange between the helium ions of the stellar wind, having a remarkably larger velocity, and the neutral helium atoms at the boundary of the planetary magnetosphere. This combination of absorption due to neutralized stellar wind ions and evaporating helium atoms from the planet can account, at least qualitatively, for the broad line profile observed during Night 1 corresponding to a velocity in excess of $\sim40\ \mathrm{km\ s}^{-1}$, and its remarkable blueshift (cf. Figure 4, top panel). On the other hand, the reduced blueshifts and/or smaller line widths during the other nights (cf. Figure 4, lower panels) can be accounted for by a significantly slower stellar wind. It is interesting to note that not only the average stellar wind velocity is relevant here to account for the blueshift of the line centroids, but also its velocity dispersion close to the planet, owing to the turbulence in the wind itself, that plays a role in the line broadening because the wind neutrals collide with the helium atoms of the planetary evaporating flow and exchange momentum with them.

Plasma stratification inside the planetary magnetosphere

The plasma evaporating from the planet and filling its magnetosphere reaches a condition of mechanical equilibrium described by the equation

where $p$ is the plasma pressure, $\rho$ its density, $\Phi$ the gravitational potential, $\mathbf{J}$ the current density, and $\mathbf{B}$ the magnetic field. The centrifugal force can be neglected because the planet is rotating slowly (cf. Appendix § 10). Therefore, the gravitational potential can be written as

where $G$ is the gravitation constant, $M_{\mathrm{p}}$ the mass of the planet, and $r$ the radial distance from the centre of the planet. We can project Eq. (E.1) along any magnetic field line having a unit vector tangent $\hat{\mathbf{s}}$, that is, $\mathbf{B} = B \hat{\mathbf{s}}$, yielding

The scalar product of the unit vectors $\hat{\mathbf{r}} \cdot\hat{\mathbf{s}} = \cos\psi$, where $\psi$ is the angle between the tangent to the magnetic field line and the radial direction. Hence $dr = ds \cos\psi$ and Eq. (E.3) becomes

Adopting an isothermal atmosphere with a uniform temperature $T$, the ideal gas law gives $p = \tilde{R}T\rho/\tilde{\mu}$, where $\tilde{R}$ is the gas constant and $\tilde{\mu}$ the mean molecular weight of the plasma. For a fully ionized, solar composition plasma, $\tilde{\mu} \sim0.6$. Equation (E.4) can be integrated analytically making use of the ideal gas law and introducing the pressure scale height $H_{\mathrm{p}} \equiv\tilde{R}T/(\tilde{\mu}g)$, where $g = GM_{\mathrm{p}}/R_{\mathrm{p}}^2$ is the acceleration of gravity at the surface of the planet

where $\rho_0$ is the density at the surface of the planet, where $r = R_{\mathrm{p}}$, and that coincides with the base of the magnetosphere in our simple model.

A lower limit for the density in the magnetosphere can be found from Eq. (E.5) in the limit $r \to\infty$, obtaining $\rho> \rho_0\exp(-R_{\mathrm{p}}/H_{\mathrm{p}})$. In the case of V1298 Tau b, assuming a solar composition plasma and an isothermal magnetosphere with $T = 10^{4}\ \mathrm{K}$, we obtain $H_{\mathrm{p}} \sim1.6R_{\mathrm{p}}$ giving a lower limit for the density inside the magnetosphere of $0.53\rho_0$. In other words, the density does not vary by more than a factor of $\sim2$ inside the planetary magnetosphere in our isothermal model. This is essentially a consequence of the high planet radial inflation yielding a very small surface gravity and therefore a large pressure scale height. On the other hand, for a planet having a stronger surface gravity $g$ and the same temperature of the magnetospheric plasma, the pressure scale height $H_{\mathrm{p}} \propto g^{-1}$ will be smaller and the reduction of the density at large distance from the planet surface will be remarkably larger because the limit density decreases as $\rho_0\exp(-R_{\mathrm{p}}/H_{\mathrm{p}})$.