The classical T Tauri star CI Tau observed with SPIRou: magnetospheric accretion and planetary formation

We report new observations of the classical T Tauri star CI Tau with the SPIRou near-infrared spectropolarimeter and velocimeter at the Canada-France-Hawaii Telescope (CFHT) in late 2019, 2020 and 2022, complemented with observations obtained with the ESPaDOnS optical spectropolarimeter at CFHT in late 2020. From our SPIRou and ESPaDOnS spectra, to which we applied Least-Squares Deconvolution, we infer longitudinal ﬁelds clearly modulated with the 9-d rotation period of CI Tau. Using Zeeman-Doppler imaging, we reconstruct the large-scale magnetic topology, ﬁrst from SPIRou data only in all three seasons, then from our 2020 SPIRou and ESPaDOnS data simultaneously. We ﬁnd that CI Tau hosts a mainly axisymmetric poloidal ﬁeld, with a 1 kG dipole slightly tilted to the rotation axis and dark spots close to the pole that coincide with the footpoints of accretion funnels linking the star to the inner disc. Our results also suggest that CI Tau accretes mass from the disc in a stable fashion. We further ﬁnd that radial velocities (RV) derived from atomic and CO lines in SPIRou spectra are both rotationally modulated, but with a much lower amplitude than that expected from the putative candidate planet CI Tau b. We conﬁrm the presence of a RV signal at a period of 23.86 d reported in a separate analysis, but detect it clearly in CO lines only and not in atomic lines, suggesting that it likely traces a non-axisymmetric structure in the inner disc of CI Tau rather than a massive close-in planet.


INTRODUCTION
In the last three decades, thousands of extra-solar planets and exoplanet systems have been detected and characterized by various methods, in particular velocimetry, measuring planet masses through the reflex motion planets induce on their host stars, and photometry, measuring planet radii from the flux reduction planets generate as they transit in front of their host stars.However, very few observations exist about nascent exoplanets, younger than, e.g., 20 Myr, with only a couple of multi-planet systems known to date (e.g., David et al. 2019;Plavchan et al. 2020;Donati et al. 2023;Finociety et al. 2023b).Yet, documenting such objects is key to provide constraints on theoretical models of how planets build up and migrate in the protoplanetary ★ E-mail: jean-francois.donati@irap.omp.euaccretion disc surrounding forming stars, and from which both stars and planets feed.
Pre-main-sequence (PMS) low-mass stars, also called T Tauri stars (TTSs), are particularly interesting in this respect, both those that are still surrounded by an accretion disc (the classical TTSs / cTTSs) and the ones that exhausted their discs (the weak-line TTSs / wTTSs).Whereas cTTSs exhibit extended discs with a wide variety of spatial structures, sometimes even featuring distant planets detected through imaging techniques (Keppler et al. 2018), wTTSs are known to host close-in massive planets, some of them even transiting in front of their host stars (David et al. 2016(David et al. , 2019)), giving the opportunity to derive estimates of the planet bulk densities and thereby hints on the planet internal structures.
Magnetic fields are known to play a key role throughout the whole process of star and planet formation (e.g., Hennebelle et al. 2020), including in the cTTS phase where the dynamo-generated large-scale field of the host star is strong enough to evacuate the central regions of the accretion disc and carve a large magnetospheric gap extending up to 10 stellar radii.The magnetic field of the central star connects to that at the inner edge of the accretion disc, controlling accretion between the disc and the star and potentially slowing down the rotation of the host star (e.g., Romanova et al. 2004;Zanni & Ferreira 2013;Blinova et al. 2016;Pantolmos et al. 2020).This star-disc interaction may also impact the fate of newborn planets as they migrate inwards within the accretion disc and pile up near its inner edge (e.g., Mulders et al. 2015), allowing them to potentially escape from being swallowed by the host star (Lin et al. 1996;Romanova & Lovelace 2006).
One such cTTS of particular interest is CI Tau, aged 2 Myr and known to host a massive accretion disc extending to several hundreds of au from the central star (e.g., Guilloteau et al. 2014), inclined at 50 • to the line of sight and featuring several density gaps suggesting ongoing planet formation (Clarke et al. 2018).More recently, interferometric K-band observations of the central regions of the accretion disc showed that CI Tau also hosts an inner disc ring, whose inner edge is located at 0.2 au from the star, and that is misaligned with the main disc and tilted at 70 • to the line of sight (Gravity Collaboration et al. 2023), presumably as a result of stardisc interactions.CI Tau is indeed known to be strongly magnetic, with a large-scale field in excess of 1 kG (Donati et al. 2020a) and a small-scale field of about 2 kG (Sokal et al. 2020).A candidate eccentric close-in massive planet of minimum mass 8.1 M at an orbital period of 9 d was reported, on the basis of CI Tau undergoing regular radial velocity (RV) variations with a semi-amplitude of about 1 km s −1 (Johns- Krull et al. 2016), making this young star even more interesting to study.However, spectropolarimetric observations of CI Tau carried out with ESPaDOnS in the optical domain at the 3.6-m Canada-France-Hawaii Telescope (CFHT) unambiguously demonstrated that the 9-d period is actually the rotation period of the central star, and that the reported RV variations are in fact mostly due to stellar activity, and in particular to the presence of cool surface features coinciding with the footpoints of accretion funnels linking the star to the inner accretion disc, and carried in and out of the observer's view by rotation (Donati et al. 2020a).
To document this interesting case further and better characterize the immediate environment of CI Tau, we carried out new observations, both with the near-infrared (nIR) SPIRou spectropolarimeter / velocimeter (Donati et al. 2020b) and with its optical companion ES-PaDOnS (Donati 2003) at CFHT.In this paper, we start by presenting the SPIRou and ESPaDOnS observations we collected (Sec.2), and briefly recall the main parameters of CI Tau (Sec.3).We then analyse these new data, first measuring the longitudinal field ℓ (Sec.4) then carrying out a full tomographic imaging of the surface of CI Tau (Sec.5) in order to further constrain how magnetospheric accretion proceeds between the inner disc and the surface of CI Tau (Sec.6).We complete this study with a detailed RV analysis of whether CI Tau may host close-in massive planets (Sec.7) before summarizing our results and discussing their implications for our understanding of star / planet formation (Sec.8).

SPIROU AND ESPADONS OBSERVATIONS
CI Tau was regularly observed from early 2019 to early 2023 with the SPIRou nIR spectropolarimeter / high-precision velocimeter (Donati et al. 2020b) at CFHT, first within the SPIRou Legacy Survey (SLS) then within the SLS consolidation and enhancement pro-gramme SPICE.Both the SLS and SPICE are large observational programmes (allocated 310 and 174 CFHT nights over 7 and 4 semesters respectively) dedicated to the study of planetary systems around nearby M dwarfs, and of star / planet formation in the presence of magnetic fields.Additional observations of CI Tau were collected in late 2020 with the ESPaDOnS optical spectropolarimeter (Donati 2003) at CFHT, within the PI program of Jerome Bouvier (run ID 20BF01).SPIRou and ESPaDOnS collect unpolarized and polarized stellar spectra, covering (in a single exposure) wavelength intervals of 0.95-2.50m at a resolving power of 70 000 for SPIRou, and of 0.37-1.00m at a resolving power of 65 000 for ESPaDOnS.Each polarization sequence from both instruments consists of 4 subexposures, associated with different orientations of the Fresnel rhomb retarders (to remove systematics in polarization spectra to first order, see Donati et al. 1997), yielding one unpolarized (Stokes ), one circularly polarized (Stokes ) spectrum, and a null polarization check called allowing one to diagnose potential instrumental or data reduction problems.
A total of 101 SPIRou polarization sequences of CI Tau were collected in 3 main observing seasons, 34 from 2019 October to 2020 February, 39 from 2020 September to 2021 January and 28 from 2022 November to 2023 January, following an unsuccessful early start in 2019 February where 6 polarized sequences were recorded but discarded due to cold weather conditions in which the preliminary Fresnel rhombs used at this time did not behave properly retardationwise 1 .Poor weather conditions also led us to discard two polarized sequences in the first season, and another two in the last season, whereas an instrumental issue with the Fabry-Perot reference channel (strongly outshining the science channels) spoiled 3 polarization sequences in 2020 January.This left us altogether with 94 usable Stokes , and spectra of CI Tau, with 29, 39 and 26 for the first, second and last season respectively, recorded over a total time span of almost 1200 d or 133 rotation cycles.An additional 13 optical spectra of CI Tau were collected with ESPaDOnS in 2020 November and December over a time frame of 16 d or 1.8 rotation cycles.The full observation log for these spectra is given in Appendix A.
Both SPIRou and ESPaDOnS data were processed with Libre ESpRIT, the nominal reduction pipeline of ESPaDOnS at CFHT, that was also adapted for SPIRou (Donati et al. 2020b).Least-Squares Deconvolution (LSD, Donati et al. 1997) was then applied to all reduced spectra, using a line mask constructed from the VALD-3 database (Ryabchikova et al. 2015) for an effective temperature eff =4250 K and a logarithmic surface gravity log =4.0 adapted to CI Tau (see Sec 3).In these masks, we selected atomic lines of relative depth larger than 10 percent for SPIRou spectra (for a total of ≃1500 lines, featuring an average wavelength and Landé factor of 1750 nm and 1.2 respectively), and of relative depth larger than 40 percent for ESPaDOnS spectra (for a total of ≃9000 lines of average wavelength and Landé factor 640 nm and 1.2).The line thresholds were set to ensure that the noise levels in the resulting Stokes LSD profiles are similar in both sets of spectra, ranging from 1.9 to 6.2 (median 2.9, in units of 10 −4 where denotes the continuum intensity) for SPIRou spectra, and 2.6 to 4.0 (median 2.9) for ESPaDOnS spectra (in the same units).In the case of SPIRou spectra, we also used a second mask restricted to CO lines between 2.2 and 2.4 m only and including in particular all CO lines of the 2.3 m CO bandhead (seen in absorption and of mostly stellar origin), which we use mainly in Sec. 7.
All phases and rotation cycles associated with our SPIRou and ESPaDOnS spectra were computed assuming a rotation period of rot = 9.01 d (see Sec. 4) and counting from an arbitrary starting BJD0 of 2458762.0 (taken from our first usable SPIRou observation, see Table A1).

FUNDAMENTAL PARAMETERS OF CI TAU
We briefly recall here the main characteristics of CI Tau, following our previous study based on older ESPaDOnS data (Donati et al. 2020a), and complement it with new results from the refereed literature.older CI Tau is a 2-Myr-old cTTS located at a distance of 160.3 ± 0.4 pc from the Sun (Gaia Collaboration et al. 2023), in the Taurus star forming region.CI Tau has a mass of 0.90 ± 0.02 M ⊙ (Simon et al. 2019), a photospheric temperature of eff = 4200 ± 50 K, a radius of 2.0 ± 0.3 R ⊙ , and rotates with a period of rot = 9.010 ± 0.023 d (see Sec. 4).This value of rot places the corotation radius cor , i.e., the distance from the centre of the star where the Keplerian orbital period in the equatorial plane equals rot , to 0.082 ± 0.001 au or equivalently 8.8 ± 1.3 ★ .With a logarithmic luminosity relative to the Sun of log( ★ /L ⊙ ) = 0.1 ± 0.1, CI Tau is still in the process of completing its phase of nearly isothermal contraction along its Hayashi track, and is presumably still fully convective (Donati et al. 2020a).
Besides, CI Tau is surrounded by a massive accretion disc inclined to the line-of-sight at an angle of 50 ± 4 • (Guilloteau et al. 2014) and extending up to a distance of several hundreds of au.It features a succession of dusty rings and gaps located at radii of ≃13, 39 and 100 au (Clarke et al. 2018) that suggests ongoing planet formation.It was recently discovered that the inner regions of the accretion disc are significantly tilted with respect to the outer regions, with the inner disc being tilted at 70 ± 1 • with respect to the line of sight (Gravity Collaboration et al. 2023).The origin of the inner-outer disc misalignment is not clear, but could be due to magnetic warping caused by star-disc magnetic interactions (e.g., Romanova et al. 2021), to exotic phenomena occurring within the disc (e.g., Demidova & Grinin 2023), and / or to the presence of a massive close-in planet embedded in the disc.In fact, such a planet was invoked for CI Tau (Johns-Krull et al. 2016) but what was assumed to be the planet orbital period finally turns out to be the rotation period of the star, whereas the reported RV semi-amplitude (of 1.08 ± 0.25 km s −1 ) thought to probe the reflex motion of the star under the gravitational pull of the planet is in fact attributable to stellar activity (Donati et al. 2020a).
We summarize in Table 1 the main parameters of CI Tau used in (or derived from) our study.In particular, we assume in the following that the rotation axis of CI Tau is inclined at = 70 • to the line of sight, consistent with the adopted sin , rot , eff and luminosity, and in line with the reported tilt angle of the inner disc.

THE LONGITUDINAL FIELD OF CI TAU
We pursued our analysis by computing the longitudinal field ℓ (i.e., the line-of-sight-projected component of the vector magnetic field averaged over the visible hemisphere) of CI Tau for each of our observing epochs, from the Stokes and LSD profiles derived in Sec. 2, and following Donati et al. (1997).The magnetic field of CI Tau being rather strong and the nIR line profiles being significantly broadened by Zeeman broadening (Sokal et al. 2020), Stokes LSD signatures were integrated over a window of ±30 km s −1 (centred on the stellar rest frame) to estimate the first moment and its error bar, whereas the equivalent width of the Stokes LSD profiles is measured through a Gaussian fit.This window is large enough to ensure that no magnetic information is lost but small enough to minimize the resulting noise on ℓ , with small changes on the window size having minimal impact on the results.As usual (e.g., Donati et al. 1997Donati et al. , 2023)), we also derive a pseudo longitudinal field from the polarization check in the exact same way we do from , which allows us to check that it is consistent with 0 within the error bars, i.e., that it yields a reduced chi-square 2 r close to 1.The ℓ values corresponding to our 94 validated SPIRou visits are listed in Table A1, and range from −170 to 120 G (median −5 G) with error bars of 10 to 30 G (median 15 G).The corresponding 2 r (with respect to ℓ = 0 G line) is equal to 22, demonstrating that the large-scale field of CI Tau is very clearly detected in the recorded Stokes Zeeman signatures; the same operation applied to yields 2 r = 0.92, indicating that no spurious contamination is observed and that our formal error bars, dominated by photon noise, are consistent with the observed dispersion.We note that our ℓ values are typically an order of magnitude weaker than the small-scale fields estimated from the Zeeman broadening of nIR lines (≃2 kG, Sokal et al. 2020), as usual for such active stars whose tangled small-scale fields generate circular polarization signatures that mostly cancel out and thereby do not contribute much to longitudinal fields (e.g., Morin et al. 2010;Donati et al. 2023).
To investigate the quasi-periodic (QP) behaviour of ℓ and its evolution with time, we carry out a Gaussian-Process Regression (GPR) fit of the ℓ curve, with a covariance function ( , ′ ) of the following type: where 1 is the amplitude (in G) of the Gaussian Process (GP), 2 its recurrence period (i.e., rot , in d), 3 the evolution timescale (in d) on which the shape of the ℓ curve changes, and 4 a smoothing parameter describing the amount of harmonic complexity needed to describe the data.We end up selecting the QP GPR fit to our ℓ points (arranged in a vector denoted ) that features the highest likelihood L, defined by: (2) where is the covariance matrix for our 94 epochs, Σ the diagonal variance matrix associated with and =2 5 the contribution from an additional white noise that we introduce as a fifth hyperparameter 5 (in case our error bars on ℓ were underestimated for some reason).The hyper-parameter domain is then explored using a Monte-Carlo Markov Chain (MCMC) process, yielding posterior distributions and error bars for all hyper-parameters.
The results of the GPR fit are shown in Fig. 1 (with a zoom on seasons 2020 and 2022 in the middle and lower panels), whereas the derived hyper-parameters are listed in Table 2.The first conclusion is that rot = 9.010 ± 0.023 d, thereby confirming unambiguously the findings of our previous study (Donati et al. 2020a).We stress in particular that our value of rot is the only one on which the MCMC systematically converged whatever prior we used and initial value we provided for this parameter.The measured period is in fact so well defined that the 1-yr aliases do not even show up in the associated corner plot.Besides, this means that the 6.6 d period quoted in several studies (e.g., Johns-Krull et al. 2016;Biddle et al. 2018Biddle et al. , 2021) can be safely rejected as a potential rotation period of CI Tau.
We also find that CI Tau shows a ℓ curve with a rather simple, nearly sinusoidal and slowly-evolving modulation pattern, hence the relatively high value of the smoothing parameter ( 4 = 0.60 ± 0.11) and of the evolution timescale ( 3 = 146 +42 −32 d) compared to the more evolved young active star AU Mic that exhibits a more complex and rapidly evolving ℓ curve (Donati et al. 2023).We finally outline that our ℓ values can be fitted within the noise level ( 2 r = 0.82) and thus that the reconstructed additional white noise is consistent with zero.
We proceeded in the same way to derive ℓ values from our ES-PaDOnS LSD profiles of CI Tau from late 2020, now integrating over a smaller window of width ±20 km s −1 (in the stellar rest frame) given that Stokes and LSD profiles derived from the optical domain are significantly narrower than those inferred at nIR wavelengths.The 13 ℓ values we obtained from photospheric lines range from −5 to 130 G (median 55 G) with error bars of 12 to 22 G (median 15 G), and are shown in Fig. 2 (left panel).We can already note that the corresponding ℓ curve is significantly different from that obtained in our previous study (with data from early 2017, Donati et al. 2020a), which showed two clear minima over a rotation period whereas that from our 2020 data only shows one.
We also derived ℓ values from the narrow emission component (NC) of the Ca infrared triplet (IRT) lines (see Fig. 2, right panel) and of the He D 3 line (not shown), both known to probe the postshock region at the footpoints of accretion funnels (in addition to the quieter, non-accreting chromosphere for the Ca IRT lines), integrating over a wide window (±30 and ±50 km s −1 for the Ca IRT and He D 3 line respectively) to take into account the large intrinsic widths of these lines 2 .
A sine fit (plus first harmonic) to the ℓ points from ESPaDOnS LSD profiles of photospheric lines yields a period of 9.0 ± 0.3 d (see Fig. 2, bottom left panel), fully consistent (though expectedly less accurate) than that derived from our SPIRou data (spanning a much longer time frame of 1200 d instead of only 16 d).Doing the same on the ℓ data derived from the NC of the Ca IRT lines is less successful, the amplitude of the modulation pattern being less clear as a result of intrinsic variability affecting more strongly accretion lines than photospheric lines (see Fig. 2, bottom right panel); we nonetheless see a modulated trend, with ℓ being stronger and more negative in the second half of the rotation cycle.ℓ values derived from the He D 3 line show no clear modulation (within error bars of ≃300 G), and stay more or less constant at −1 kG.This is again different from the results of our previous study where the corresponding ℓ curve displayed a very obvious modulation, varying from 0 to -2 kG in a more or less sinusoidal fashion (Donati et al. 2020a).
One can note (see, e.g., Fig. 2 and Table A2) that ℓ values from ESPaDOnS LSD profiles are discrepant with those from SPIRou LSD profiles collected a few weeks before and after the ESPaDOnS observations, the SPIRou ℓ curve reaching down to −130 G whereas the ESPaDOnS ℓ points are almost all positive.Similarly, ℓ values from the Ca IRT and He D 3 lines are constantly negative, unambiguously indicating that the accreting regions on the visible hemisphere at the surface of the star are associated with magnetic patches of negative polarity.Despite this discrepancy, we can see that the epoch at which ℓ reaches its minimum value when measured from ESPaDOnS LSD profiles (−3.6 ± 16.1 G, on 2020 November 20 at BJD ≃ 2459175.94,see Table A2) coincides well with that at which ℓ is expected to reach its strongest negative value at SPIRou wavelengths (see Fig. 1 middle panel); this is also the time at which ℓ from Ca IRT lines is most negative (see Fig. 2, bottom right panel).Altogether, it implies that all three ℓ curves are in phase though different in amplitude, shape and vertical offset.
These apparent inconsistencies can in fact be exploited to our advantage in order to more finely characterize what happens at the surface of CI Tau.We propose they directly reflect (i) that photospheric brightness is non-uniform at the surface of CI Tau, (ii) that accreting regions of negative polarity (bright at chromospheric level) tend to coincide with cool regions at photospheric level, that are darker in the optical than in the nIR (as a result of the wavelength dependence of the spot-to-photosphere contrast), and (iii) these dark regions face the observer around phase ≃0.8 in our ephemeris (see Sec. 2) and are located at high latitudes on the visible hemisphere given the low amplitude modulation (<1 km s −1 ) of the RV curves from both emission line NCs (see, e.g., Fig. 2 top right panel for the Ca IRT).This preliminary conclusion is qualitatively consistent with the findings of our previous work on CI Tau (Donati et al. 2020a) as well as of similar spectropolarimetric studies of cTTSs (e.g., Donati et al. 2011).Our combined SPIRou and ESPaDOnS data set of CI Tau from 2020 September to 2021 January is thus quite unique in this respect, yielding an unprecedented richness of complementary diagnostics on the large-scale magnetic field, on the photospheric brightness and on the ongoing embedded accretion flows.We elaborate along these lines in a more quantitative way in Sec.5.3, by simultaneously modelling the Stokes and profiles of our three sets of lines (SPIRou LSD profiles, ESPaDOnS LSD profiles, the NC of the Ca IRT) for this specific observing season.A2).The dashed lines are sine fits to the phase-folded data points, including the first harmonic in the left column and the fundamental only in the right one.

ZEEMAN-DOPPLER IMAGING OF CI TAU
From time series of Stokes and LSD profiles, one can reconstruct the magnetic field and brightness distribution at the surface of an active star like CI Tau using Zeeman-Doppler Imaging (ZDI), a tomographic technique inspired from medical imaging that allows one to invert phase-resolved sets of line profiles into maps of the large-scale vector magnetic field and of the photospheric brightness at the surface of the star.

Zeeman-Doppler Imaging
To achieve this goal, ZDI proceeds iteratively, starting from a small magnetic seed and a featureless brightness map, and progressively building up information by exploring the parameter space using a variant of the conjugate gradient method (e.g., Skilling & Bryan 1984;Brown et al. 1991;Donati & Brown 1997;Donati et al. 2006; Kochukhov 2016).At each iteration, ZDI adds features on the maps by comparing the observed Stokes profiles with the synthetic ones associated with the current images.The aim is to obtain a given agreement with the data (i.e., a given 2 r , usually 1) with minimal information in the maps to ensure that whatever is reconstructed is mandatorily required by the data.This reflects the fact that the problem is ill-posed and allows an infinite number of solutions of variable complexity, among which ZDI selects the simplest one that matches the data at the requested level of 2 r .When optical data are also available (as for our 2020 observing season), one can constrain the large-scale field and brightness map by simultaneously using Stokes and LSD profiles of photospheric lines from both domains, explicitly taking into account that the spot-to-photosphere brightness contrast differs at optical and nIR wavelengths.Furthermore, one can also reconstruct at the same time with ZDI a surface map of the accretion-induced excess emission in specific emission lines, the NC of the Ca IRT in our case, from time-series of their profiles (e.g., Donati et al. 2010Donati et al. , 2011)); such maps inform us on where, at the surface of the star, the local line emission, boosted by accretion, gets much larger than the quiet chromospheric emission.
In practice, we start by dividing the surface of the star into a grid of a few thousand (typically 5000) cells so that synthetic Stokes and profiles at each observation epoch can be computed by summing up the contributions of all cells in the reconstructed images, taking into account the assumed stellar parameters, in particular the inclination of the rotation axis to the line of sight (set to 70 • , see Table 1), the line-of-sight-projected equatorial velocity at the surface of the star sin (set to 9.5±0.5 km s −1 following Donati et al. 2020a), and the linear limb darkening coefficient (set to 0.3 and 0.7 at nIR and optical wavelengths respectively).We also assume that the surface of CI Tau rotates as a solid body, consistent with previous findings that fully-convective TTSs similar to CI Tau show weak levels of surface differential rotation (e.g., Finociety et al. 2023a).
Local Stokes and contributions from each cell are derived using Unno-Rachkovsky's analytical solution of the polarized radiative transfer equation in a plane-parallel Milne Eddington atmosphere (Landi degl'Innocenti & Landolfi 2004), assuming a Landé factor and average wavelength of 1.2 and 1750 nm for the SPIRou photospheric LSD profiles, and 1.2 and 640 nm for the ESPaDOnS ones.For the Doppler width D of the local profile, we assume D = 3 km s −1 , both for SPIRou and ESPaDOnS data, a value that yields synthetic LSD profiles consistent with those of CI Tau and other slowly rotating TTSs.The brightness maps that ZDI recovers consist of sets of independent pixels describing the local brightness (relative to the quiet photosphere and denoted ) over all surface grid cells.When using SPIRou and ESPaDOnS data simultaneously, is assumed to depend on wavelength in a way where S at 1750 nm can be simply inferred from E at 640 nm by S = E .If the spot brightness obeys Planck's law, we should have ≃ 0.40; if however the spot contrast depends much more steeply on wavelength (as for, e.g., the wTTS V410 Tau, Finociety et al. 2021), can end up being much smaller.
In the case of the Ca IRT profiles, we proceed as in previous similar papers (e.g., Donati et al. 2010Donati et al. , 2011Donati et al. , 2012)), i.e., by applying to the data a filtering procedure whose aim is to retain rotational modulation only and discard intrinsic variations, hence helping the convergence of the imaging code.The local profile we use to describe this line has a wavelength 850 nm, a Doppler width D = 7 km s −1 (ensuring synthetic profiles consistent with those of CI Tau) and a Landé factor 1.0.The quantity ZDI recovers, again handled as a set of independent pixels, describes the fractional area of each surface pixel in which accretion occurs, and where the local line emission flux, boosted by accretion, is assumed to be larger than that in the quiet surrounding chromosphere by an accretion-induced emission enhancement factor .As in our previous papers (e.g., Donati et al. 2010Donati et al. , 2011Donati et al. , 2012)), we set = 10 for the current study, with the exact value of having no more than a global scaling effect on the recovered images (ensuring that the integrated fractional accretion area times is constant).
Last but not least, the magnetic field at the surface of the star is described as a spherical harmonics (SH) expansion, using the formalism of Donati et al. (2006) in which the poloidal and toroidal components of the vector field are expressed with 3 sets of complex SH coefficients, ℓ, and ℓ, for the poloidal component, and ℓ, for the toroidal component 3 , where ℓ and note the degree and order 3 In this paper, we use modified and more consistent expressions for the field of the corresponding SH term in the expansion.In the case of CI Tau and given the low sin , we only use SH terms up to ℓ = 5, which in practice is enough for describing the relatively simple large-scale field topology that the star hosts.We further assume that the field is mostly antisymmetric with respect to the centre of the star (i.e. that odd SH modes dominate, as in, e.g., Donati et al. 2011), so that accretion funnels linking the inner disc to the star are anchored at high latitudes, in agreement with our observations (see Sec. 4).In practice, this is achieved by penalizing even SH modes with respect to odd ones in the entropy function.Finally, we assume that only a fraction of each grid cell (called filling factor of the large-scale field, equal for all cells) actually contributes to Stokes profiles, with a magnetic flux over the cell equal to (and thus a magnetic field within the magnetic portion of the cells equal to / ).Similarly, we assume that a fraction of each grid cell (called filling factor of the small-scale field, again equal for all cells) hosts small-scale fields of strength / (i.e., with a small-scale magnetic flux over the cell equal to = / ).This simple model implies in particular that the small-scale field locally scales up with the largescale field (with a scaling factor of / ), which is likely no more than a rough approximation, but at least ensures that the resulting Zeeman broadening from small-scale fields is consistent with the reconstructed large-scale field.For our study, we set ≃ 0.8 for all lines, consistent with actual measurements (Sokal et al. 2020), whereas we set ≃ 0.4 (again for all lines), slightly larger though still consistent with the value used in our previous study (i.e., 0.3, Donati et al. 2020a) and yielding optimal fits to the data.

ZDI from SPIRou data only
We start by applying ZDI to our SPIRou LSD profiles of CI Tau only, for each of our 3 observing seasons, i.e., 2019 (2019 October to December), 2020 (2020 September to 2021 January) and 2022 (2022 November to 2023 January).Note that the 6 profiles recorded in 2020 February (see Table A2) were left out of the 2019 season set, as they were collected 2 months after the bulk of the other 23 profiles, at a time where the large-scale field already started to evolve significantly according to the ℓ curve (see Fig. 1 top panel).The resulting ZDI fits to the LSD Stokes and profiles are shown in Fig. 3 whereas the resulting maps of the large-scale field are presented in Fig. 4. Stokes profiles are fitted at an average 2 r level of about 1.2 for all seasons, i.e. slightly larger than 1 as a likely result of the field being subject to intrinsic variability over each observing season, both on short (e.g., accretion) and medium (e.g., convection) timescales (see, e.g., Fig. 1).
Although brightness maps are reconstructed at the same time as magnetic maps when ZDI simultaneously fits Stokes and data, we find that CI Tau exhibits no obvious brightness surface feature that is large enough and / or with a high enough contrast in the nIR with respect to the surrounding photosphere to generate clear Stokes profile distortions (see Fig. 3) and thereby to show up in our reconstructed brightness images (not shown).This is in contrast with the brightness images in the optical domain that we reconstructed in our previous study, that showed high-contrast cool spots associated with strong magnetic field regions (Donati et al. 2020b).
From the reconstructed magnetic maps, we conclude that CI Tau hosts a large-scale magnetic field of average strength ≃0.8-0.9 kG components, where ℓ, is replaced with ℓ, + ℓ, in the equations of the meridional and azimuthal field components (see, e.   over the star at all times, reaching a maximum intensity of 0.9, 1.1 and 1.3 kG in seasons 2019, 2020 and 2022.Taking into account the assumed scaling factor / ≃ 2 between the small-scale and large-scale field (see Sec. 5.1), this implies a small-scale field of order 1.6-1.8kG, reaching up to 2.6 kG, in reasonable agreement with literature results (Sokal et al. 2020).The large-scale field we reconstruct is almost fully poloidal and axisymmetric, and mainly consists of a 0.8 kG dipole inclined at 11 • (in 2019) to 18 • (in 2020 and 2022) to the rotation axis, towards phases 0.1, 0.9 and 0.4 in seasons 2019, 2020 and 2022 respectively.The quadrupole and octupole components are weak, at most 0.3 kG for the latter in 2022.We summarize the magnetic properties of the reconstructed magnetic topologies in Table 3.
What it suggests already is that the large-scale field of CI Tau significantly changed since our previous study (Donati et al. 2020a), with the dipole field being weaker by potentially up to a factor of 2 and the octupole component having almost vanished (weaker by an order of magnitude).We cannot entirely rule out that part of this difference reflects the different data used in both study (SPIRou nIR data here and ESPaDOnS optical data in Donati et al. 2020a).However, looking at the new set of ESPaDOnS data we collected in season 2020 and in particular the ℓ curves (see Sec. 4 and Fig. 2), we can already confirm that significant changes occurred in the largescale of CI Tau (see Sec. 4).The following section gives further support in this direction.

ZDI from SPIRou and ESPaDOnS data
For our 2020 data set, we also applied ZDI to our combined SPIRou and ESPaDOnS data, using now three independent sets of lines (LSD photospheric profiles from both SPIRou and ESPaDOnS spectra, and the Ca IRT from the ESPaDOnS spectra) to characterize in a more self-consistent way the large-scale and small-scale field of CI Tau in this specific season.As pointed out in Sec. 4, the longitudinal field curves associated with our three sets of lines, despite varying in phase, are apparently inconsistent, with ESPaDOnS photospheric lines showing mainly positive longitudinal fields, the NC of Ca IRT lines exhibiting negative fields only, and SPIRou photospheric Table 3. Properties of the large-scale and small-scale magnetic field of CI Tau for our 3 observing seasons, using our SPIRou data only (columns 2 to 6) or SPIRou and ESPaDOnS data simultaneously (in 2020 only, columns 7 to 11).We list the average reconstructed large-scale field strength < > (columns 2 and 7), the maximum small-scale field strength (columns 3 and 8), the polar strength of the dipole component d (columns 4 and 9), the tilt / phase of the dipole field to the rotation axis (columns 5 and 10) and the amount of magnetic energy reconstructed in the poloidal component of the field and in the axisymmetric modes of this component (columns 6 and 11).Error bars are typically equal to 5-10% for field strengths and percentages, and 5-10  lines displaying both positive and negative longitudinal fields (see Figs. 1 and 2).The simplest qualitative explanation for this apparent discrepancy is that dark spots are present at the surface of the star at ESPaDOnS wavelengths, obscuring some of the most magnetic regions and preventing them from significantly contributing to Stokes signatures at optical wavelengths.However, as we have seen in Sec.5.2, brightness images reconstructed simultaneously with magnetic maps from SPIRou data show virtually no feature at the surface of the star, either because these features have a very low contrast with respect to the photosphere at SPIRou wavelengths, and / or because they are minor contributions to the line profile distortion (with respect to, e.g., magnetic fields).We explore both options below.
We start by assuming that spot brightness varies with wavelength as expected from Planck's law, and thus that , the exponent with which relative brightness at SPIRou wavelengths S can be inferred from that at ESPaDOnS wavelengths E , is equal to ≃ 0.4 (see Sec. 5.1).In this context, we are able to obtain a convincing fit to all three sets of lines (see Fig. B1).The corresponding ZDI maps are shown in Fig. 5, where we present the reconstructed brightness map at both ESPaDOnS and SPIRou wavelengths, along with the map of accretion-induced excess emission in the NC of the Ca IRT lines.We obtain that the reconstructed large-scale field reaches an averaged strength of 840 G and a maximum intensity of 1.2 kG, i.e., slightly larger than the values obtained from SPIRou data only (see Sec. 5.2).The derived magnetic topology is again mostly poloidal and axisymmetric, as previously indicated by SPIRou data.The dipole component of the large-scale field reaches almost 1.1 kG, i.e., ≃35 percent larger than our estimate from SPIRou data alone (thanks to the additional constraints from the Stokes and ESPaDOnS data) and is tilted at ≃15 • to the rotation axis towards phase 0.90.We also confirm that the quadrupole and octupole components of the large-scale field are weak, smaller than 150 G in season 2020.
The map of accretion-induced excess emission in the NC of the Ca IRT lines (see Fig. 5, bottom right) demonstrates that accretion occurs mostly towards a region that is centred on the pole and extends to intermediate latitudes, consistent with a low-amplitude modulation of emission peaking around phase 0.9 (see Fig. 2 middle right panel), i.e., the phase towards which the dipole component of the reconstructed field is tilted.It is also the phase at which the longitudinal field as measured in SPIRou LSD profiles of photospheric lines (see Fig. 1 and Table A2) and in the NC of the Ca IRT lines (see Fig. 2 and Table A2) is strongest, at least when Ca IRT ℓ values are averaged over the two rotation cycles we monitored with ES-PaDOnS.This is consistent with the footpoints of accretion funnels linking the stellar surface to the inner accretion disc being anchored at high latitudes on CI Tau, close to the region of maximum magnetic field.It also matches the longitudinal field measured in the NC of the He D 3 line (≃-1 kG, see Sec. 4), equal to the projection of the field vector in the accretion region where the average field is inclined at ≃55 • to the line of sight.
The reconstructed brightness image of CI Tau shown in Fig. 5 is predicted to generate photometric variations with a full amplitude of 4.5% in the ESPaDOnS domain ( band), and 1.8% in the SPIRou domain ( band), both of which would be hardly detectable amidst accretion-induced photometric perturbations causing intrinsic photometric fluctuations at a level of ≃0.1 mag (i.e., 10%) RMS for CI Tau (see Fig. 2 middle left panel, and Manick et al, submitted).It explains why the rotation period is hard to detect from a simple light-curve periodogram, especially in seasons where the large-scale magnetic topology, and thereby the brightness distribution at the surface of the star, is close to being axisymmetric, as was the case in our 2020 observing season.The derived maps are also consistent, by construction, with the observed RV variations in the various spectral lines used for the ZDI imaging, that do not exceed a peak-to-peak amplitude of 0.5 km s −1 for photospheric lines and 1 km s −1 for the NC of the Ca IRT (see, e.g., Fig. 2).We come back in Sec.7 for a more extensive discussion on the SPIRou RVs of CI Tau.As a sanity check, we also verified that the maps reconstructed from our 2020 ESPaDOnS data alone (using both Stokes and LSD profiles of photospheric lines and those of the NC of the Ca IRT) are consistent with those using both SPIRou and ESPaDOnS data (shown in Fig. 5).
Assuming now that is much smaller than 0.4 to ensure that the reconstructed brightness map at 1750 nm is as featureless as those derived in Sec.5.2 from SPIRou data alone, we find in practice that needs to be as small as 0.01 to obtain this result.This would mean in particular that brightness features at the surface of CI Tau at SPIRou wavelengths are drastically less contrasted than, and cannot be simply derived with Planck's law from those at ESPaDOnS wavelengths, which resembles the case of another TTS recently studied with SPIRou (namely V410 Tau, Finociety et al. 2021).The magnetic map we derive in this case (not shown) is similar to that obtained when assuming ≃ 0.4 (although with a ≃20% weaker dipole component), and so is the map of accretion-induced excess emission in the NC of the Ca IRT lines.The brightness map at ESPaDOnS wavelength is also similar to that shown in Fig. 5 (bottom left plot), whereas that at SPIRou wavelength is featureless (by design).The predicted relative photometric variations corresponding to both, featuring peak-to-peak amplitudes of only 6% and 0.5% respectively, would again be hardly detectable in the observed light curves of an actively accreting cTTS like CI Tau.
One way to distinguish between this second option (where most of the distortions in the LSD profiles of photospheric lines at SPIRou wavelengths are attributed to magnetic fields, as in Sec.5.2) and the first one (with ≃ 0.4, and where brightness spots also contribute to the LSD profile distortions) is to look at LSD profiles of CO molecular lines, in particular those of the 2.3 m CO bandhead, that are not sensitive to magnetic fields.Although the SNR in the LSD profiles associated with CO lines only (not shown) is about twice lower than that when all atomic lines are used, one can nonetheless distinguish profile distortions and variability that probe the presence of surface temperature features.For this reason, we conclude that our first option, corresponding to = 0.4 and to the maps shown in Fig. 5, is the most likely.This implies in particular that dark spots do contribute to LSD Stokes profile distortions, though only marginally in the SPIRou domain, hence why they do not show up in the brightness maps derived from SPIRou data only.We come back on this point in Sec. 7.
We also confirm our preliminary conclusion that the largescale field of CI Tau significantly changed since our first study (Donati et al. 2020a), and can further claim that the brightness distribution drastically evolved between the two epochs, from a dark region at photospheric level that is strongly off-centred with respect to the pole in early 2017 (generating a 4-5 km s −1 amplitude RV signature of the LSD Stokes profiles, e.g., Donati et al. 2020a) to a much more polar-centred spot configuration (inducing 10× smaller peakto-peak RV variations of 0.5 km s −1 , see Fig. 2 and Sec. 7).The same conclusion applies to the reconstructed maps of accretion-induced emission that is also much more centred on the pole during our 2020 observing season.Besides, our results unambiguously confirm that the RV signature with a 9-d period that CI Tau exhibited in late 2017 can be safely attributed to the activity of the star itself rather than to a close-in massive planet (Johns-Krull et al. 2016;Biddle et al. 2018Biddle et al. , 2021)), as further discussed in Sec. 7.

VEILING AND MASS ACCRETION RATE OF CI TAU
Veiling in the spectrum of CI Tau and other cTTSs observed with SPIRou was recently discussed at length in Sousa et al. (2023).Their conclusions from the same CI Tau spectra in seasons 2019 and 2020 are that average veiling in the , and bands is of order 0.4-0.5, whereas average veiling in the band is significantly larger, of order 0.9, with only limited temporal variations.They find in particular that the nIR veiling reflects dust heating in the inner disc, presumably boosted by accretion, hence why veiling is maximum in the K band over the SPIRou domain.The veiling we derive at optical wavelengths from our ESPaDOnS spectra in season 2020 is in the range 0.8-1.4(see Fig. 2 middle left plot), showing that nIR and optical veiling are different in nature and do not trace the same physical effects and spatial regions in the magnetosphere of CI Tau.We find that the optical veiling of CI Tau, tracing accretion at the surface of the star, does not peak when the accretion funnel crosses the line of sight (at phase 0.9, see Sec. 5.3); although surprising, this situation can happen when the accretion region is close to the pole (e.g., on DN Tau in 2012, Donati et al. 2013), as a likely result of the high level of intrinsic variability inherent to accretion processes.
Regarding the mass accretion rate of CI Tau, Sousa et al. ( 2023) conclude, from the equivalent widths of the Pa and Br lines, that log = −7.5 ± 0.2 M ⊙ yr −1 , consistent with the measurement from H and H lines in our previous ESPaDOnS campaign (log = −7.6±0.3 M ⊙ yr −1 , Donati et al. 2020a).For optical emission lines that show both a NC and a broad component (BC), namely the He 3 and the Ca IRT lines, the NC, tracing the material reaching the stellar surface at the footpoints of accretion funnels, has a much weaker equivalent width (by typically an order of magnitude, Donati et al. 2020a) than the BC, which is thought to trace material further out in the magnetosphere, in the inner disc or in the stellar wind.In the particular case of CI Tau where the BC of the He 3 line is significantly blue shifted with respect to, and uncorrelated with, the NC (Donati et al. 2020a), it further suggests that the BC actually traces a hot wind coming from the star itself and / or the inner disc, rather than in-falling material within the funnel flows (Beristain et al. 2001).In such a case, one may argue that only the NC of these lines, i.e., the one associated with funnel flows, should be taken into account to derive an estimate of , and doing so yields an average value of log = −8.5 ± 0.2 M ⊙ yr −1 for both our 2020 ESPaDOnS data and our older ones from early 2017 (using the same method as Donati et al. 2020a).Admittedly, empirical relations between emission line fluxes and accretion luminosities were derived from the whole lines and not from subcomponents, like the NC or BC of the He 3 and the Ca IRT lines; however, these relations may also overestimate in cases like CI Tau where the BC of both lines, and thus the Ba, Pa and Br lines as well, are potentially contaminated by a hot wind component from the star and / or the inner disc.It may for instance suggest that most of the inner disc material of CI Tau, potentially contributing to the BC of the He 3 and the Ca IRT lines and to the Ba, Pa and Br lines, is either ejected through magnetically-driven disc winds (e.g., Shu et al. 1994;Nisini et al. 2004) or magnetospheric ejection events (as in, e.g., Romanova et al. 2004;Ustyugova et al. 2006;Zanni & Ferreira 2013;Pantolmos et al. 2020), with only a small fraction reaching the surface of the host star and thereby showing up in the NC of the He 3 and the Ca IRT lines.The large mass loss rate recently reported from outflows and jets in the inner disc of CI Tau (log = −7.0M ⊙ yr −1 , Flores-Rivera et al. 2023), consistent with the mass accretion rate derived from the Pa and Br lines, comes as further support for our hypothesis.This picture, qualitatively consistent with the theoretical model of Lesur (2021), is however still speculative at this stage.In the following we simply consider that the two cases mentioned above correspond to an upper and a lower limit on the mass accretion rate at the surface of CI Tau.
From the simulation results of Bessolaz et al. (2008) regarding the radius of the magnetospheric gap mag that the large-scale magnetic field of CI Tau is able to carve at the centre of the inner accretion disc at the epoch of our observations, and using our estimate of the dipole component of the large-scale field of CI Tau as determined from our 2020 SPIRou and ESPaDOnS data, we find that mag / cor ≃ 0.35 for log ≃ −7.5 M ⊙ yr −1 and mag / cor ≃ 0.67 for log ≃ −8.5 M ⊙ yr −1 , with typical uncertainties on mag / cor of order 10%.Very similar results are obtained when using expressions of the magnetospheric radius from other studies (e.g., Blinova et al. 2016Blinova et al. , 2019;;Pantolmos et al. 2020).In the first case, this ratio is clearly too small for the large-scale magnetic field to trigger the observed stable poleward accretion, and succeed in spinning down CI Tau to its slow rotation rate (according to simulations, e.g., Zanni & Ferreira 2013;Blinova et al. 2016;Pantolmos et al. 2020); in the second case however, the large-scale field is just strong enough, yielding a better agreement with observations.Besides, recent simulations suggest that tilts in the large-scale field of the host star, like those we report for CI Tau, are able to induce warps of the inner accretion disc, but not a significant change in the magnetospheric size (e.g., Romanova et al. 2021).Reconciling our upper limit on the mass accretion rate with simulation predictions would require CI Tau to host a magnetic field with a dipole component as strong as 3.5 kG, which is clearly inconsistent with existing observations.Another option is that simulations are grossly off, e.g., lacking one critical physical ingredient; however, it seems unlikely given their relative success at reproducing a variety of observational cases (e.g., Blinova et al. 2019).
Looking at the temporal variability of the 1083 nm He , Pa and Br lines of CI Tau over the full set of SPIRou observations yields further information about how accretion proceeds within the magnetosphere of CI Tau.The 1083 nm He line features a more or less centred emission component with a full-width-at-half-maximum (FWHM) of ≃250 km s −1 , along with a blue-shifted absorption component extending down to -350 km s −1 (see Fig. 6 top left panel).An additional red-shifted pseudo-absorption component (always staying above a relative flux level of 1.0) also shows up at velocities 0-120 km s −1 .All three components are found to vary with time.The corresponding 2D periodogram (see Fig. 6 bottom left panel) features a clear peak at a period of about 9.05 d, i.e., slightly larger than that found from ℓ data, in conjunction with the red-shifted absorption component.This red-shifted component most likely traces the outer part of the accretion funnel linking the inner disc to the visible polar regions at the stellar surface of CI Tau, as it periodically crosses the line of sight.Since no clear 9-d periodicity shows up in association with the blue-shifted absorption component of the He line, which would argue in favour of a stellar wind, we conclude that it presumably probes a wind from the inner regions of the accretion disc (despite the overall shape of the He line may suggest a stellar wind, e.g., Kwan et al. 2007).In fact, the 2D He periodogram only shows power in the blue wing for periods larger than ≃17 d (see Fig. C1 left panel), suggesting a structured disc wind from radii 0.13 au and beyond.
The Pa line consists mainly of a central emission of FWHM ≃250 km s −1 (as for the central emission of the 1083 nm He line) with a pseudo-absorption feature centred at velocities of about 0-120 km s −1 (see Fig. 6 top right panel).The 2D periodogram shows that it apparently varies with the same period as the red-shifted absorption component of the 1083 nm He line (see Fig. 6 bottom right panel) and thereby likely traces as well the accretion funnel linking the inner disc to the star as it crosses the line of sight.Another weaker blue-shifted component at velocities of about -200 km s −1 also varies with a period of about 8.95 d, i.e., slightly shorter than that found from the ℓ data, and possibly with a period of 9.05 d as well (although some of the corresponding power, but not all, may reflect aliasing).As for He , the 2D periodogram shows most power in the blue wing of Pa for periods larger than ≃17 d (see Fig. C1 right panel), suggesting again a structured disc wind from radii 0.13 au and beyond.Finally, we find that Br (not shown) exhibits an FWHM and temporal behaviour similar to those of Pa but with more noise, making it harder to diagnose the main components of the line profile and their corresponding periodicities.
Altogether, these lines indicate that the outer portions of accretion funnels (generating the red-shifted absorption component) are modulated on a timescale that is slightly longer than the rotation period measured from ℓ data, possibly because they are anchored in a region of the inner disc that rotates a bit more slowly than the star and thereby forces accretion funnels to trail stellar rotation; however, this would contradict our previous finding that mag is unlikely to be larger than cor given the derived magnetic field and mass accretion rates.Another option is that CI Tau exhibits surface differential rotation, with the polar regions, where accretion funnels are anchored, rotating more slowly than the average rotation rate at the surface of the star (as measured from ℓ ).That the weak blue-shifted pseudo-absorption component of Pa , presumably tracing a stellar wind, is also modulated, with two different periods that are slightly shorter and slightly longer than the average stellar rotation period respectively, may further support this second option and suggest the presence of a weak accretion-powered stellar wind emerging from both low and high latitudes at the surface of CI Tau, in a way sim- ilar to coronal streamers and coronal holes in the case of the solar wind.This option would imply a surface differential rotation of order 10 mrad d −1 at the surface of CI Tau, consistent with the low differential rotation latitudinal shear measured in other fully convective TTSs (e.g., Finociety et al. 2023a) and justifying independently our choice of assuming solid body rotation in the imaging section of our study (see Sec. 5).
More generally, the 9-d periodicity that the He and Pa lines exhibit on a time scale of 1200 d in the red wing (and the blue wing as well in the case of Pa ) demonstrates that magnetospheric accretion on CI Tau proceeds in a relatively stable manner along a more or less stable magnetic topology, rather than in an unstable fashion through chaotic and rapidly evolving accretion tongues crossing field lines towards the stellar equator (e.g., Blinova et al. 2016), as no clear period would otherwise emerge from the 2D periodograms.This further supports our previous conclusion that mag / cor must be larger than 0.6 in CI Tau.
We show in Fig. D1 the stacked H and H profiles of CI Tau collected during our late 2020 ESPaDOnS run, which also show clear pseudo-absorption features in the red and blue wings of the profiles, qualitatively similar to what we see in Pa .

RADIAL VELOCITY MODELING CI TAU
Finally, we examine RVs of CI Tau derived from Gaussian fits to the Stokes LSD profiles of our SPIRou spectra.We end up with 94 values covering an overall timescale of about 1200 d, both for our main LSD profiles extracted from atomic lines throughout the whole SPIRou domain, and for those derived using the CO bandhead mask (see Sec. 2).We use simple Gaussian fits to LSD profiles to extract RVs in this study, rather than the line-by-line approach (Artigau et al. 2022) shown to behave optimally for stars whose spectra vary modestly with respect to the median, but to be less suited for accreting cTTSs that exhibit large spectral variations with time (including veiling, see Sec. 6).We find that CI Tau features a median RV of 16.7 km s −1 with an RMS dispersion of 0.3 km s −1 in atomic lines, and a median RV of 17.1 km s −1 with an RMS dispersion of 0.5 km s −1 in CO bandhead lines.The small RV shift between the two sets of lines may reflect that on average, atomic lines are not formed at the same atmospheric depth as the CO bandhead lines, and therefore do not probe the same average convection velocities.In particular, we speculate that the downward flows in the falling lanes of convective cells are relatively brighter for CO lines than for atomic lines (with respect to the upward flows in the rising granules), thereby inducing a small red shift between the two sets of lines.The larger RV dispersion of CO lines largely reflects the lower SNR of the corresponding LSD profiles (with respect to those from atomic lines).
We analyse these RVs in the same way as in our previous studies (e.g., Donati et al. 2023), looking for the impact of modulated activity on the RVs but also for the potential contribution of a close-in massive planet that may be orbiting around CI Tau (at a different period than that originally proposed by Johns-Krull et al. 2016, as we now know that the 9-d period traces stellar rotation).Following Manick et al. (submitted) who find evidence for a close-in massive planet on a 24-d orbit around CI Tau, we look whether an RV signal may be present in the data, in particular at orbital periods close to 24 d.We achieve this by adjusting our RV measurements using GPR to model the contribution from stellar activity (as in Sec.4), and a sine wave (or a Keplerian model if needed) to describe the reflex motion induced by a potential close-in planet.We then run a Monte Carlo Markov Chain (MCMC) experiment to derive optimal values, error bars and posterior distributions for all model parameters, in the cases where only activity is taken into account and where both activity and a close-in planet are considered.In both cases, the evolution timescale ( 3 ) and the smoothing parameter ( 4 ) are fixed to typical values of 200 d and 0.5, to limit the flexibility of GPR given that the error bars on both sets of RVs are only a few times smaller than their RMS dispersion.The results of the MCMC runs are detailed in Table 4.
We find that stellar activity is contributing to the RV curve at a semi-amplitude level of 0.30 km s −1 for atomic lines and 0.20 km s −1 for CO bandhead lines, demonstrating that magnetic effects indeed dominate profile distortions of atomic lines, as anticipated in Sec. 5.It also demonstrates that, although smaller, brightness variations at the surface of CI Tau contribute to profile distortions as well, in particular for CO bandhead lines that are not sensitive to magnetic fields.This further supports our conclusion that CI Tau hosts cool spots at its surface that are affecting photospheric line profiles, mostly at ESPaDOnS wavelengths but also in the SPIRou domain, consistent with our previous conclusion (see Sec. 5.3).We note that the period with which the activity jitter is modulated (equal to 8.92±0.02d in the case of atomic lines, and a similar value with a larger error bar in the case of the CO bandhead lines, see Table 4) is slightly but significantly lower than that derived from ℓ data.As in Sec.6, we speculate that this difference reflects latitudinal differential rotation at the surface of CI Tau, with regions contributing most to the activity jitter being those located close to the equator, i.e., at a lower latitude than those shaping the ℓ curve.Incidently, we find that the FWHMs of the LSD profiles of atomic (varying from 21.4 to 25.0 km s −1 ) are also rotationally modulated (the GPR fit yielding an RMS of 0.5 km s −1 ), with a clear period of 8.93 ± 0.04 d, consistent with that derived from the associated RVs, and a semi-amplitude of 0.7 ± 0.2 km s −1 .Besides, we find that these RVs correlate reasonably well with the first derivative of the FWHMs given the significant level of noise on both quantities (correlation coefficient ≃ 0.5), further demonstrating that magnetic activity is indeed causing the 9-d RV modulation.
We also examined the possible presence of a close-in massive planet through the reflex motion and the RV signature it may induce on CI Tau itself, in particular at the orbital period at which a tentative signal has recently been reported (i.e.,24 d,Manick et al.,submitted) after being initially spotted in photometric data but finally labelled as a likely alias (Biddle et al. 2018).The marginal logarithmic likelihood log L of a given solution is computed using the approach of Chib & Jeliazkov (2001) as described in Haywood et al. (2014), and the significance of the RV signature of a putative planet is estimated from the difference in log L , i.e., the logarithmic Bayes Factor log BF, with respect to our reference model with no planet.We find that atomic lines show no more than a hint of a signal (at a 2 level) at this period, with an amplitude of only = 0.06 +0.04 −0.03 km s −1 (corresponding to a minimum planet mass of sin = 0.79 +0.52 −0.39 M if the RV signal is caused by a planet) and a log BF of only 1.9, i.e., clearly too small for this signal to be even qualified as tentative.Carrying out the same experiment on the RVs from the CO bandhead lines, we obtain a different result, this time with a clear 5.6 signal of semi-amplitude = 0.28 +0.06 −0.05 km s −1 (corresponding to a minimum planet mass of sin = 3.70 +0.80 −0.66 M ) at a well-determined orbital period of = 23.86 ± 0.04 d, and associated with a log BF of 13.1 apparently indicating a firm detection.The resulting fit to the raw RVs is shown in Fig. 7 (top panel) along with the candidate planet signal in the filtered RVs (middle) and the residual RVs (bottom, RMS of 0.32 km s −1 ), whereas the RV curve phase-folded on the recovered orbital period is displayed in Fig. 8. Assuming a Keplerian rather than a circular orbit to model the planet signature yields an eccentricity consistent with 0 (with an error bar of about 0.25) and does not improve log BF in a significant way (Δ log BF ≃ 0.1), justifying a posteriori our hypothesis of a circular orbit.We point out that, although similar at first glance, the characteristics of the detected RV signal inferred from our study significantly differ from those of Manick et al. (submitted), whose results we were not able to reproduce, even from their own data.We finally note that no clear signal shows up in the 2D periodograms of He , Pa and Br lines at the period of the detected RV signature, apart from a weak one in the blue wing of He (see Fig. C1 for the He and Pa lines).
At this point, it is not clear why the detected signal, if induced by an orbiting planet, would be detected in the RVs of the CO bandhead lines but not in those of atomic lines, despite the latter being more precise and less dispersed around the GPR fit (0.18 km s −1 RMS) than the former (0.32 km s −1 RMS).If probing the presence of a close-in planet, the RV signal should indeed show up in the same way in both sets of lines, and not just in one of them.Admittedly, CO bandhead lines are less affected by activity, in particular magnetic fields; however, the RV precision that they yield is twice lower than that from atomic lines, implying that a putative planet signal that shows up in CO lines should also be detected with a similar semiamplitude in atomic lines.We speculate that the detected RV signal rather probes the presence of a non-axisymmetric density structure in the accretion disc, located at an average distance of ≃0.16 au where the Keplerian period is consistent with that of the detected RV signal (23.86 ± 0.04 d), and distorting the Stokes LSD profiles of atomic and CO bandhead lines of CI Tau in a different way as a result of, e.g., the disc structure being much cooler than the photosphere of CI Tau.Although the Keplerian velocity at a distance of 0.16 au (i.e., 57.4 km s −1 ) is in principle 2.5× and 3.8× larger than the FWHM of the LSD Stokes profiles of CI Tau in atomic and CO bandhead lines respectively, some of the disc material may rotate at a strongly sub-Keplerian velocity as a result of, e.g., magnetic fields (Donati et al. 2005;Shu et al. 2007), thereby contributing to the disc profile with a component whose FWHM is closer to that of the stellar lines.As this non-axisymmetric disc structure orbits around the star, its spectral contribution, adding up to the stellar spectrum, is expected to modulate RVs of LSD profiles in a way that depends on the selected lines, qualitatively similar to what we observe.The nature of this axisymmetric structure is unclear, but could reflect the presence of strong magnetic fields in the disc, as in the case of FU Ori (Donati et al. 2005); it may also relate to the presence of a vortex (Varga et al. 2021), and / or a nascent planet, with a mass low enough not to be detected in atomic lines but large enough to induce a non-axisymmetric density structure in the accretion disc.
We note that the beat periods between the 23.86-d RV signal and the 9.01-d rotation period of CI Tau, equal to 6.5 and 14.5 d, fall close to the secondary peaks in the periodogram of the 2017 K2 light curve of CI Tau as reported by Biddle et al. (2018, located at 6.6 and 14.3 d).It may suggest that the putative disc structure at 0.16 au that we invoke to explain the RV signal in CO bandhead lines may also participate to star-disc interactions, e.g., generating enhanced accretion each time it falls in conjunction with the accretion funnels that rotate with the star (e.g., as in the pulsed accretion model of Teyssandier & Lai 2020).
Table 4. MCMC results for the 4 studied cases (no planet and planet b, using RVs from atomic lines and from CO bandhead lines).For each case, we list the recovered GP and planet parameters with their error bars, as well as the priors used whenever relevant.The last 4 rows give the 2 r and the RMS of the best fit to our RV data, as well as the associated marginal logarithmic likelihood log L and marginal logarithmic likelihood variation Δ log L with respect to the corresponding model without planet.

SUMMARY AND DISCUSSION
In  2018,2021), suggesting that this alternate period may relate to some other phenomenon of the complex star-disc interactions taking place in CI Tau.We observe that ℓ curves differ in both wavelength domains, with nIR photospheric lines probing longitudinal fields varying from -170 to 120 G over our three observing seasons (including 2020), whereas optical photospheric lines yield ℓ values only fluctuating from -5 to 130 G in late 2020.It indicates the presence of brightness features at the surface of CI Tau with a contrast that depends on wavelength, thereby affecting differently optical and nIR spectra.Modeling the Stokes and profiles of SPIRou spectra with ZDI allows us to reconstruct the distribution of brightness and magnetic features at the surface of CI Tau, indicating that CI Tau hosts a mainly poloidal and axisymmetric large-scale field whose average strength ranges from 0.8 to 0.9 kG, and that essentially consists of a slightly tilted dipole with a polar strength of 0.8 kG.Adding in the ESPaDOnS data for the ZDI modeling of our 2020 observations provides stronger constraints on the field topology, the associated brightness map and the distribution of the accretion-induced excess emission in Ca IRT lines.It tells that accretion mainly proceeds towards the pole at the surface of CI Tau, with the footpoints of accretion funnels showing up as bright features at chromospheric level and dark features at photospheric level.The updated magnetic topology is similar to that derived from SPIRou data only, except for a stronger dipole component (to compensate for the relative darkness of polar regions) of strength 1.1 kG.We also report a clear change of both the magnetic topology and brightness map of CI Tau with respect to our previous ESPaDOnS observations in early 2017 (Donati et al. 2020a), with the dark features at the footpoints of accretion funnels being now much less off-centred with respect to the pole and thereby generating much smaller profile distortions and RV variations (in both optical and nIR domains).
Both the 1083 nm He and Pa lines show a clear modulation at a 9-d period in the red wing of their profile, tracing the base of the accretion funnel close to the inner disc crossing the line of sight as the star rotates.The blue wing of Pa is also modulated by rotation, possibly probing a weak accretion-powered stellar wind escaping from both low and high latitude regions at the surface of the star.These modulation periods, slightly different from that of ℓ data, suggest that surface differential rotation is weak at the surface of CI Tau, consistent with previous results on fully convective TTSs (e.g., Finociety et al. 2023a).Altogether, it indicates that accretion is stable enough at the surface of CI Tau to generate a clear signal in the 2D periodograms of both lines, since unstable accretion as described by, e.g., Blinova et al. (2016), with disc material penetrating the magnetosphere through chaotic equatorial accretion tongues, is unlikely to produce a coherent signal on a timescale of 1200 d.The mass accretion rate inferred from the Pa and Br lines of the same SPIRou spectra by Sousa et al. (2023) is equal to log = −7.5 ± 0.2 M ⊙ yr −1 , similar to that derived from Ba lines in optical spectra (Donati et al. 2020a).In the particular case of CI Tau however, the BC of the He 3 line is significantly blue shifted, suggesting that it probes a hot wind from the star and / or the inner accretion disc rather than in-falling material in funnel flows (Beristain et al. 2001).Considering the NC of both features only yields log = −8.5 ± 0.2 M ⊙ yr −1 for CI Tau in 2020, which we consider as a lower limit to the mass accretion rate, as opposed to the upper limit derived from Ba, Pa and Br lines, or from the combined NC and BC of the He 3 and Ca IRT lines.We find that this upper limit implies a magnetospheric size (relative to the corotation radius) of about mag / cor ≃ 0.35 given the strength we measure for the large-scale dipole field of CI Tau, whereas the lower limit yields mag / cor ≃ 0.67.The former estimate appears inconsistent with requirements for stable accretion as derived from numerical simulations (e.g., Blinova et al. 2016;Pantolmos et al. 2020), whereas the latter is in better agreement.Our results may thus imply that only a small fraction of the disc material traced with Ba, Pa and Br lines reaches the magnetosphere of CI Tau (as in, e.g., Lesur 2021) and makes it to the surface of the star in a stable accretion pattern, where it shows up as the NC of the He 3 and Ca IRT lines, whereas the BC of both lines as well as most of material traced with Ba, Pa and Br lines would probe a wind from the inner disc and / or the star.The large mass loss rate recently reported from outflows and jets in the inner disc of CI Tau (Flores-Rivera et al. 2023), consistent with the mass accretion rate derived from the Pa and Br lines, comes as further support for this hypothesis.We however stress that this picture is still speculative and requires confirmation from further observations.
We also studied RVs from SPIRou spectra of CI Tau, analyzing both LSD profiles of atomic lines and those of CO bandhead lines, and using GPR to model the activity jitter.We find that both sets of RVs are modulated with a period slightly smaller than that found from the ℓ curve, again suggesting that weak surface differential rotation is present at the surface of the star.Besides, the semi-amplitude of the modulation, respectively equal to 0.3 and 0.2 km s −1 for atomic and CO lines, is inconsistent with the 1 km s −1 semi-amplitude reported for the putative close-in planet (Johns-Krull et al. 2016), further demonstrating that the 9-d RV signature is indeed of stel-lar rather than planetary origin.Looking for RV signatures potentially adding up on top of the activity jitter at a different period, and in particular at the 24 d period recently reported to show up in CI Tau photometric light curves and SPIRou and ESPaDOnS data (Manick et al, submitted), we find that atomic lines exhibit no more than a weak, barely significant, 2 signature of semi-amplitude = 0.06 +0.04 −0.03 km s −1 , whereas CO lines show a clear 5.6 signal of semi-amplitude = 0.28 +0.06 −0.05 km s −1 , the latter at a well-defined period of 23.86±0.04d.Atomic lines being more affected by activity than CO lines (hence the 1.5× larger semi-amplitude associated with the activity jitter), one can expect the former to be less sensitive to the planet signal than the latter; however, RVs from atomic lines are about twice more precise than those from CO lines, which should compensate for the lower sensitivity.For this reason, we speculate that the observed 23.86-d RV signal rather probes the presence of a non-axisymmetric structure in the accretion disc, possibly linked with strong magnetic fields, a vortex, and / or of a lower-mass planet, at an approximate distance of 0.16 au from the centre of the host star.This is at least qualitatively consistent with the latest interferometric measurements with GRAVITY (Gravity Collaboration et al. 2023), suggesting the presence of a gap ranging from 0.1 to 0.2 au in the inner disc of CI Tau, and with the study of CO lines in the M band (at 4.7 m, Kozdon et al. 2023), also arguing for a non-axisymmetric disc structure and a potential gap at a distance of about 0.14 au.
Our study brings us a step forward in the understanding of the complex magnetospheric accretion processes at work in the stardisc interactions of CI Tau, and more generally on their role in star / planet formation.More analyses and observations of the same type are nonetheless needed to progress further (e.g., Bouvier et al. 2023), and determine more accurately how much of the inner disc material is being transferred to the magnetosphere or expelled outwards through disc winds.Similarly, we need to understand in a more quantitative way the physical nature of the non-axisymmetric disc structure inducing the observed RV signal detected in CO lines but not in atomic lines of CI Tau.Looking in more details at the profiles of CO lines and how they vary with time, as in Flagg et al. (2019) but for the RV signal we detect in the SPIRou data, should tell us more about what this enigmatic disc structure consists of; this is already planned for a forthcoming companion paper.Our study also demonstrates the extreme value of collecting (quasi) simultaneous optical and nIR spectropolarimetric data of cTTSs.This should be feasible soon with the VISION project at CFHT, that will allow both SPIRou and ESPaDOnS to observe the same object with the two instruments, maintaining accurate polarimetric information in both channels.Similarly, coordinated spectropolarimetric and interferometric observations, like those recently carried out with GRAVITY (Gravity Collaboration et al. 2023) but with a more extensive monitoring of the rotation cycle, should also bring key complementary information about how magnetospheric accretion proceeds in CTTSs like CI Tau.

Figure 1 .
Figure1.Longitudinal magnetic field ℓ of CI Tau (red dots) as measured with SPIRou throughout our campaign, and QP GPR fit to the data (cyan full line) with corresponding 68% confidence intervals (cyan dotted lines).The residuals are shown in the bottom plot of each panel.The top panel shows the whole data set, whereas the lower 2 panels present a zoom on the 2020 and 2022 data respectively.The RMS of the residuals is about 14 G, consistent with our median error bar of 15 G, yielding 2 r = 0.82 (whereas the 2 r with respect to the ℓ = 0 G line is 22).

Figure 2 .
Figure 2. Variability and modulation of the ESPaDOnS LSD profiles (left) and of the NC of the Ca IRT (right) lines, as a function of rotation phase (computed as mentioned in Sec. 2).Each panel shows the RVs (top), the veiled EWs (or veiling in the case of LSD profiles, middle) and the ℓ values (bottom).The red circle, green squares and blue triangles depict measurements obtained during rotation cycles 45, 46 and 47 (see TableA2).The dashed lines are sine fits to the phase-folded data points, including the first harmonic in the left column and the fundamental only in the right one.

Figure 3 .
Figure 3. Observed (thick black line) and modelled (thin red line) LSD Stokes (top panels) and (bottom) profiles of CI Tau from SPIRou data, for seasons 2019 (left), 2020 (middle) and 2022 (right).Rotation cycles (counting from 0, 38 and 124 for the 2019, 2020 and 2022 seasons respectively, see TableA1) are indicated to the right of all LSD profiles, while ±1 error bars are added to the left of Stokes signatures.

Figure 4 .
Figure 4. Reconstructed maps of the large-scale field of CI Tau (left, middle and right columns for the radial, azimuthal and meridional components in spherical coordinates, in G), for seasons 2019, 2020 and 2022 (top to bottom row respectively), derived with ZDI from the SPIRou Stokes and LSD profiles of Fig. 3.The maps are shown in a flattened polar projection down to latitude −60 • , with the north pole at the centre and the equator depicted as a bold line.Outer ticks indicate phases of observations.Positive radial, azimuthal and meridional fields respectively point outwards, counterclockwise and polewards.

2020Figure 5 .
Figure 5. Same as Fig. 4 for the 2020 season, using now all three sets of lines from ESPaDOnS and SPIRou data and assuming = 0.4 (see Sec. 5.3).Below the three components of the magnetic field (top row), we show the reconstructed brightness map as seen at ESPaDOnS wavelengths (left) and SPIRou wavelengths (middle), along with the map of accretion-induced excess emission in the NC of the Ca IRT (bottom right).Darker color levels indicate darker photospheric features for the bottom left and middle plots, and brighter chromospheric features for the bottom right plot.

Figure 6 .
Figure 6.Spectra of CI Tau in the region of the He triplet (top left) and Pa line (top right), in the stellar rest frame, and the corresponding 2D periodograms (bottom).In the two upper panels, the curves show the superposition of all individual spectra, whereas the red vertical dotted lines depict the location of the He triplet and of the Pa line.The horizontal white line in the bottom panels indicate the rotation period measured from ℓ data.

Figure 7 .
Figure 7. Raw (top), filtered (middle) and residual (bottom) RVs of CI Tau from CO bandhead lines (red dots) over our SPIRou observing campaign of 1200 d.The top plot shows the MCMC fit to the RV data, including a QP GPR modeling of the activity and the RV signatures of a putative close-in planets of orbital period 23.86 ± 0.04 d (cyan), whereas the middle plot shows the planet RV signature alone once activity is filtered out.The RMS of the RV residuals is 0.32 km s −1 .

Figure 8 .
Figure 8. Filtered (top plot) and residual (bottom plot) RVs for the 23.86-d RV signal detected in CO bandhead lines of CI Tau.The red dots are the individual RV points with the respective error bars, whereas the black stars are average RVs over 0.1 phase bins.As in Fig. 7, the dispersion of RV residuals is 0.32 km s −1 .

Figure B1 .
Figure B1.Observed (thick black line) and modelled (thin red line) LSD Stokes (top panels) and (bottom) profiles of CI Tau, for the photospheric lines in the ESPaDOnS domain (left), the Ca IRT lines (middle) and the photospheric lines in the SPIRou domain (right) for our 2020 observing season.Rotation cycles (counting from 38, see TableA2) are indicated to the right of all LSD profiles, while ±1 error bars are added to the left of Stokes signatures.

Figure C1 .Figure D1 .
Figure C1.2D periodograms of the 1083 nm He (left) and Pa (right) lines for periods ranging from 6 to 30 d.The white horizontal lines depict periods the rotation period derived from our ℓ data (9.01 d), the period of the RV signal detected from CO bandhead lines (23.86 d), and the additional periods inferred from the K2 light-curve (i.e., 6.6 and 14.3 d Biddle et al. 2018, 2021)    .

Table 1 .
Parameters of CI Tau used in / derived from our study

Table 2 .
Results of our MCMC modeling of the ℓ curve of CI Tau with QP GPR.For each hyper-parameter, we list the fitted value, the corresponding error bar and the assumed prior.The knee of the modified Jeffreys prior is set to , i.e., the median error bar of our ℓ measurements (i.e., 15 G).We also quote the resulting 2 r and RMS of the final GPR fit.
• for field inclinations.
(Johns-Krull et al. 2016;bservations of the cTTS CI Tau collected with the SPIRou nIR spectropolarimeter / velocimeter at CFHT in the framework of the SLS and SPICE large programmes, complemented with PI observations recorded with the ESPaDOnS optical spectropolarimeter.A total of 94 and 13 validated Stokes and spectra of CI Tau were collected with SPIRou and ESPaDOnS, over a period of about 1200 d spanning three observing seasons for SPIRou(2019, 2020 and 2022), and over 16 d in late 2020 for ESPaDOnS.LSD profiles were computed for all spectra, yielding ℓ values that are unambiguously modulated by the rotation of the star, equal to 9.010 ± 0.023 d.We can thus safely exclude the alternate rotation period (of 6.6 d) that was initially proposed as the rotation period of CI Tau in several studies(Johns-Krull et al. 2016; Biddle et al.