Monitoring the young planet host V1298 Tau with SPIRou: planetary system and evolving large-scale magnetic field

We report results of a spectropolarimetric monitoring of the young Sun-like star V1298~Tau based on data collected with the near-infrared spectropolarimeter SPIRou at the Canada-France-Hawaii Telescope between late 2019 and early 2023. Using Zeeman-Doppler Imaging and the Time-dependent Imaging of Magnetic Stars methods on circularly polarized spectra, we reconstructed the large-scale magnetic topology of the star (and its temporal evolution), found to be mainly poloidal and axisymmetric with an average strength varying from 90 to 170 G over the ~3.5 years of monitoring. The magnetic field features a dipole whose strength evolves from 85 to 245 G, and whose inclination with respect to the stellar rotation axis remains stable until 2023 where we observe a sudden change, suggesting that the field may undergo a polarity reversal, potentially similar to those periodically experienced by the Sun. Our data suggest that the differential rotation shearing the surface of V1298 Tau is about 1.5 times stronger than that of the Sun. When coupling our data with previous photometric results from K2 and TESS and assuming circular orbits for all four planets, we report a $3.9\sigma$ detection of the radial velocity signature of the outermost planet (e), associated with a most probable mass, density and orbital period of $M_e=0.95^{+0.33}_{-0.24} \ \rm M_{\rm jup}$, $\rho_e=1.66^{+0.61}_{-0.48}$ $\rm g\,cm^{-3}$ and $P_e=53.0039\pm0.0001 \ \rm d$, respectively. For the 3 inner planets, we only derive 99\% confidence upper limits on their mass of $0.44\ \rm M_{\rm jup}$, $0.22\ \rm M_{\rm jup}$ and $0.25\ \rm M_{\rm jup}$, for b, c and d, respectively.


INTRODUCTION
The detection and characterization of exoplanets around young Sunlike stars is key for deriving observational constraints on theoretical models of the early phases of star / planet formation.In particular, it is essential to estimate the radii and masses, and thereby the bulk densities, of young planets to refine the mass-radius relations for planets younger than 30 Myr (Mann et al. 2016).In addition, characterizing the orbital configuration of young close-in giants will bring new clues on the inward migration that giant planets are subject to as they form within the protoplanetary discs of their host stars, and thereby on the impact migrating giant planets have on the dynamical architecture of early planetary systems (Baruteau et al. 2014).However, only few such planets have yet been detected, partly due to the high activity levels that Pre-Main Sequence (PMS) stars ex-★ E-mail: benjamin.finociety@irap.omp.euhibit, especially due to their short rotation period.Such active stars indeed generate photometric and radial velocity (RV) signatures that are much larger than those expected from transiting planets, even in the case of close-in giant ones.To date, only 9 planets younger than 25 Myr have well-measured radii thanks to the detection of photometric transits whose depths are typically 10 times smaller than the photometric variations caused by activity: the Super-Neptune K2-33 b (David et al. 2016); the two warm Neptunes AU Mic b and c (Plavchan et al. 2020;Martioli et al. 2021;Szabó et al. 2022); the four Jupiter-size planets V1298 Tau b, c, d and e (David et al. 2019b;Feinstein et al. 2022); the hot Jupiter HIP 67522 b (Rizzuto et al. 2020); and the warm Jupiter HD 114082 b (Zakhozhay et al. 2022).Constraining the mass of these planets is important to refine the mass-radius relations at such early stages of evolution but also to better characterize these young planets in terms of density and composition.However, detecting these planets through velocimetric monitoring is still challenging given that the amplitude of the RV fluctuations induced by the star itself can reach more than ∼10× the amplitude of the planet signatures (e.g.Donati et al. 2023).
This stellar activity is strongly related to the magnetic field hosted by PMS stars.It is therefore essential to characterize the large-scale magnetic field of such stars to better understand the underlying dynamo processes at work in the stellar interior, but also to analyse star-planet interactions (e.g.Strugarek et al. 2015).In particular, Zeeman-Doppler Imaging (ZDI; Semel 1989;Donati et al. 2006), an efficient tomographic method, already allowed one to reconstruct the large-scale magnetic topology of active PMS stars from time series of circularly polarized spectra (e.g.Donati et al. 2016Donati et al. , 2023;;Yu et al. 2017Yu et al. , 2019;;Finociety et al. 2021Finociety et al. , 2023)).New tomographic methods, such as Time-dependent Imaging of Magnetic Stars (TIMeS; Finociety & Donati 2022), are under development to take into account the temporal evolution of stellar magnetic fields, which would help to unveil the long-term evolution of the field and to identify whether these stars feature a magnetic cycle or not.
In this paper we focus on V1298 Tau, an early-K post T Tauri star with no reported accretion (David et al. 2019a), known to be the youngest solar-mass star with a multiplanetary system for which the radii of Jupiter-like planets have been precisely measured to date.V1298 Tau is located in the Taurus constellation at a distance of 108.0 ± 0.2 pc (Gaia Collaboration et al. 2023).From Gaia measurements, the effective temperature and the logarithmic surface gravity are found to be equal to  eff = 4941 +31 −16 K and log  = 4.227 +0.010 −0.008 (Gaia Collaboration et al. 2023), respectively, in agreement with previous values derived by Suárez Mascareño et al. (2021) and David et al. (2019a), of 5050 ± 100 K / 4.25 ± 0.20 and 4970 ± 120 K / 4.246 ± 0.034, respectively.Suárez Mascareño et al. (2021) also outlined that the effective temperature is consistent with that inferred from the   photometry obtained with the two Micron All-Sky Survey (2MASS) and from the Johnson  photometry collected with the AAVSO Photometric All-Sky Survey (APASS) when assuming an extinction  ( − ) = 0.061, corresponding to an extinction   = 0.19 ± 0.10 following the standard extinction law.A previous independent estimate of the color excess from the multicolor Maidanak observations yields consistent value [ ( −   ) = 0.06 and   = 0.24 mag; Grankin 2013].
Given the stellar rotation period and the line-of-sight projected equatorial velocity we obtained (Sec.3.2), one can infer a relatively accurate estimate of the stellar radius, equal to  * = 1.43 ± 0.03 R ⊙ , assuming the stellar rotation axis is perpendicular to the line of sight.Using the Gaia estimate for the log , one therefore infers a stellar mass  * = 1.26 ± 0.06 M ⊙ .From  eff and  * , we also estimate the luminosity relative to the Sun to be equal to  * / ⊙ = 1.094 ± 0.050.Our updated values of the stellar parameters are summarized in Table 1.These estimates are consistent with those previously determined by Suárez Mascareño et al. (2021), equal to 1.170 ± 0.060 M ⊙ , 1.278 ± 0.070 R ⊙ and 0.954 ± 0.040, for the mass, radius and relative luminosity, respectively, from absolute magnitudes and evolutionary models.From its position in the temperature-luminosity diagram, using our updated stellar parameters, we find that the age of V1298 Tau is ∼10-15 Myr (using the models of Siess et al. 2000 or Baraffe et al. 2015; Fig. 1), making V1298 Tau younger than previously estimated (e.g.David et al. 2019a;Suárez Mascareño et al. 2021).More recent studies focussing on the kinematics of V1298 Tau, based on Gaia DR3, suggest that the star may belong to the D2 (Krolikowski et al. 2021) or D3 (Gaidos et al. 2022) subgroups of the Taurus star forming region, both being younger than 10 Myr.This comes as further argument in favour of V1298 Tau being younger than initially estimated.The evolutionary models of Siess et al. (2000) and Baraffe et al. (2015) also predict that the  (d) 2.91 period used to phase data  rot (d) 2.910 ± 0.005 from  ℓ  rot (d) 2.909 ± 0.009 from RVs  sin  (km s −1 ) 24.9 ± 0.5 from ZDI Ω eq (mrad d −1 ) 2185.1 ± 0.5 from ZDI dΩ (mrad d −1 ) 82.0 ± 2.0 derived from ZDI age (Myr) 10 ± 1 from of Baraffe et al. (2015) age (Myr) 14 ± 1 from of Siess et al. (2000) Figure 1.Location of V1298 Tau in the HR diagram, using the updated  eff and  * from Table 1.The black lines correspond to evolutionary tracks from models of Baraffe et al. (2015) for the stellar mass indicated above the curves while isochrones are depicted as magenta dashed lines.
depth of the convective envelope of V1298 Tau is about 40% of the stellar radius, similar to that of AB Dor ( * = 1 M ⊙ ) or LQ Hya ( * = 0.95 M ⊙ ), two young solar-mass stars slightly older than V1298 Tau (40 -50 Myr; Donati et al. 2003a).V1298 Tau is a target of great interest to study planetary systems at very early stages of their evolution.In the last few years, many studies therefore focused on this system to further characterize the planet properties.In particular, from transit events observed in the K2 and TESS light-curves, the orbital periods of the three innermost planets (b, c and d) were precisely measured (David et al. 2019b;Feinstein et al. 2022).However, and with only 2 transits reported so far, the orbital period of planet e is less well determined, with 17 possible values ranging from 42.7 to 62.8 d (Feinstein et al. 2022), 7 of them being more plausible, between 43.3 and 55.4 d (Sikora et al. 2023).From velocimetric measurements, Suárez Mascareño et al. (2021) estimated the mass of planets b and e to be equal to 0.64 ± 0.19 and 1.16 ± 0.30 M , respectively, while Sikora et al. (2023) found (i) a lower mass for planet e (0.66 ± 0.26 M ), (ii) only an upper limit for planet b's mass of 0.50 M and (iii) a lowsignificance 2.2 candidate detection of planet c with a mass of 0.062 +0.029  −0.028 M .A more recent study of Turrini et al. (2023), based on the same RV dataset as in Suárez Mascareño et al. (2021), reports a slightly larger mass for planet e, equal to 1.23 +0.35  −0.48 M , more consistent with the estimate of Suárez Mascareño et al. (2021).In addition, the obliquity of the orbital axis of planets b and c with respect to the stellar rotation axis is found to be low, suggesting that these planets underwent a smooth migration inwards within the protoplanetary disc (Feinstein et al. 2021;Gaidos et al. 2022).V1298 Tau also offers the opportunity to study atmospheres of young planets, by measuring or constraining the atmospheric evaporation of the innermost planets resulting from the high activity of the host star (Poppenhaeger et al. 2021;Vissapragada et al. 2021;Maggio et al. 2022).
V1298 Tau has been monitored between late 2019 and early 2023 with the near-infrared (NIR) spectropolarimeter SPIRou in the framework of the SPIRou Legacy Survey (SLS), a large programme allocated 310 nights at the Canada-France-Hawaii Telescope (CFHT), within a PI programme, and most recently, within the SPICE large programme (174 nights allocated between late 2022 and mid 2024) at CFHT, aiming at consolidating and enhancing the results provided by the SLS.We start this paper with a detailed description of our spectropolarimetric data set (Sec. 2).We then present the results concerning the large-scale magnetic field of V1298 Tau, reconstructed thanks to ZDI and TIMeS (Sec.3).In Sec. 4, we analyse our velocimetric measurements to constrain the mass of the planets hosted by V1298 Tau and the orbital period of planet e.In Sec. 5, we focus on the stellar activity using three activity proxies in the NIR domain (the He i triplet at 1083 nm and the Paschen  and Brackett  lines).We finally summarize and discuss our results in Sec. 6.

OBSERVATIONS
We extensively observed V1298 Tau with SPIRou between late 2019 and early 2023, collecting 181 high-resolution spectra ranging from 950 to 2500 nm at a spectral resolving power of ∼70,000 (Donati et al. 2020) (Donati et al. 1997).We show the journal of observations in Appendix A.
Data reduction was first carried out with a version of the Libre-ESpRIT pipeline, orginally developped for ESPADOnS data (Donati et al. 1997), adapted for SPIRou observations (Donati et al. 2020).In addition, we corrected the spectra from telluric lines using a PCA approach similar to the one developped in Artigau et al. (2014), providing us 176 usable telluric-corrected spectra in both unpolarized (Stokes ) and circularly polarized (Stokes ) observations, leaving out 5 spectra whose quality is too low for our spectropolarimetric analysis (Table A1).The signal-to-noise ratio (SNR) per pixel in the  band of the usable spectra range from 80 to 253 (median of 212).We then applied Least-Squares Deconvolution (LSD; Donati et al. 1997) on the 176 spectra, using a mask generated with the VALD-3 database (Ryabchikova et al. 2015), featuring an effective temperature and a log  of 5000 K and 4.5, respectively, and containing moderate to strong atomic lines only (relative depth > 10% of the continuum).The resulting Stokes  LSD profiles show clear Zeeman signatures whose typical peak-to-peak amplitude is of 0.2%, with SNRs ranging from 1400 to 5700 (median of 4770).As for the other similar young planet-hosting star AU Mic (Donati et al. 2023), we however do not observe clear distorsions of the Stokes  LSD profiles, suggesting that our spectroscopic data are only weakly affected by brightness features at the surface of the star, in agreement with the small amplitude of the photometric fluctuations in the TESS light curve (2021 Sep 16 -Nov 06) of the order of a few percent1 .We will therefore only focus on the reconstruction of the large-scale magnetic field from our LSD profiles in the following (Sec.3).
We also processed the 181 observations with the nominal SPIRou pipeline, called APERO (Cook et al. 2022), much better optimised for RVs than Libre-ESpRIT.Applying the line-by-line (LBL) method (Artigau et al. 2022) on each reduced spectrum yielded 174 accurate RV values, with a typical error bar of 14 m s −1 (Table A2).Seven spectra were rejected (on 2019 Dec 08, 2019 Dec 13, 2020 Jan 26, 2020 Jan 27, 2020 Jan 28, 2020 Nov 06 and 2022 Feb 01), affected by instrumental effects and/or bad weather conditions, which prevented us to obtain accurate RVs at these dates.

Longitudinal field
We computed the longitudinal field  ℓ , i.e. the line-of-sightprojected component of the vector magnetic field averaged over the visible hemisphere, as the first moment of our 176 Stokes  LSD profiles (Donati et al. 1997).We find that  ℓ ranges from −115 to 89 G between 2019 and 2023, with a typical uncertainty of 9.5 G.As we expect  ℓ to be rotationally modulated, we modelled these values using a quasi-periodic (QP) Gaussian Process (GP; Rasmussen & Williams 2006), following the definition given by Rajpaul et al. (2015): where  and  ′ correspond to the dates of two observations. 1 is the amplitude of the GP,  2 is the exponential decay time-scale (estimating the typical lifetime of the features generating the signal, e.g.spots),  3 is the recurrence period (close to the stellar rotation period  rot ) and  4 is the smoothing parameter setting the amount of short-term variations included in the fit.We also added another parameter, called  5 , representing a potential excess of uncorrelated noise, not taken into account in our measured error bars (e.g.due to  instrinsic variability).We then want to maximize the log likelihood function log L defined in Eq. ( 2) : where  is the vector of the measured  ℓ , whose length is equal to  0 (i.e. to the number of observations).K corresponds to the covariance matrix associated with the quasi-periodic kernel,  represents a diagonal matrix containing the variance of our measurements and S = 2 5 I with I being the identity matrix.We used a MCMC approach, thanks to the emcee python module (Foreman-Mackey et al. 2013), to explore the parameter domains and sample the posterior distribution of our 5 parameters, based on 20000 iterations of 150 walkers.We then estimated the optimal set of parameters (maximizing log L) from these posterior distributions (see Table 2), after removing a burn-in period of 5000 iterations.
We find that the recurrence period is well constrained ( 3 = 2.910± 0.005 d) and consistent with the estimate of the stellar rotation period provided by Suárez Mascareño et al. (2021). ℓ data are fitted down to  2  = 0.61, the RMS of the residuals being smaller than the typical error bar (see Fig. 2), suggesting that these data are not significantly affected by intrinsic variability as demonstrated by a value of  5 smaller than the typical uncertainty and compatible with 0 within ∼2.
As expected by the low value of the decay time-scale, close to one month ( 2 = 36 ± 5 d, i.e. about 4 times smaller than a typical observing season), we observe a clear evolution of the longitudinal field between late 2019 and early 2023, and within each season.In particular, we note that the full amplitude of the curve has a minimum in 2020 (∼30 G) and a maximum in 2021 (∼130 G).This evolution already indicates that we need to split the dataset of each season into several subsets to apply ZDI that assumes a static configuration of the magnetic topology (Sec.3.2).
Running GP regression on individual observing season, we find that the decay time-scale is much less constrained and equal to 68 We note that we were not able to measure the surface field of V1298 Tau, unlike what has been done for the other young star AU Mic, given the large line-of-sight projected equatorial velocity ( sin ) value (∼25 km s −1 , see Sec. 3.2) that prevented us from obtaining reliable measurements from rotationally broadened unpolarized line profiles.

Zeeman-Doppler Imaging
In order to reconstruct the large-scale magnetic field of V1298 Tau, we used ZDI (Semel 1989;Brown et al. 1991;Donati & Brown 1997;Donati et al. 2006), a tomographic method allowing one to invert Stokes  and  LSD profiles into a large-scale magnetic topology at the surface of an active star.As this problem is ill-posed (infinite number of topologies compatible with the data at a given  2  level), ZDI uses the principle of maximum entropy to select the simplest solution, i.e. the one with the minimal amount of information needed to fit the data down to the requested  2  level, by iteratively adding magnetic features at the surface of the star and comparing the associated LSD profiles with the observed ones.
In practice, we divide the stellar surface into a grid of 3000 cells, initially with no magnetic field.At each iteration, the local Stokes  and/or  LSD profiles are computed from each cell using the Unno-Rachkovsky's solution of the polarized radiative transfer equations in a plane-parallel Milne Eddington atmosphere (e.g.Landi Degl'Innocenti & Landolfi 2004), and then integrated over the visible stellar hemisphere, providing the synthetic LSD profiles corresponding to the reconstructed image.As mentioned in Finociety et al. (2023), we replaced the built-in prescription for the limb-darkening law by a linear relation for the continuum only, choosing a limbdarkening coefficient  = 0.3 compatible with  eff = 5000 K and log  = 4.5 in the  band (Claret & Bloemen 2011).
The large-scale magnetic field is described by the sum of a poloidal and a toroidal component, both described as spherical harmonic expansions using the formalism of Donati et al. (2006) 2 .In particular, the poloidal component is characterized by sets of complex coefficients noted  ℓ, and  ℓ, while the toroidal component is fully described by a set of complex coefficients  ℓ, , where ℓ and  denote the degree and order of the spherical harmonic mode in the expansion.We limited our study to the spherical harmonic modes up to a degree ℓ max = 10 as adding higher degree modes does not significantly improve the reconstruction.
ZDI also allows one to reconstruct the distribution of dark and bright features at the surface of active stars from distorsions of the Stokes  LSD profiles.However, we find that, for V1298 Tau, brightness inhomogeneities have only a minor impact on our data (and our reconstructed magnetic maps) as only very weak variations, due to spots and plages, are observed in our Stokes  LSD profiles.
We usually assume that the star rotates a solid body in a first step, then include surface differential rotation (DR) in the reconstruction as described in Sec.3.3.All results shown below were obtained taking into account the mean DR we estimated for V1298 Tau, following the process detailed in Sec.3.3.
Before applying this method to our dataset, we first computed the rotation cycles  of our observations using the following ephemeris: where the initial date (BJD 0 = 2459500) was arbitrarily chosen and  rot = 2.91 d, consistent with our estimates from  ℓ and RVs (see Sec. 4) and that of Suárez Mascareño et al. (2021).Given the short decay time-scale of magnetic features, we were forced to split our dataset into 10 subsets, gathering between 12 and 22 LSD profiles spread over a maximum of 40 d (see Table A1 for the definition of the subsets 3 ).We then applied ZDI on each of these subsets, assuming a line model featuring a mean wavelength, Doppler width and Landé factor of 1750 nm, 3.4 km s −1 and 1.2, respectively.Assuming that the orbital axis of the transiting planets coincides with the stellar rotation axis , we should have  close to 90 • .However, we set  = 80 • in ZDI to reduce the degeneracy between the northern and southern hemispheres in the reconstruction process.Using the magnetic topology of 2020 Aug-Sep (see Fig. 3), we simulated dataset with the same phase coverage, SNRs and  sin  for an inclination of  = 80 • and 90 • .Applying ZDI with  = 80 • for both datasets shows that the impact on the reconstruction of the 10 • difference for  is minimal with respect to the original distribution, and helps to remove part of the north/south degeneracy, therefore justifying our choice.We derive  sin  = 24.9 ± 0.5 km s −1 as part of the imaging process.
For each of the 10 subsets, both Stokes  and  LSD profiles were fitted down to a unit  2  ; we show the reconstructed magnetic topology in Fig. B1 (for 2019) and Fig. 3 (for 2020-2023) while the brightness maps are shown in Fig B2 .As expected from the weak distorsions of the Stokes  LSD profiles, the brightness maps show only a few low contrasted features, associated with a low spot coverage of about 3%, with an episodic spot showing up at some epochs in the the polar region of the northern hemisphere.We note that the SNRs of the Stokes  LSD profiles associated with the 2019 SPIRou data are about twice smaller than those of the other epochs (partly due to a lower exposure time), thereby reducing our ability to clearly detect Zeeman signatures, especially in 2019 Nov-Dec (see Fig. B3).It suggests that the associated reconstructed images are less reliable than those at later epochs.In the following, we will therefore only discuss magnetic reconstructions derived from data collected between 2020 Aug and 2023 Feb.
Despite a rapid evolution of the longitudinal field, we find that the magnetic topology of this star remains mainly poloidal, this component enclosing between 65% (2021 Oct) and 90% (2020 Aug-Sep) of the reconstructed magnetic energy, more than 45% of which being concentrated in axisymmetric modes, until 2023 Jan (see Table 3).This field component mainly consists of a dipole, whose inclination with respect to the rotation axis is low (< 35 • ) until 2023 Jan and whose strength ranges from 85 to 245 G.In 2023 Jan, we see a change in the topology, especially for the radial component of the field, with the dipole component of the poloidal field being now tilted at 75 • , potentially suggesting the beginning of a polarity reversal.To investigate this point further and find out whether this evolution is real or rather caused by, e.g., limited data sampling, we carried out the following experiment.Using the ZDI map of 2022 Nov, we simulated a dataset with the same phase coverage and SNR as that of our 2023 Jan data, and applied ZDI to these simulated data.As the reconstructed large-scale magnetic topology is close to the input field, we can conclude that the strong tilt of the magnetic topology that we report from out 2023 Jan data is likely real.
The large-scale magnetic topology of V1298 Tau also features a significant toroidal field, enclosing between 10% (2020 Aug-Sep) and 35% (2021 Oct) of the reconstructed magnetic energy, and found to be mainly axisymmetric and dipolar.We also find that the quadraticallyaveraged large-scale magnetic flux over the stellar surface, noted <  >, varies by ∼80% between 2020 Aug and 2022 Nov (from 100 to 180 G). 3 The temporal distribution of the observations was too sparse at some epochs, preventing us to use 16 out of the 149 LSD profiles (e.g.isolated observations in 2020 Nov).
Table 3. Properties of the large-scale magnetic field of V1298 Tau between 2020 and 2023, from ZDI reconstruction (including differential rotation) of 8 datasets.Columns 1-2 show the number of the dataset and the associated epoch.Columns 3 to 5 list the quadratically-averaged large-scale magnetic flux over the stellar surface, the strength and tilt, with respect to the stellar rotation axis, of the dipole (of the poloidal component).The two last columns gather the fraction of reconstructed magnetic energy enclosed by the poloidal/toroidal components of the field and in the axisymmetric modes of these components.The typical error bars on these values are of 5-10%.

Differential rotation
From ZDI reconstructions, we see that our data show some variability within each observing season.Part of this variability can potentially be due to surface differential rotation.One can use ZDI to estimate DR by assuming that the surface shear follows a solar-like law given by: where Ω eq and dΩ correspond to the parameters of the DR law, i.e. the rotation rate at the equator and the pole-to-equator rotation rate difference, respectively, and  represents the colatitude.
As our subsets cover up to ∼40 d, we need to merge some of them to diagnose subtle temporal variations of the Stokes  LSD profiles under the effect of DR.We thus estimated DR during the 2020 Aug -2021 Jan season from subsets #3 to #5 (49 Stokes  profiles spanning 135 d, i.e. ∼46 rotation cycles), and during the 2021 Sep -2022 Feb season from subsets #6 to #8 (47 observations spanning 109 d, i.e. ∼37 rotation cycles).The Zeeman signatures collected between 2019 Oct and Dec are both too noisy and too sparse to provide us with reliable estimates of DR, while mediocre phase coverage in 2023 and the sudden evolution of the overall topology between the first half and the second half of the 2022 Nov -2023 Feb season prevents us to estimate DR at this time.
In practice, we fitted these two subsets of Stokes  LSD profiles at a given amount of information (i.e.magnetic energy), over a grid of DR parameters.This process yielded  2  maps from which we estimated the optimal values for both Ω eq and dΩ and their associated error bars by adjusting a 2D paraboloid close to the minimum (Donati et al. 2000;Petit et al. 2002;Donati et al. 2003b;Finociety et al. 2021Finociety et al. , 2023)).
While both estimates of Ω eq are consistent, estimates of dΩ differ by more than 3, suggesting that the DR may vary at the surface of V1298 Tau between the two consecutive observing seasons (see Fig. 4).In addition, for both observing seasons, the confidence contours are slightly distorted, reflecting that the  2  maps are not perfect 2D paraboloids, especially when moving away from the minimum of the maps, which likely increases the size of the error bars.
The higher level of DR measured from 2020 Aug -2021 Jan data is consistent with the evolution of the longitudinal field.We indeed estimate a lower decay time-scale of  ℓ during this observing season, indicating a faster evolution of the large-scale field that may partly reflect a stronger DR.
We chose a unique set of parameters to describe the DR at the surface of V1298 Tau, computed as the weighted means of the estimates provided by both considered observing seasons, i.e.Ω eq = 2185.1 ± 0.5 mrad d −1 and dΩ = 82.0 ± 2.0 mrad d −1 .These values imply that the equator of the star rotates with a period of 2.875 ± 0.001 d while the pole rotates in 2.988 ± 0.003 d, the equator lapping the pole by one cycle in ∼76 d, i.e. about half an observing season.Using these values, we find that both our  ℓ (Sec.3.1) and RV data (Sec.4.2), with a period of 2.910±0.005and 2.909±0.009d, respectively, probe the same latitudes of about 30 • (between 31-37 • and 28-39 • , respectively).
Splitting each data sets in two subsets, each covering about half the observing season, and carrying the same analysis yield discrepant DR parameters, with error bars up to 3 times larger than those mentioned above, demonstrating that a large number of Stokes  profiles collected over the whole observing season are needed to reliably estimate the DR parameters.

Time-dependent Imaging of Magnetic Stars
As our data are spread over several years, we also applied the new tomographic method named Time-dependent Imaging of Magnetic Stars (TIMeS) outlined in Finociety & Donati (2022).As for ZDI, this method aims at finding the simplest large-scale magnetic topology consistent with the data (i.e.Stokes  LSD profiles) but this time allowing it to evolve with time.TIMeS uses sparse approximations to identify as few spherical harmonic modes (associated with  ℓ, ,  ℓ, and  ℓ, coefficients) as possible to reconstruct the magnetic topology and GPs to model the time-dependence of each identified coefficients.This method can only be applied in the cases where the magnetic field is not too strong, typically with a magnetic flux lower than ∼3 kG for a  sin  = 25 km s −1 .Given the  sin  of V1298 Tau and the values of <  > and of the local field inferred from the ZDI reconstructions in the previous sections (≪ 3 kG), we therefore assume that TIMeS can be applied to our data.
In practice, thanks to the principle of sparsity, we first look for  maps computed from Stokes  LSD profiles over a grid of parameters Ω eq , the rotation rate at the equator, and dΩ the difference in rotation rate between the pole and the equator.Red contours show the 3 confidence interval, while the magenta crosses depict the optimal value with its associated error bars, derived from the fit of a 2D paraboloid close to the minimum of the grid.The right panel shows the contours associated with 1 and 3 confidence intervals obtained from both observing seasons (2020 Aug -2021 Jan: blue; 2021 Sep -2022 Feb: red).The magenta cross depicts the mean parameters.
the simplest linear combination of known magnetic topologies that provides a satisfactory fit to subsets of  consecutive profiles, assuming that the large-scale magnetic field does not significantly evolve over the time span covered by the  profiles.We did not consider the subsets for which 2 consecutive profiles are separated by more than twice the decay time-scale of the longitudinal field.As mentioned in Finociety & Donati (2022), the  profiles should sample the rotation cycle while typically covering 10 to 20% of the decay time-scale of the longitudinal field, explaining why we set  = 6 in this paper (corresponding to ∼20% of the  ℓ decay time-scale).This first step yields the spherical harmonic modes that significantly contribute to the overall topology over the whole dataset, with a value of the corresponding  ℓ, ,  ℓ, and  ℓ, coefficients for each subset.TIMeS then models these time series of coefficients using GPs with a squared exponential kernel, allowing one to compute the magnetic maps and associated Stokes  LSD profiles for each observed rotation cycle.
We applied TIMeS to all the Stokes  LSD profiles covering 2020 Aug -2023 Feb (i.e.subsets #3 to #10) at once, allowing all spherical harmonic modes up to ℓ max = 10 (as in ZDI, see Sec. 3.2) and assuming  = 70 • to better identify axisymmetric structures4 (e.g.dipole) that are missed when using  = 80 • .We were able to fit the data down to  2  = 1.49, revealing that the model does not succeed in fitting the data down to the noise level, likely due to the limited available resolution at the surface of the star, the model only featuring 11 modes up to ℓ = 4, even though ℓ max = 10.We show the reconstructed magnetic maps during the 2020 Aug -2021 Jan season at specific rotation cycles in Figs. 5, C6.Despite some loss of spatial resolution, we find that the reconstructed maps are quite similar to those obtained with ZDI at the same epochs, with only a subtle evolution within each subset considered in the previous sections (as expected from the unit  2  obtained with ZDI).In particular, we identify the same topology for the radial field and the similar large features in the azimuthal field (e.g.around phases 0.35-0.45 in late 2020).
This first application of TIMeS to real data therefore demonstrates that this method is efficient at reconstructing reliable magnetic topologies even if the data span several years, with only a small amount of information that is not reconstructed (especially the smallest features).In addition, we find that the fraction of the magnetic energy enclosed in the poloidal field is globally consistent with that derived from ZDI applied on each subset, except for subsets #6 and #10 which likely results from the loss of spatial resolution with TIMeS (Fig. 8).
Using Lomb-Scargle periodograms (Zechmeister & Kürster 2009) with the astropy python module5 on each of the  ℓ, ,  ℓ, and  ℓ, coefficients found with TIMeS does not reveal any common periodicity (see Fig. C7), suggesting that we did not yet detect a magnetic cycle in our data, which is not surprising given that only a hint of a first polarity switch has been observed up to now.Future observations of V1298 Tau will allow us to conclude whether the putative change in the poloidal field polarity is confirmed, and is part of a cycle much longer than the time span of our current observations.

Orbital period of planet e
We analysed the RVs obtained from the data reduction performed with APERO using the LBL method (Artigau et al. 2022) in order to look for the RV signatures induced by the 4 transiting planets hosted by V1298 Tau.The orbital period and transit times of these planets were accurately determined by Feinstein et al. (2022), except for planet e for which the orbital period is still unknown (but longer than 40 d; Suárez Mascareño et al. 2021;Feinstein et al. 2022;Sikora et al. 2023).Using these values (see Table 4) and assuming circular orbits for all planets (consistent with the results of Suárez Mascareño et al. 2021 andSikora et al. 2023), we only need to retrieve the semiamplitude of the planet RV signatures (  ,   ,   and   for planet b, c, d and e, respectively) as well as the orbital period of planet e, noted   , from our data.
In practice, we used a model featuring a QP GP to model the RV signal induced by stellar activity, 4 Keplerian curves associated with the planet signatures (assuming circular orbits) and an additional white noise (e.g.due to intrinsic variability).As for the longitudinal field, we sampled the posterior distribution of each parameter (planet and GP) using a MCMC method to estimate their optimal value.
We note that the posterior distribution of   is multimodal if we impose a uniform prior on the orbital period of planet e (between 42 and 65 d, following the list of probable periods of Feinstein et al. 2022), showing 1 global maximum at about 53.5 d and three local maxima at 42.7 d, 46.4 d and 60.6 d (see Fig. D1).We therefore run our process again, each time with a narrow Gaussian prior centred on one of these periods with a standard deviation of 1 d to investigate which one is the most likely.We then computed the marginal logarithmic likelihood log L M following the method described by Chib & Jeliazkov (2001), in order to estimate the significance of one solution with respect to the others from the difference in log L M (i.e.logarithmic Bayes Factor / log BF).Considering the model featuring an orbital period for planet e at 53.5 d as the reference, we find that the log BF for the 3 other models (associated with an orbital period of 42.7, 46.4 and 60.6 d, respectively) is equal to −1.4, −2.4 and −4.1, suggesting that the model featuring a prior centred on an orbital period of 53.5 d is more likely than the three others.From the MCMC approach, we therefore estimate the most likely orbital period to be   = 53.5 ± 0.4 d and we thus use this model to further characterize the 4 planets in the following.

Planet masses
With our reference model, we are able to fit our data down to  2  = 5.97 with the residuals exhibiting a RMS dispersion of 34.8 m s −1 .This value for the  2  suggest that discrepancies exist between the model and our measurements and that our error bars (estimated from photon noise) are likely underestimated, e.g.due to intrinsic variability that has not be taken into account.More generally, our detections of the planet-induced RV variations are limited by the high level of stellar activity compared to the recovered planet RV signatures (GP amplitude of ∼120 m s −1 , i.e. 2.6 to 28 times larger than the semi-amplitude of the planet signatures) and of intrinsic variability (modeled by the additional white noise, equal to 46 ± 5 m s −1 , i.e. 3× larger than the typical error bar of 14 m s −1 ) of V1298 Tau.In particular, we find that the decay time-scale is equal to 32 +6 −5 d, compatible with the one derived from  ℓ measurements, indicating that the surface of the star evolves rapidly.The GP is also modulated by the stellar rotation with a period equal to 2.909 ± 0.009 d, similar to the one found from  ℓ measurements.All these results clearly illustrate the need for robust filtering methods given how difficult the detection of planets around young active stars is, even when their ephemerides (orbital periods, transit times) are well determined.
Despite this intense activity, we find that planet e is best detected, at a 3.9 level, with a RV semi-amplitude   = 47 +16 −12 m s −1 , associated with a mass   = 1.01 +0.35  −0.26 M (for an orbital period of 53.5 ± 0.4 d).Assuming a radius of 0.89 ± 0.04 R , following Feinstein et al. (2022), the planet density is found to be   = 1.79 +0.66  −0.51 g cm −3 .For the three innermost planets, we do not obtain clear detections as we find   = 4.2 ± 10.7 m s −1 ,   = 4.8 ± 6.2 m s −1 and   = 4.9 ± 6.7 m s −1 , corresponding to masses of   = 0.07 ± 0.18 M ,   = 0.05 ± 0.07 M and   = 0.05 ± 0.08 M , and densities of   = 0.14 ± 0.36 g cm −3 ,   = 0.76 ± 0.99 g cm −3 and   = 0.40 ± 0.66, respectively.Based on the 99% confidence level interval associated with the posterior distribution of each planet mass derived with the MCMC approach (after converting the sampled semi-amplitudes into masses), we derive upper limits on these masses of   < 0.44 M ,   < 0.22 M and   < 0.25 M .We summarize all the planet parameters in Table 4 for the 4 best periods and we show the associated corner plots in Appendix D. We also show the best fit (i.e. for   = 53.5 d) and the associated phasefolded filtered RVs in Figs. 9 and 10, respectively.
Assuming an eccentric orbit for planet e6 (in the case of the most likely orbital period of 53.5 d), yields no improvement in confidence   In black, we show the fraction of the reconstructed magnetic energy enclosed in the poloidal field, estimated using TIMeS.The red points correspond to the value reconstructed from ZDI applied on each subsets, assuming typical error bars of 10% (see Table 3).The vertical blue dotted lines depict the limits between each subsets used in ZDI (with their number at the top of the plot) while the magenta dashed line denotes the limit above which the magnetic field is mainly poloidal (i.e. more than 50% of the reconstructed magnetic energy is enclosed in the poloidal component).level (log BF = −0.7 with respect to our best model) and an eccentricity consistent with 0 (0.03 ± 0.05), indicating a posteriori that our previous assumption of circular orbits was justified to model our data.
In addition, we computed the Lomb-Scargle periodograms of the raw RVs, filtered RVs and residuals RVs to look for potential additional planets that may have not been detected before (e.g.nontransiting planets).Looking at the raw RVs, we do not see any significant peak at 1 yr (with a FAP lower than 10%), suggesting that the contamination of our data from telluric lines is insignificant for our purpose.When focussing on the filtered RVs (in which only activity was removed), we see a significant power at the orbital period of planet e, but no clear peak in the periodogram of the residual RVs (see Fig. 11).We simulated RV datasets with the same temporal distribution as our actual SPIRou data using our best model on which we added the signature of a potential fifth planet beyond the orbit of V1298 Tau e.Using these datasets, we estimate that only planets more massive than 2.7 and 5.9 M at 0.7 and 1.4 au (corresponding to an orbital period of 6 months and 1.5 yr), respectively, could be reliably detected with our dataset.

Additional constraints from previous studies
From the transit times of planet e in the TESS and K2 light curves, Feinstein et al. (2022) derived 17 probable orbital periods for this planet.We therefore modeled our data, assuming circular orbits, this time imposing a Gaussian prior on the orbital period of this planet centered on each of these values with a standard deviation of 0.0001, corresponding to the typical error bar reported by Feinstein et al. (2022).
From the log L M , we find that the most likely period is   = 53.0039± 0.0001 d, yielding a fit to our RV data associated with a  2  = 6.30 and a RMS dispersion of the residuals of 35.8 m s −1 .For this period, we report a 4 detection of the RV signature of planet e (  = 44 +15 −11 m s −1 ) associated with a mass of   = 0.95 +0.33 −0.24 , while the 3 other planets remain undetected.
Using this value as a reference, we computed the log BF between this model and the 16 others.From this criterion, we cannot definitely rule out any of the 17 periods provided by Feinstein et al. (2022) as all periods yield a | log BF| < 6.4.We however distinguish 3 groups: the most likely periods (| log BF| < 2.5), the plausible periods (2.5 < | log BF| < 5) and the least likely periods (| log BF| ≥ 5).In addition to   = 53.0039d, 3 other values belong to the most likely group (  = 54.2085,46.7681 and 45.8687 d, associated with log BF of −1.1, −2.0 and −2.4 respectively).The intermediate plausible group contains 6 orbital periods (  = 45.0033, 47.7035, 48.6770, 59.6294, 61.1583 and 62.7678 d), that are still considered as potential values, even though the derived mass for planet e can be more than twice smaller than that obtained in Sec.4.2 or in the most likely group.The least likely group therefore gathers the 7 remaining values (  = 44, 1699, 49.6911, 50.7484, 51.8516, 55.4692, 56.7899 and 58.1750 d).We summarized the results for the 4 most likely orbital periods in Table 5, while those associated with the plausible and least likely groups are given in Tables D1 and D2, respectively.
As for the previous case in Sec.4.2, estimating the eccentricity of planet e as part of the process does not signifcantly improve the results (e.g.log BF = 0.2).

CHROMOSPHERIC ACTIVITY
We focussed on three lines known to probe stellar activity in the NIR domain, namely the He i triplet at 1083 nm, and the hydrogen lines Paschen beta (Pa) at 1282 nm and Brackett gamma (Br) at 2165 nm (Zirin 1982;Short & Doyle 1998).We show the individual spectra and the median one in Fig. 12.We note that the extreme blue wing of Pa does not reach the continuum because of the nearby Ca i line that blends with it.However, this feature does not vary more than the continuum and should thereferore have no significant impact on our analysis.
We quantified the variations of the lines with respect to the median computing the equivalent width variations (EWVs) as done for several other young stars (e.g.Finociety et al. 2021Finociety et al. , 2023;;Donati et al. 2023).In practice, we divided each telluric-corrected spectra by the median one and we computed the equivalent width of the residual spectra by integrating between −40 and +40 km s −1 in the stellar rest frame.The EWVs show a dispersion of 2.70, 1.26 and 1.15 km s −1 , for the He i, Pa and Br lines, respectively.We find that these indices are only weakly correlated with the stellar rotation.Fitting the EWVs of each line using a QP GP, we note that only the He i EWVs are rotationnaly modulated with a period  = 2.907 ± 0.018 d, close to the values found with our  ℓ and RV measurements.For this indicator, the excess of uncorrelated noise, tracing the intrinsic variability, is about 22× larger than the formal photon noise error bar for each indicator (median of 0.10 km s −1 ), reflecting the intense activity triggered by V1298 Tau, as already suggested by our RV fit.In particular, we note that the amplitude of this weak modulation (i.e. of the GP) is equal to 2.23 km s −1 i.e. about 65% that of the excess of white noise.

SUMMARY AND DISCUSSION
We presented the analysis of NIR spectropolarimetric observations of V1298 Tau collected with SPIRou between 2019 Oct 02 and 2023 Feb 06, aimed at characterizing the large-scale magnetic topology of Table 4. MCMC results for the analysis of the RV data. 1 st and 2 nd columns list the parameters and their associated priors corresponding to a model featuring an activity signal and four planet-related RV signatures.The 4 last columns show the best value for each parameter, in the case of a narrow prior centred on 3 different values of   0 for the orbital period of planet e.The knee of the modified Jeffreys priors is set to the typical RV uncertainty, noted   .Note that the log BF is computed with respect to the marginal likelihood associated with   0 = 53.5 d.

Parameter
Prior the host star and at further constraining the masses and densities of the 4 transiting planets.

Evolution of the magnetic field of V1298 Tau
Our Stokes  LSD profiles show clear Zeeman signatures demonstrating that V1298 Tau hosts a strong large-scale magnetic field.The longitudinal field  ℓ , tracing this large-scale component, is modulated by the stellar rotation with a period of 2.910±0.005d, consistent with, though much more accurate than the estimate of 2.91 ± 0.05 d inferred from RV measurements by Suárez Mascareño et al. (2021). ℓ evolves rapidly, over a time-scale of 36 ± 5 d (the full amplitude of the pattern varying from ∼30 G in 2020 to ∼130 G in 2021), suggesting that the large-scale magnetic topology changes both within an observing season and between two consecutive seasons.
Given the fast evolution of  ℓ , we split our data into several subsets per season on which ZDI was applied independently to reconstruct the large-scale magnetic topology of V1298 Tau.The large-scale magnetic topology of V1298 Tau is found to be mainly poloidal and axisymmetric with a toroidal component becoming less axisymmetric from 2020 to 2023.Comparing the magnetic reconstructions of this star with those of two young solar-mass stars similar to V1298 Tau, namely AB Dor (e.g.Donati et al. 2003b) and LQ Hya (e.g.Lehtinen et al. 2022), we find that V1298 Tau hosts a simpler large-scale magnetic field, with a mostly axisymmetric radial field mainly characterized by the dipole component of the poloidal field unlike both other stars.All three stars also continuously show azimuthal features at their surface, that most likely indicate that the underlying dynamo processes operate throughout the entire convective zone of these stars, as suggested by Donati et al. (2003b).In addition, ZDI reconstructions from our latest observing epoch show a rapid change in the radial field, with the dipole component (of the poloidal field) becoming much more tilted in early 2023 (from 20-30 • to 75 • ), and possibly tracing a polarity reversal, as that recently reported for LQ Hya (Lehtinen et al. 2022).These changes suggest that the mag-Table 5. MCMC results for the analysis of the RV data for the 4 most likely planet e orbital periods derived by Feinstein et al. (2022). 1 st and 2 nd columns list the parameters and their associated priors corresponding to a model featuring an activity signal and four planet-related RV signatures.The 4 last columns show the best value for each parameter, in the case of a narrow prior centred on different values of   0 for the orbital period of planet e.The knee of the modified Jeffreys priors is set to the typical RV uncertainty, noted   .Note that the log BF is computed with respect to the marginal logarithmic likelihood associated with   0 = 53.0039d.

Parameter
Prior netic field of V1298 Tau follows a long-term evolution potentially similar to the 11-year solar cycle.Part of the observed magnetic variability is likely attributable to DR at the surface of the star, found to be equal to 82.0±2.0 mrad d −1 , i.e. ∼1.5 times that of the Sun.This value implies a typical lap time of ∼76 d for the equator to lap the pole by one rotation, implying that DR clearly distorts the regular modulation of spectropolarimetric data over a typical observing season spanning about 4.5 months (i.e.∼135 d, see Table A1).This level of DR is similar to that of young solar-mass stars whose stellar structure resembles that of V1298 Tau, such as AB Dor (dΩ ABDor of the order of 70 mrad d −1 ; Jeffers et al. 2007) and LQ Hya (dΩ LQHya of the order of 100 mrad d −1 ; Donati et al. 2003b).In addition, our results suggest that the DR at the surface of V1298 Tau may be varying with time, passing from 90.6 ± 2.8 mrad d −1 in 2020 Aug -2021 Jan to 73.5 ± 2.8 mrad d −1 in 2021 Sep -2022 Feb.This variation of ∼15 mrad d −1 is larger than the one observed for the Sun (< 1 mrad d −1 ) but is in good agreement with the observations of AB Dor and LQ Hya for which large variations of DR have also been reported from year to year, with, in particular, dΩ ABDor ranging from 60 to 75 mrad d −1 (Jeffers et al. 2007) and dΩ LQHya ranging from 14 to 200 mrad d −1 (Donati et al. 2003b;Kővári et al. 2004;Lehtinen et al. 2022).The global increase in the strength of the magnetic field of V1298 Tau between 2020 Aug -2021 Jan and 2021 Sep -2022 Feb (see the evolution of <  > in Sec.3.2) may help to slow down the DR and may therefore contribute to the observed decrease in dΩ between both seasons.The DR parameters of the Sun vary periodically with a period corresponding to that of the solar cycle and this peridocity is correlated with solar activity (Poljančić Beljan et al. 2022).However, with only two measurements of DR at the surface of V1298 Tau so far and no reported magnetic cycle, we cannot conclude whether DR parameters for V1298 Tau follow the same trend as those of the Sun.
We used, for the first time on actual data, the new tomographic method TIMeS (Finociety & Donati 2022) allowing one to recon- struct the long-term temporal evolution of magnetic fields without splitting the whole data set.We fitted all the Stokes  LSD profiles down to  2  = 1.49, suggesting that discrepancies still exist between the synthetic model and the observed data.This is not surprising as Finociety & Donati (2022) mentioned in the original paper that the method is likely to miss the smallest features that only contribute little to the observed data.One major difference between TIMeS and ZDI is that the model derived with TIMeS only uses some spherical harmonics modes that strongly contribute to the magnetic topology over the whole time span of the observations to describe the magnetic field while the model derived with ZDI uses all spherical harmonic modes up to a maximum degree ℓ max .We therefore suppose that the TIMeS method failed at fitting the data down to a unit  2  because some high degree modes do not have a strong contribution to the magnetic field over all epochs (e.g.small features that may have a small lifetime) and they were therefore not selected for the model (only modes up to ℓ = 4 were indeed identified by the method).Future improvements of the method to increase the flexibility in the choice of the modes should help to reduce these discrepancies (e.g. through new criteria of selection, implementation of DR in the model), but the results we obtained with this first application on real data are already quite promising.We indeed find an evolution of the magnetic topology with TIMeS similar to that reconstructed from the application of ZDI on data subsets.In addition, TIMeS allows one to access to the temporal evolution of the large-scale field and its decomposition into poloidal and toroidal components.Comparing the fraction of reconstructed magnetic energy enclosed in the poloidal field derived from ZDI and TIMeS, we see that we have a global consistency, except at some epochs at the beginning or end of a season (typically 2021 Oct, 2023 Jan).To investigate whether we can distinguish a magnetic cycle within our data set, we computed periodograms of the coefficients reconstructed with TIMeS.We however do not see any significant period shared by all coefficients, and therefore conclude that, if this star has a magnetic cycle, it must be longer than the time span of our 2020-2023 observations (897 d, i.e. ∼2.5 years), as already suggested by the fact that only the beginning of a potential reversal has been identified in our ZDI reconstruction.Thorough spectropolarimetric monitorings of V1298 Tau are thus clearly needed in the coming years, to help us identify whether the dynamo processes at work in this star are similar to those operating in the Sun (e.g.Ω processes) or more complex ones (e.g. 2 ,  2 Ω processes, see Rincon 2019 for a review).

Constraining the planet parameters
Using the LBL method (Artigau et al. 2022), we computed 174 accurate RV values, over the 3.5 years of SPIRou monitoring, from which we estimated the mass and orbital period of the outermost planet (e) and upper limits on the mass of the three innermost planets.Assuming circular orbits, our SPIRou RV measurements yield 4 probable orbital periods for planet e at 42.7 ± 0.3 d, 46.4 ± 0.4 d, 53.5 ± 0.4 d and 60.6 ± 0.7 d.Computing the log L M associated with each model, we find that, among these solutions, the most likely model is the one featuring   = 53.5 ± 0.4 d.We note that our best value disagrees with those found found by Sikora et al. (2023) and Turrini et al. (2023), equal to 46.768131±000076 d and 45.46 +1.10  −0.23 d, respectively, but these values remain compatible with one of our less likely solutions (46.4 ± 0.3 d) within 1 and 2.4, respectively.Although the log BF suggests that the period   = 42.7 ± 0.3 d is the second most likely, we consider it to be less probable as this value corresponds to the lower limit provided by the TESS monitoring (Feinstein et al. 2022).We also consider that our largest period (60.6±Using our most likely model, i.e. the one featuring a period of 53.5 d, we obtain the most significant detection of the RV signature of planet e since the new constraints on its orbital period provided by the TESS monitoring.We indeed report a 3.9 detection of the RV signature of planet e, yielding a mass of 1.01 +0.35  −0.26 M .This value is intermediate between those obtained by Sikora et al. (2023), Suárez Mascareño et al. (2021) and Turrini et al. (2023), equal to 0.66 ± 0.26, 1.16 ± 0.30 M and 1.23 +0.35  −0.48 M , respectively, but still compatible within 1 with all values.From our mass value and the updated radius of planet e provided by Feinstein et al. (2022), we estimate the density of this planet to be equal to 1.79 +0.66 −0.51 g cm −3 , slightly lower than those given in Feinstein et al. (2022) and Turrini et al. (2023), using similar planet radius, but still compatible within 1.As already noted by previous studies, this density is suprisingly high for a young planet, the implications of this value being discussed in the next Section.
The RV signatures of the three innermost planets are undetected from our SPIRou measurements.We therefore estimate the following upper limits on these masses, based on the 99% confidence interval of their posterior distributions sampled through the MCMC approach:   < 0.44 M ,   < 0.22 M and   < 0.25 M .These results illustrate the need for (i) powerful filtering methods to get rid of activity jitter in RV measurements to clearly unveil the planet prop-erties and (ii) more observations to better sample the RV signature induced by each planet of this system.
Our results are generally consistent with those provided by Sikora et al. ( 2023) but less so with those of Suárez Mascareño et al. (2021).We outline that the masses derived by Suárez Mascareño et al. (2021) may be unreliable as a result of a potential overfit of their RV data as recently suggested by Blunt et al. (2023).Although our results do not allow us to strongly constrain the mass of the three innermost planets, we note that our estimate of   = 4.2 ± 10.7 m s −1 disagrees at a ∼3 level with the value of 41 ± 12 m s −1 previously reported in Suárez Mascareño et al. (2021), suggesting that planet b may be several times less massive.This result is further strengthened by the upper limit of 0.08 M provided by the HST transmission spectra of V1298 Tau b (Barat et al., submitted, 2023), fully consistent with our estimate of   = 0.07 ± 0.17 M .
Our estimate of planet c's mass (  = 0.05±0.07M ) is similar to the estimate of 19.8 +9.3 −8.9  ⊕ (i.e.0.062 +0.029 −0.028 M ) derived by Sikora et al. ( 2023), but with error bars 2.3× larger, both results being fully consistent within 1.Our RV data does not allow us to bring stronger constraints on the mass of planet d with respect to the previous studies, our upper limit being ∼1.6× larger than that of Sikora et al. (2023) but still ∼1.6×lower than that of Suárez Mascareño et al. (2021).The derived upper mass limits yield upper limits for the planet densities of   < 0.89 g cm −3 ,   < 2.94 g cm −3 and   < 1.89 g cm −3 .Future observations of V1298 Tau will help to refine these densities and suggest what their internal structures could be.To compare the results from both previous velocimetric studies and ours, we show the position of the four planets in a mass-radius diagram in Fig. 13.
We further confirm that the assumption of circular orbits is best adapted to model our RV data, in particular for planet e for which we derive an eccentricity consistent with 0 (with an error bar of 0.05).
Using the transit times of planet e derived from the K2 and TESS light curves Feinstein et al. (2022) listed 17 probable values for   with a precision of 0.0001 d.We therefore used these values as a prior on   for the modeling of our RV data.As outlined in Sec.4.3, the models associated with these periods can be sorted according to the value of their log L M .As expected from our main analysis, we find that, among these 17 periods,   = 53.0039± 0.0001 d is the most likely, this value being the closest to our best period and consistent within 1.3.Three periods yield a log L M slightly smaller than that of the model featuring   = 53.0039d (i.e.|Δ log L M | < 2.5).Among these 3 periods, one remains close to 53.5 d (i.e.  = 54.2085d) while the 2 others (  = 46.7681and 45.8687 d) are in good agreement with one of the other solutions of our main analysis, i.e.   = 46.4± 0.4.These results thus strengthen our initial conclusion of an orbital period close to 53 d for the outermost planet of the V1298 Tau system.Future photometric transits of planet e are needed to further pinpoint whether its actual orbital period is ∼46 or ∼53 d.
With an orbital period of   = 53.0039d, the RV semi-amplitude of planet e is found to be equal to   = 44 +15 −11 m s −1 , corresponding to a mass and density of 0.95 +0.33  −0.24 M and   = 1.66 +0.61 −0.48 g cm −3 , respectively.These values are therefore fully consistent with and only sligthly lower than those derived in our main analysis, as expected from the small difference in the orbital period of planet e between both analyses.

Formation and evolution of the V1298 Tau system
According to the orbital periods of the four planets, the V1298 Tau system is likely not in a resonant chain with only planet b, c and d having respective period ratios close to 3:2 (for c vs d) and 2:1 (for d vs b).A resonant chain would indeed favour an orbital period of planet e close to 48 or 49 d rather than 53 d (or even 46 d), which seems unlikely given our data (see Tables D1 and D2).Our finding is consistent with previous studies focussing on the formation of this planetary system (Tejada Arevalo et al. 2022;Turrini et al. 2023).This implies that either the V1298 Tau system did not form in a resonant chain, or if it did, it must have undergone a rapid phase of dynamical instabilities that broke the resonant chain after the disc dispersal (i.e. in 10-15 Myr given the young age of V1298 Tau).
Although theoretical models indicate that disc migration in a multiplanetary system tends to favour the planets to be in a resonant chain, several processes, such as planet-disc interactions, may prevent the formation of such chains (e.g.interaction between a planet and the wake of another one; Baruteau & Papaloizou 2013).We also outline that the turbulence of the circumstellar disc can play a role on the disc-planet interactions that eventually disrupt the strict mean motion resonance between the planets (Pierens et al. 2011;Paardekooper et al. 2013).Turrini et al. (2023) discussed several mechanisms that may break a resonant chain after its formation.The first one is the eccentricity damping through tidal dissipation (Batygin & Morbidelli 2013) but this scenario seems unlikely as the timescale on which it would occur for planet c is 2 orders of magnitude larger than the age of the system (Turrini et al. 2023).Photoevaporation of planets c and d may also contribute to the break of resonance (Goldberg et al. 2022) but, once again, such a phenomenon would occur over more than 100 Myr, inconsistent with the age of the star (Turrini et al. 2023).More plausible explanations involve planetesimal scattering or planet-planet scattering.However, planetesimal scattering seems unlikely for this system as simulations show that a planetesimal belt of 50  ⊕ would not have a significant impact on the V1298 Tau system and a much more massive belt would be required to break the resonance (Turrini et al. 2023).Turrini et al. (2021) suggest that the most likely scenario to explain the non-resonance between planets b and e relates to planet-planet scattering, involving the presence of one or more planets on an eccentric orbit beyond the orbit of planet e.We however consider that this scenario is quite unlikely given that the eccentricity of planet e we derived from our data is consistent with 0. In addition, such additional planets would be difficult to detect through photometric transits and RV monitoring due to their long orbital periods.In particular, no hint of a fifth planet in this system has been reported so far and, using GAIA and SPHERE observations, Turrini et al. (2023) indicate that no planet more massive than 2 M should orbit beyond 50 au from the host star (i.e. with an orbital period longer than 300 yr).From simulations, we also estimated that our data excludes the presence of planets more massive than 2.7 M and 5.9 M orbiting at 0.7 and 1.4 au (i.e.corresponding to orbital periods of 6 months or 1.5 yr), respectively, in good agreement with the results of Turrini et al. (2023) for similar orbital periods.
According to our best model, planet e has a suprisingly high density given its young age 7 .In particular, the mass and density we derived for this planet are consistent with the plausible composition inferred by Sikora et al. (2023) from models of Fortney et al. (2007) accounting for the irradiation from a Sun-like host star at a distance of 0.1 au 8 , i.e. a 100  ⊕ core composed of a 50/50 by mass ice-rock mixture surrounded by a H/He envelope.In addition, even though our derived mass is slightly lower than the one inferred from their recent study, Turrini et al. (2023) suggest that the ice-rock core may be more massive than 100  ⊕ with about 2/3 of planet e composed of heavy elements and 1/3 of hydrogen and helium.This would mean a relative core mass ∼6× larger than that of Jupiter according to formation models (Venturini & Helled 2020), i.e. a planet composition closer to that of Neptune than that of Jupiter.This result is in good agreement with the models of Baraffe et al. (2008), indicating that planet e should be composed of more than 50% of heavy elements (up to 90%) to match the measured radius (0.89 R ) at this young an age (10-15 Myr).However, these models predict that, assuming a metallicity of 0.10 (as in the host star; Suárez Mascareño et al. 2021), the planet should be older than 10 Gyr given the measured radius, which is clearly unrealistic.
To date, the cause of such a high density of planet e is still speculative.One explanation relates to the critical mass for the onset of runaway gas accretion onto a protoplanetary core.In the case of planet e, this mass would be larger than that generally admitted (10 7 Assuming solar metallicity and no irradiation effects, the models of Baraffe et al. (2008) predict that a 1 M planet of 10-15 Myr orbiting around a solarmass star should have a radius of ∼1.3 R , implying a mean density of about 0.56 g cm −3 (i.e.3.2× smaller than the one we derived for V1298 Tau e).Accounting for the irradation from a Sun at 0.045 au, the predicted density is even smaller, close to 0.45 g cm −3 . 8The semi-major axis of the planets hosted by V1298 Tau are equal to 0.172 ± 0.003, 0.084 ± 0.001, 0.110 ± 0.002 and 0.293 ± 0.005 au, for planet b, c, d and e, respectively, assuming the orbital periods of Feinstein et al. (2022) for the 3 inner planets and an orbital period of 53.5 ± 0.4 d for planet e.Assuming an orbital period of 53.0039 ± 0.0001 d for planet e yields a consistent semi-major axis of 0.291 ± 0.005 au.
to 20  ⊕ ) possibly due to uncertainties on the hypotheses on which models rely such as the size, mass and composition of the dust, having a direct impact on the opacity and on the cooling of the gas that can be accreted by the planet.Alternative explanations for the high density of V1298 Tau e may be related to a high accretion rate of pebbles, enriching the content of heavy elements of the planet, through different mechanisms.One generally assumes that a planet on a circular orbit reaches its pebble isolation mass (PIM, i.e. the mass above which the planet cannot accrete more pebbles) for  ∼10 − 20  ⊕ depending on the disc's aspect ratio around the planet (e.g., Ataiee et al. 2018).However, in some conditions, a planet on an eccentric orbit can reach a much higher PIM (> 100  ⊕ for moderate eccentricity, e.g. > 0.04, and a turbulent visosity  = 0.005; Chametla et al. 2022), which could be the case of V1298 Tau e, according to the estimate we derived (0.03 ± 0.05).The accretion of dust and gas through a vertical (or meridional) circulation (Binkert et al. 2021;Szulágyi et al. 2022) is another process that can partly contribute to the increase in the heavy elements content in V1298 Tau e even if the planet has reached its PIM, that deserves to be further investigated.Turrini et al. (2023), focussing on the evolution of the V1298 Tau system, also discussed other possible explanations for the high density observed for planet e.One is that this planet formed beyond the CO 2 snowline and then underwent a large-scale migration within a massive circumstellar disc in such a way to accrete heavy elements along its pathway (as probably for HAT-P-20 b; Thorngren et al. 2016), e.g. through the accretion of an intense planetesimal flux preventing runaway accretion of the disc gas to occur (Turrini et al. 2023).As also noted by Turrini et al. (2023), another option is that this density results from an accretion of vapor coming from an enhancement of the volatile content of the circumstellar disc due to drifting and evaporating pebbles when they reach their ice lines (e.g.Booth & Ilee 2019;Schneider & Bitsch 2021).Characterizing the atmosphere of planet e through transit observations would help to discriminate between both hypotheses as the accretion of planetesimal is expected to yield a lower refractory-to-volatile ratio as well as a different C/O ratio than the second process (Turrini et al. 2021;Schneider & Bitsch 2021;Pacetti et al. 2022).We however caution that these two processes may not be consistent with what we know of V1298 Tau.In particular, the metallicity of the host star is close to solar (0.10 ± 0.15; Suárez Mascareño et al. 2021), and does not hint at a high content of heavy elements in the circumstellar disc that would lead towards the second hypothesis.Though much older (∼1.5 Gyr; Ment et al. 2018), HAT-P-2 b ( = 8.62 +0.39  −0.55 M ,  = 0.951 +0.039 −0.053 R ; Loeillet et al. 2008) is another example of a hot Jupiter with an even higher density ( = 12.5 +2.6 −3.6 g cm −3 ; Loeillet et al. 2008) that may be related to a large content of heavy elements (up to 500  ⊕ ; Baraffe et al. 2008).The most likely scenario to explain the formation of HAT-P-2 b is that the planet is actually the result of several collisions with planetary embryos (Baraffe et al. 2008).In this scenario, the absolute mass of heavy elements in the planet increases and the gas is ejected from the planet, therefore increasing the relative mass of heavy elements with respect to the initial mass of the planet.This mechanism may be an alternative explanation for the formation of V1298 Tau e, which deserves to be also investigated in future studies.Maggio et al. (2022) mentioned that the two outermost planets (e and b) should not be affected by a strong photoevaporation given their large mass with respect to their radius.However, this conclusion was drawn from the mass of planet b inferred by Suárez Mascareño et al. (2021), which is larger than the upper limit that we report in the current paper and that of Sikora et al. (2023), suggesting that the conclusion may be revised taking into account our updated mass upper limit.In addition, Turrini et al. (2023) indicate that, during its formation, planet b may have undergone a migration with heavy element enrichment (as for planet e).However, this conclusion was also drawn from a density derived assuming the mass obtained by Suárez Mascareño et al. (2021) and thus likely overestimated and inconsistent with our data.In particular, Turrini et al. (2023) suggest that planet b may contain up to 90  ⊕ (i.e.∼0.28 M ) of heavy elements, representing 64% of our mass upper limit (0.44 M ), which thus seems unlikely given our upper limit on the density of planet b (< 0.89 g cm −3 ).The scenario for the formation of planet b can therefore be revised when the actual mass of this planet is further constrained.
According to Maggio et al. (2022), planets c and d may undergo strong atmospheric evaporation if their masses is lower than ∼40  ⊕ (∼0.13 M ).Our mass estimates are thus in favour of strong atmospheric photoevaporation, strengthening the conclusion of Sikora et al. (2023) about planet c.The V1298 Tau system, hosting 4 planets younger than 30 Myr, is an ideal laboratory to study star-planet interactions at early stages of the evolution.Future monitoring of this target will therefore help us to better understand the impact of such interactions on the formation and evolution of young planets.
In a forthcoming work, we plan to couple all existing RV measurements (i.e.collected with HARPS-N, MAROON-X, SPIRou and HIRES; Suárez Mascareño et al. 2021;Sikora et al. 2023;Blunt et al. 2023) to perform a consistent analysis and hopefully to refine the planet masses.Such a study should include two potentially coupled GPs to simultaneously and consistently model the activity signals in the optical and NIR domains.In parallel, we will pursue our spectropolarimetric monitoring of V1298 Tau thanks to the SPICE Large Programme in the next semesters.In light of the results reported in the present paper, these future observations will help us to fulfill two major goals.The first one is to reconstruct the large-scale magnetic topology of the star using ZDI and TIMeS from new data to better characterize the dynamo process at the origin of the observed field and its evolution (e.g.potential magnetic cycle subject to Ω processes).The second goal is obviously to refine estimates of the mass of the four planets, but also to firmly ascertain the orbital period of planet e, a quite challenging task given the amplitude of the stellar activity jitter more than 2.6× larger than the semi-amplitude of the presumably most massive planet (i.e.V1298 Tau e).Obtaining accurate measurements of the mass of the four planets of this system is an essential step to further constrain the mass-radius relation at early stages of evolution, as only few other planets aged of less than 30 Myr have precise measured radius and mass (AU Mic b and c, HD 114082 b and HIP 67522 b; Donati et al. 2023;Zakhozhay et al. 2022;Rizzuto et al. 2020), but also to better understand the evolution of young planets and the associated star-planet interaction.Coupling spectroscopic and high-precision photometric observations (e.g.TESS, PLATO), especially during transit events, will clearly help us to refine the orbital parameters, the atmospheric escape, if any, (e.g.study of the He i line during transit events) and the mass of the 4 planets.More generally, we would need multi-instruments campaigns (e.g.combination of optical/NIR spectropolarimeters, such as ESPaDOnS, SPIRou, CRIRES+) would provide huge benefits to better characterize the magnetic topology of V1298 Tau and the properties of its planetary system.2023).We also plot 3 isodensities as thin grey lines, corresponding to a density of 0.50 (solid), 1.00 (dashed) and 2.00 (dotted) g cm −3 .
MNRAS 000, 1-24 (2023)         MNRAS 000, 1-24 (2023) Table D1.MCMC results for the analysis of the RV data for the 6 planet e orbital periods derived by Feinstein et al. (2022) identified as plausible from our data.1 st and 2 nd columns list the parameters and their associated priors corresponding to a model featuring an activity signal and four planet-related RV signatures.The 6 last columns show the best value for each parameter, in the case of a narrow prior centred on different values of   0 for the orbital period of planet e.The knee of the modified Jeffreys priors is set to the typical RV uncertainty, noted   .The log BF is computed with respect to the marginal logarithmic likelihood associated with   0 = 53.0039d (see Table 5).
. Altogether, 97 of these observations were collected in the framework of the SLS (2019 Oct 02 -2021 Jan 08, 2022 Feb 01), more specifically within the WP3 package focusing on the magnetic topology of PMS stars, 52 within the PI Programme of Benjamin Finociety (run ID 21BF21, 2021 Sep 17 -2022 Jan 30) and the 32 remaining ones within the SPICE large programme (2022 Nov 02 -2023 Feb 09).More specifically, the observations are spread over four observing seasons: (i) ∼140 d from 2019 Oct to 2020 Feb, (ii) ∼135 d from 2020 Aug to 2021 Jan, (iii) ∼137 d from 2021 Sep to 2022 Feb; and (iv) ∼99 d from 2022 Nov to 2023 Feb.Each observation is composed of a sequence of 4 subexposures of ∼500 s taken at different azimuths of the polarimeter retarder in order to get rid of potential spurious signals in the polarization and systematic errors at first order

Figure 2 .
Figure 2. Longitudinal field measurements of V1298 Tau obtained from SPIRou data over our observing period (top panel) with a zoom on the 2020 Aug -2021 Jan (middle panel) and 2021 Sep -2022 Feb (bottom panel) seasons.For each panel, the top plot shows he raw measurements as black dots with their associated error bars while the red curve represents the quasi-periodic GP fit with the 1 confidence area in light orange.The bottom plots show the residuals between the raw measurements and the fit, showing a RMS dispersion of 8.7 G, i.e. smallest than the typical error bar of 9.5 G, and a  2  = 0.61.
Aug -2021 Jan, 2021 Sep -2022 Feb and 2022 Nov -2023   Feb, respectively, reflecting a slower evolution of the magnetic field during the 2021-2022 observing season.

Figure 3 .
Figure 3. ZDI maps of the radial (left), meridional (middle) and azimuthal (right) magnetic field components for the 8 subsets associated with the following epochs: 2020 Aug-Sep, 2020 Sep-Oct 2020, 2020 Dec -2021 Jan, 2021 Oct, 2021 Nov-Dec, 2022 Jan, 2022 Nov and 2023 Jan.The star is shown in a flattened polar view down to −60 • , with the pole at the center, the equator depicted as a bold circle and the −30 • , 30 • and 60 • latitude parallels represented by dashed circles.Positive radial, meridional and azimuthal fields are shown in red and point outwards, polewards and counter-clockwise, respectively.The ticks around the star refer to the phases of the SPIRou spectropolarimetric observations used in ZDI.These maps were obtained assuming the mean differential rotation derived in Sec.3.3 and are shown at mid-time throughout the corresponding subset.MNRAS 000, 1-24 (2023)

Figure 4 .
Figure 4. Surface differential rotation of V1298 Tau as measured from SPIRou data during the 2020 Aug -2021 Jan (upper left) and 2021 Sep -2022 Feb (bottom left) seasons. 2 maps computed from Stokes  LSD profiles over a grid of parameters Ω eq , the rotation rate at the equator, and dΩ the difference in rotation rate between the pole and the equator.Red contours show the 3 confidence interval, while the magenta crosses depict the optimal value with its associated error bars, derived from the fit of a 2D paraboloid close to the minimum of the grid.The right panel shows the contours associated with 1 and 3 confidence intervals obtained from both observing seasons (2020 Aug -2021 Jan: blue; 2021 Sep -2022 Feb: red).The magenta cross depicts the mean parameters.
6 and 7 while those associated with 2021 Sep -2022 Feb and 2022 Nov -2023 Feb data can be seen in Appendix C (Figs. C1 to C5).The synthetic Stokes  LSD profiles obtained with TIMeS are also shown in Fig.

Figure 5 .
Figure 5. Magnetic maps in 2020 Aug-Sep reconstructed with ZDI without differential rotation (first row) and TIMeS (second to fifth rows).As TIMeS allows the reconstruction of a map at each observed date, we only show four specific epochs indicated on the left.See Fig. 3 for a detailed description of the Figure.

Figure 8 .
Figure 8. Evolution with time of the poloidal component of the field between 2020 and 2023.In black, we show the fraction of the reconstructed magnetic energy enclosed in the poloidal field, estimated using TIMeS.The red points correspond to the value reconstructed from ZDI applied on each subsets, assuming typical error bars of 10% (see Table3).The vertical blue dotted lines depict the limits between each subsets used in ZDI (with their number at the top of the plot) while the magenta dashed line denotes the limit above which the magnetic field is mainly poloidal (i.e. more than 50% of the reconstructed magnetic energy is enclosed in the poloidal component).

Figure 9 .
Figure 9. RVs of V1298 Tau between end 2019 and early 2023.The first panel shows the raw RVs (black dots) along with the best model (including the activity jitter and the planet sginatures) in red and the 1 confidence level in light orange.2 nd to 5 th panels show the retrieved RV signatures of planets b, c, d and e, respectively.The last panel shows the residuals RVs, having a RMS dispersion of 24.3 m s −1 , i.e. close to the typical error bar of our measurements.

Figure 10 .
Figure 10.Phase-folded filtered (top panels) and residuals (bottom panels) RVs for planets b (1 st panel), c (2 nd panel), d (3 rd panel) and e (4 th panel).The 174 measured RVs are shown as black dots (with their error bar) while the red stars correspond to RVs averaged over 0.1 orbital cycle.The RV signature (derived from our MCMC process) associated with each planet is shown as a blue solid line.The RMS dispersion of the residuals is equal to 34.8 m s −1 (see Fig. 9).

Figure 11 .
Figure 11.Lomb-Scargle periodograms of the raw RVs (top), filtered RVs (middle) and residuals (bottom) computed using the astropy python module.The horizontal dotted, dashed and solid red lines depict the FAPs of 0.1, 1 and 10%, respectively.The vertical dashed magenta line corresponds to the stellar rotation period while the dashed cyan lines show the orbital periods of planets e, b, c and d (from left to right).The orange curve depicts the window function, the orange solid and dotted vertical lines correspond to the synodic period of the Moon and its 1-yr aliases, respectively, and the orange dashed vertical line outlines the 1-yr period.

Figure 12 .
Figure 12.Spectra of V1298 Tau for the He i triplet (top), Pa (middle) and Br (bottom).For each panel, the top curve shows the median profile with the dispersion in each velocity bin (in light grey) while the bottom curves depict the superposition of all individual telluric-corrected spectra.The red vertical dashed lines show the location of the lines.The blue wing of Pa is likely blended with a Ca i line.

Figure 13 .
Figure 13.Mass-radius diagram.Top panel: Posterior distributions of the planet masses derived from the MCMC approach.V1298 Tau b, c, d and e are shown in magenta, orange, green and blue, respectively.Bottom panel: The light grey dots correspond to confirmed exoplanets with radius and masses known with a precision of 30% or lower.The full squares show the position of the 4 planets according to the mass and error bars derived in the present paper (assuming the radii of Feinstein et al. 2022), while the left-pointing triangles depict the upper limit associated with the 99% confidence level.The open up-pointing triangles show the estimates and 3 upper limits from Suárez Mascareño et al. (2021) while the open stars correspond to recent estimates and 2 upper limits from Sikora et al. (2023).Upper limits are identified by an arrow pointing to the left.The same color code as for the upper panel is used to identify the 4 planets.For comparison, the black dots correspond to the position of the two young planets AU Mic b and c from the estimates of Donati et al. (2023).We also plot 3 isodensities as thin grey lines, corresponding to a density of 0.50 (solid), 1.00 (dashed) and 2.00 (dotted) g cm −3 .MNRAS 000, 1-24 (2023)

Figure B1 .Figure B2 .
Figure B1.Maps reconstructed with ZDI from SPIRou data collected in 2019 Oct (top) and 2019 Nov-Dec (bottom).See Fig. 3 for a detailed description.These maps are less reliable than those of Fig. 3 due to a SNR of the Stokes  LSD profiles about twice smaller than in the other observing epochs.

Figure C7 .
Figure C7.Reconstructed coefficients with TIMeS using Stokes  LSD profiles between 2020 and 2023.Left column: temporal evolution of the 10 coefficients describing the reconstructed magnetic topology, with  ℓ, and  ℓ, characterizing the poloidal component and  ℓ, representing the toroidal component (ℓ and  being the degree and order of the associated spherical harmonic mode in the expansion).Right column: Periodograms associated with the time series of the reconstructed coefficients, computed using the astropy python.The dotted, dashed and solid red lines show the FAP levels associated with a probability of 1, 0.1 and 0.01%, respectively.The orange curve depicts the window function, the orange solid and dotted vertical lines correspond to the synodic period of the Moon and its 1-yr aliases, respectively, and the orange dashed vertical line outlines the 1-yr period.

Figure D1 .
Figure D1.Corner plot of the posterior distributions of the hyperparameters of quasi-periodic GP and the planet parameters obtained through a MCMC approach assuming a uniform prior between 42 and 65 d for   .The red lines show the medians of the posterior distributions, while the orange dashed lines correspond to the value that maximizes the log likelihood function, chosen as the optimal parameter.This plot was created using the corner python module (Foreman-Mackey 2016).

Figure D2 .Figure D4 .
Figure D2.Same as Fig. D1 with a narrow prior centered on 42.7 d for   .

Table 1 .
Stellar properties used in our study.

Table 2 .
Optimal value obtained for the GP modeling of the longitudinal field  ℓ from a MCMC approach.The knee of the modified Jeffreys priors is set to the typical uncertainty on the the longitudinal field, noted   ℓ .

Table A1 .
Spectropolarimetric observations collected between 2019 and 2023, reduced with the Libre-ESpRIT pipeline adapted for SPIRou data.For each observation, we list the date (day, hour), the barycentric Julian date (BJD), the rotation cycle (as defined in Sec.3.2), the SNR in the  band, the noise level in the LSD profile, the longitudinal field  ℓ and the number of the subset used in ZDI (see Sec. 3.2).

Table A2 .
Artigau et al. (2022)bservations collected between 2019 and 2023, reduced with the APERO pipeline.For each observation, we list the date, the barycentric Julian date (BJD) and the associated RV computed from the LBL approach ofArtigau et al. (2022).