SPIRou spectropolarimetry of the T Tauri star TW Hydrae: magnetic ﬁelds, accretion and planets

In this paper we report near-infrared observations of the classical T Tauri star TW Hya with the SPIRou high-resolution spectropolarimeter and velocimeter at the 3.6-m Canada-France-Hawaii Telescope in 2019, 2020, 2021 and 2022. By applying Least-Squares Deconvolution (LSD) to our circularly polarized spectra, we derived longitudinal ﬁelds that vary from year to year from –200 to +100 G, and exhibit low-level modulation on the 3.6 d rotation period of TW Hya, despite the star being viewed almost pole-on. We then used Zeeman-Doppler Imaging to invert our sets of unpolarized and circularly-polarized LSD proﬁles into brightness and magnetic maps of TW Hya in all 4 seasons, and obtain that the large-scale ﬁeld of this T Tauri star mainly consists of a 1.0–1.2 kG dipole tilted at about 20 ◦ to the rotation axis, whereas the small-scale ﬁeld reaches strengths of up to 3-4 kG. We ﬁnd that the large-scale ﬁeld is strong enough to allow TW Hya to accrete material from the disc on the polar regions at the stellar surface in a more or less geometrically stable accretion pattern, but not to succeed in spinning down the star. We also report the discovery of a radial velocity signal of semi-amplitude 11 . 1 + 3 . 3 − 2 . 6 m s − 1 (detected at 4.3 𝜎 ) at a period of 8.3 d in the spectrum of TW Hya, whose origin may be attributed to either a non-axisymmetric density structure in the inner accretion disc, or to a 0 . 55 + 0 . 17 − 0 . 13 M X candidate close-in planet (if orbiting in the disc plane), at an orbital distance of 0 . 075 ± 0 . 001 au.


INTRODUCTION
It is now well established, through documented observations collected with various instruments over the last decades, that stars and their planets form at the same time, following the collapse of large gravitationally unstable turbulent molecular clouds.The cloud collapse results in an accretion disc that feeds the central protostar ★ E-mail: jean-francois.donati@irap.omp.eufrom its inner regions, and where growing protoplanets form from merging planetesimals, giving birth, in the outer disc regions, to massive planets that can migrate inwards into close-in hot Jupiters (e.g., Baruteau et al. 2014).Magnetic fields play a key role in this process in many different ways, e.g., by hampering fragmentation within the disc (Hennebelle & Teyssier 2008), by evacuating the central regions of the accretion disc and forcing disc material to flow along discrete magnetospheric funnels linking the inner disc to the surface of the host star (e.g., Bouvier et al. 2014), by ex-tracting angular momentum outwards through outflows and jets and forcing the central star to slow down through star / disc interactions (e.g., Romanova et al. 2002;Zanni & Ferreira 2013), or by making inward-migrating giant planets pile-up at the outer edge of the magnetosphere and thereby saving them from falling into their host stars (e.g., Lin et al. 1996;Romanova & Lovelace 2006;Mulders et al. 2015).
T Tauri stars (TTSs), and in particular classical TTSs (cTTSs 1 ) that still accrete material from their accretion discs, are ideal objects to investigate these critical steps of star / planet formation, and to yield observational constraints on the complex physics at work, especially in the inner regions that drive star / disc / planet interactions and where the most energetic phenomena take place.However, although accretion discs of cTTSs are presumably actively forming planets given the radial structuring of their density profiles (e.g., Clarke et al. 2018), detecting planets around cTTSs has proven quite complex, either because of limited performances in angular resolution and contrast for direct imaging searches, or due to the extreme level of intrinsic variability cTTSs are subject to (e.g., Cody et al. 2014;Sousa et al. 2016) that drastically limit the precision of indirect velocimetric measurements.So far, distant planets have only been reliably detected around a single cTTSs (Keppler et al. 2018;Haffert et al. 2019), and claims of close-in hot Jupiters detected through velocimetry often turned out to be false positives attributable to activity (Setiawan et al. 2008;Huélamo et al. 2008;Johns-Krull et al. 2016;Donati et al. 2020a).Moreover, planet detection through transit events looks hopeless given the small size of expected transits, even for massive planets, with respect to the huge intrinsic variability induced by ongoing accretion processes between the disc and the central star.
Located at the heart of the TW Hydra (TWA) association about 60 pc away from the Sun, TW Hya, aged about 8-10 Myr (e.g., Torres et al. 2008;Luhman 2023), is the closest and most studied cTTS (e.g., Rucinski & Krautter 1983;Kastner et al. 2002;Herczeg et al. 2023) that hosts a large and massive accretion disc which survived longer than the typical disc dissipation timescale (of about 3 Myr, e.g., Kraus et al. 2012), and features rings and gaps suggesting ongoing planet formation (van Boekel et al. 2017).TW Hya is also known to harbour a strong large-scale magnetic field (Yang et al. 2007;Donati et al. 2011) that interacts with the inner accretion disc, as well as small-scale fields locally reaching up to several kG (e.g., Sokal et al. 2018;Lavail et al. 2019;López-Valdivia et al. 2021).TW Hya is thus an ideal laboratory to scrutinize magnetospheric accretion processes and more generally star / planet formation and interactions in the inner discs of cTTSs.
In this paper, we report extended observations of TW Hya with the near-infrared (nIR) high-resolution cryogenic spectropolarimeter / velocimeter SPIRou installed at the Cassegrain focus of the Canada-France-Hawaii Telescope (CFHT), carried out as a monitoring program over 4 consecutive observing seasons (2019, 2020, 2021 and 2022).After describing the observational material we collected (in Sec. 2) and briefly summarizing the latest estimates for the main stellar atmospheric parameters (in Sec.3), we present our measurements of the large-scale magnetic field of TW Hya (in Sec. 4) and the magnetic modeling that we derive from our spectropolarimetric observations using tomographic imaging techniques (in Sec. 5).We then outline our radial velocity (RV) measurements of TW Hya and their modeling in a Bayesian framework (in Sec.6), and describe the characteristics of the nIR emission lines tradition-1 All abbreviations used in the paper are listed in Sec.A for easier reference.
ally probing accretion in cTTSs, as well as their temporal behaviour (in Sec. 7).We finally conclude our study and discuss its implications for our understanding of star / planet formation in magnetized cTTSs (in Sec. 8).

SPIROU OBSERVATIONS
TW Hya was observed over several successive seasons with the SPIRou nIR spectropolarimeter / high-precision velocimeter (Donati et al. 2020b) at CFHT, first within the SPIRou Legacy Survey (SLS) in 2019, 2020 and 2021, then within the PI program of Lisa Lehmann in 2022 (run IDs 22AF14 and 22AF96).SPIRou collects unpolarized and polarized stellar spectra, covering a wavelength interval of 0.95-2.50m at a resolving power of 70 000 in a single exposure.For the present study, we concentrated on circularly polarized (Stokes ) and unpolarized (Stokes ) spectra of TW Hya only.Each polarization observation consists of a sequence of 4 subexposures, associated with different azimuths of the Fresnel rhomb retarders in order to remove systematics in polarization spectra (to first order, see, e.g., Donati et al. 1997).Each sequence yields one Stokes and one Stokes spectrum, as well as a null polarization check (called ) allowing one to diagnose potential instrumental or data reduction issues.
A total of 84 polarization sequences of TW Hya were collected in 4 main seasons, 11 in 2019 (April), 14 in 2020 (February to May), 30 in 2021 (February to May), and 29 in 2022 (March to May).A single polarization sequence was recorded in most nights; however, in a few cases (on 2021 March 28, April 26, April 28 and 2022 March 23), a second sequence was collected when data quality was lower than usual in the first one.Two spectra were discarded due to very low signal to noise ratios (SNRs), one in 2021 (April 28) and another one in 2022 (May 18).It finally yielded a total of 82 usable Stokes , and spectra of TW Hya, with 11, 14, 29 and 28 of them in 2019, 2020, 2021 and 2022 respectively, spanning in each case time slots of 12, 100, 70 and 70 d, and altogether covering a temporal window of 1131 d.The full log of our observations is provided in Table B1 of Appendix B.
Our SPIRou spectra were processed with Libre ESpRIT, the nominal reduction pipeline of ESPaDOnS at CFHT, optimized for spectropolarimetry and 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 =4000 K and a logarithmic surface gravity log =4.0 adapted to TW Hya (see Sec 3).Atomic lines of relative depth larger than 10 percent were selected, for a total of ≃1300 lines, featuring an average wavelength and Landé factor of 1750 nm and 1.2 respectively.The noise levels in the resulting Stokes LSD profiles range from 1.1 to 3.8 (median 1.6, in units of 10 −4 where denotes the continuum intensity).We also applied LSD with a mask containing the CO lines of the CO bandhead (at 2.3 m) only, to obtain veiling estimates in the K band, in addition to those for the whole spectrum derived from LSD profiles of atomic lines (see Sec. 3).Phases and rotation cycles were derived assuming a rotation period of rot = 3.606 d (see Sec. 4) and counting from an arbitrary starting BJD0 of 2458488.5 (i.e., just prior to our first SPIRou observation).
Our data were also processed with APERO (version 0.7.288), the nominal SPIRou reduction pipeline (Cook et al. 2022) optimized for RV precision.The reduced spectra were first analyzed for Zeeman broadening (see Sec. 4), then for RVs (see Sec. 6) with the line-by-line (LBL) technique (version 0.63, Artigau et al. 2022).It yielded 79 nightly RVs (corrected from instrumental drifts monitored with the SPIRou RV reference module, Donati et al. 2020b) and associated error bars (median 2.2 m s −1 ), listed in Table B1.

FUNDAMENTAL PARAMETERS OF TW HYA
In this section we recall the main parameters of TW Hya from the literature, including our previous study from ESPaDOnS spectra (Donati et al. 2011).TW Hya is a cTTS located at a distance of 59.96 +0.37  −0.11 pc from the Sun (Gaia Collaboration et al. 2023), in the TWA association, and hosts the protoplanetary disc closest to the Solar System.According to Gaia again, TW Hya features a photospheric temperature of eff = 3850 ± 10 K, a logarithmic surface gravity of log = 4.05 ± 0.01 dex (cgs units) and a metallicity relative to the Sun of [M/H] = −0.50 ± 0.05 dex.We suspect that, despite the small error bars, these estimates, and in particular [M/H], are likely affected by surface spots induced by magnetic activity and by intrinsic variability caused by accretion from the surrounding disc (e.g., Herczeg et al. 2023).Although Herczeg & Hillenbrand (2014) also suggest, based on low-resolution spectra, that eff ≃ 3810 K (corresponding to a M0.5 spectral type) looks adequate, several studies involving highresolution spectra rather conclude that eff is higher and more consistent with a K7 spectral type (e.g., Torres et al. 2003;Yang et al. 2005).
Applying our own spectral characterization tool on selected atomic lines in the ESPaDOnS spectra of TW Hya from our previous study yields eff = 4060 ± 50 K and log = 4.2 ± 0.1 dex assuming solar metallicity, as appropriate for nearby young stars (e.g., Padgett 1996).When applying ZeeTurbo (a tool specifically developed for SPIRou spectra of M dwarfs, Cristofari et al. 2022aCristofari et al. ,b, 2023a,b) ,b) to our SPIRou spectra of TW Hya, and assuming again solar metallicity, we obtain eff = 3970 ± 50 K, including the effect of magnetic fields (see Sec. 4) but fixing log to 4.2 and the veiling in the YJH and K bands to 0.20 and 0.25 respectively (as measured, see below) to minimize correlation between fitted parameters.This is slightly cooler than (though still consistent with) our estimate from ESPaDOnS spectra, likely as a result of cool magnetic spots at the surface of TW Hya (Huélamo et al. 2008;Donati et al. 2011, see also the following sections) that affect nIR lines more than optical ones (the lower brightness contrast between cool spots and the photosphere in the nIR implying a larger relative contribution of these spots to nIR lines).This is why we did not include the CO lines of the CO bandhead in the fit, as these lines get stronger with decreasing temperature and thereby further bias the determination of eff towards spot temperatures.Our measurements are consistent with those of another study from high-resolution nIR spectra of TW Hya (Sokal et al. 2018, also taking into account magnetic fields and assuming solar metallicity), yielding eff = 3800 ± 100 K and log = 4.2 ± 0.1 dex.The agreement is best for log and less so for eff , presumably for the reason already outlined above, i.e., the presence of cool star spots affecting temperature determination (see also Gully-Santiago et al. 2017).
In the following, we use the estimates derived from our optical data, less affected by cool surface spots.Comparing with log vs eff synthetic tracks from the evolutionary models of Baraffe et al. (2015) yields for TW Hya a mass of ★ = 0.80 ± 0.05 M ⊙ (in good agreement with the dynamical mass derived from ALMA observations of the almost pole-on accretion disc, equal to 0.81 ± 0.16 M ⊙ , Teague et al. 2019), an age of 7.5 ± 2.5 Myr, a radius of ★ = 1.16 ± 0.13 R ⊙ consistent with interferometric measurements (1.29 ± 0.19 R ⊙ , Gravity Collaboration et al. 2020), and a logarithmic luminosity with respect to the Sun of log( ★ /L ⊙ ) = −0.48± 0.10.These evolution models also predict that TW Hya already started to develop a radiative core of mass 0.2 ± 0.1 M ⊙ .
The rotation period of TW Hya was unambiguously determined to be about 3.6 d from the reported periodic changes in RV, that were first erroneously attributed to the reflex motion of a putative massive close-in planet (Setiawan et al. 2008) then to magnetic activity inducing rotational modulation (Huélamo et al. 2008;Donati et al. 2011;Sicilia-Aguilar et al. 2023).Our new spectropolarimetric data confirm this value, with the line-of-sight component of the large-scale magnetic field integrated over the visible stellar hemisphere, called longitudinal magnetic field and denoted ℓ , steadily varying with a period of 3.606 ± 0.010 d throughout the 1131 d of our observations (see Sec. 4).The corresponding corotation radius, at which the Keplerian rotation rate equals that at the surface of the star, is equal to cor = 0.043±0.001au (7.9±0.3 ★ ).We note that photometric observations of TW Hya (e.g., Rucinski et al. 2008;Siwak et al. 2014Siwak et al. , 2018)), including the most recent ones collected with TESS in March 2019, 2021 and 2023 (each lasting about 25 d, and the first two contemporaneous with our SPIRou spectra), rarely exhibit rotational modulation, but rather a spectrum of unstable periods of order a few days likely probing intrinsic variability triggered by unsteady accretion from the inner regions of the accretion disc.For instance, Fig. C1 shows stacked periodograms of the March 2019 and 2021 binned TESS light curve, that exhibit transient periodic signals, with only little power at the rotation period (see also Sicilia-Aguilar et al. 2023).
Given this rotation period and the stellar radius derived above (1.16 ± 0.13 R ⊙ ), we can conclude that the line-of-sight projected equatorial rotation velocity sin of TW Hya must be small if the rotation axis of the star is close to the line of sight, as is that of the accretion disc ( = 5.8 +4.0 −1.7 for the disc, Teague et al. 2019).For instance, sin = 5.8 km s −1 (Sokal et al. 2018) or sin = 8.4 km s −1 (López-Valdivia et al. 2021) would imply inclination angles of the stellar axis to the line of sight of ≃ 20 • and 30 • respectively, i.e., significantly larger than that of the accretion disc.We adopt here a value of sin = 3 ± 1 km s −1 corresponding to an inclination of the stellar rotation axis ≃ 10 • , slightly larger than, though still consistent with, that of the outer disc.In addition to yielding a better match to the observed Stokes profiles than the larger estimates mentioned above (see Sec. 5), this sin is more consistent with our previous measurement (i.e., sin = 4 ± 1 km s −1 , Donati et al. 2011), with that of other independent studies (e.g., Sicilia-Aguilar et al. 2023), and with the inclination of the outer disc.
As mentioned above, the accretion disc of TW Hya, extending to a few hundred au's, features dusty rings and gaps (including one at 1 au from the star, Andrews et al. 2016;Nomura et al. 2016;van Boekel et al. 2017;Nomura et al. 2021) as well as spiral structures (Teague et al. 2019), possibly tracing the presence of planets forming and migrating throughout the disc.The detection of moving shadows at the surface of the outer disc (Debes et al. 2017) also suggests that the inner disc regions may not be coplanar with the outer ones (Teague et al. 2022;Debes et al. 2023), hence potentially supporting that the rotation axis of the central star is slightly misaligned with that of the disc, and / or that planets are indeed present in the innermost regions of the disc.Although low, mass accretion is still observed at the surface of the star (at an average To double check veiling estimates, we used our LSD Stokes profiles of TW Hya, both for atomic lines over the whole spectrum and for the CO lines of the CO bandhead in the K band, and compared them with the corresponding median LSD Stokes profiles of the weak-line TTS TWA 25, also observed with SPIRou and that we take as a reference for an unveiled spectrum of similar spectral type (as in Sousa et al. 2023).We find that the median veiling is about 20% over the whole spectrum, and only marginally larger (25%) in the K band (CO bandhead), slightly stronger than (though still consistent with) the estimates of Sousa et al. (2023) and in agreement with measurements derived from older nIR spectra of TW Hya (Lavail et al. 2019).Besides, we obtain that the median veiling per season is about the same (within a few %) for all four seasons.We also confirm that veiling in individual spectra is about 2.5× more variable with time in the K band (about 15% rms about median veiling) than in the rest of the spectrum, again qualitatively consistent with Sousa et al. (2023).
Table 1 summarizes the main parameters of TW Hya used in, or derived from, the present study.

THE LONGITUDINAL FIELD AND ZEEMAN BROADENING OF TW HYA
The next step in our analysis is to derive the longitudinal field ℓ of TW Hya following Donati et al. (1997), from each of the Stokes and LSD profiles derived in Sec. 2. In practice, we computed the first moment of the Stokes profile and its error bar, whereas the equivalent width of the Stokes LSD profiles is measured through a Gaussian fit.Stokes LSD signatures were integrated over a window of ±40 km s −1 in the stellar rest frame, given the strong magnetic broadening of line profiles (Yang et al. 2005; Sokal et al. ), with the exact integration width having little impact on the result.We proceeded in the same way with the polarization check to verify that the derived pseudo longitudinal field is consistent with 0 within the error bars, i.e., associated with a reduced chi-square 2 r close to unity.The inferred ℓ values computed from our 82 Stokes profiles are listed in Table B1 and range from −195 to 77 G (median −37 G) with error bars of 10 to 34 G (median 15 G), yielding a 2 r (with respect to the ℓ = 0 G line) equal to 32.3 for (1.09 for ).This demonstrates that the magnetic field is unambiguously detected in the Stokes signatures of TW Hya, that no spurious pollution is observed in and that our analytical error bars are consistent with the observed dispersion within 5%.Unsurprisingly, our ℓ values are an order of magnitude weaker than the small-scale fields estimated from the Zeeman broadening of nIR lines (≃3 kG, Yang et al. 2005;Sokal et al. 2018;Lavail et al. 2019;López-Valdivia et al. 2021, see also below), as usual for cool active stars harboring small-scale tangled fields whose circular polarization signatures mostly cancel out.
We then investigate the temporal behaviour of our ℓ data, arranged in a vector denoted y, using quasi-periodic (QP) Gaussian-Process Regression (GPR), with a covariance function ( , ′ ) of 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 ℓ modulation changes, and 4 a smoothing parameter describing the amount of harmonic complexity needed to describe the data.We then select the QP GPR fit that features the highest likelihood L, defined by: where C is the covariance matrix for our 82 epochs, the diagonal variance matrix associated with y, and S = 2 5 J (J being the identity matrix) the contribution from an additional white noise source that we introduce as a fifth hyper-parameter 5 (in case our error bars on ℓ were underestimated for some reason).The hyperparameter 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 whereas the derived hyper-parameters are listed in Table 2.The first surprising conclusion is that ℓ is changing sign from epoch to epoch, being mostly positive in 2019 and 2022 but negative in 2020 and 2021.As we will see in Sec. 5, this does not reflect an overall polarity switch of the magnetic field, but rather results from small changes of the large-scale topology and / or of the surface brightness distribution in an almost pole-on viewing configuration.This is fairly different from the results of our earlier optical observations with ESPaDOnS, where the longitudinal field in photospheric lines was always positive, reaching several hundred G, and that from accretion lines (Ca infrared triplet and He 3 lines in particular, Donati et al. 2011) was always negative.What we detect with SPIRou is in between, reflecting that magnetic regions are spatially associated with dark surface features, with a spot-to-photosphere brightness contrast that is smaller in the infrared than in the optical (see Sec. 5) as was the case for CI Tau (Donati et al. 2024).
Our second main result is that we are able to detect rotational modulation of ℓ , especially in 2020, despite the amplitude of the modulation being small as a result of the close to pole-on viewing configuration.The rotation period we measure is equal to rot = 3.606 ± 0.015 d, slightly larger than that derived from optical RV data (i.e., 3.5683 ± 0.0002 d, Huélamo et al. 2008).This suggests that weak differential rotation is present at the surface of TW Hya (at a level of only a few mrad d −1 ), as further discussed in Sec. 6.We also find that the rotational modulation of the ℓ curve is simple enough for the GPR fit to yield a smoothing parameter that is large and weakly constrained by the data, which we thus fix at its optimal value ( 4 = 2.5).The evolution timescale ( 3 = 232 +64 −50 d) is about 3× longer than that of the more evolved young active star AU Mic, whose ℓ curve is also more complex (Donati et al. 2023).We finally outline that ℓ can almost be fitted down to the noise level ( 2 r = 1.13) and thus that the additional white noise term 5 is only slightly larger than zero.
We also carried out a Principal Component Analysis (PCA) of our Stokes profiles (as outlined in Lehmann & Donati 2022;Lehmann et al. 2024) and find that the first PCA eigenvector is capable of reproducing most of the observed variations of the meansubtracted Stokes profiles (see Fig. D1).The second PCA eigen- vector, encoding wavelength shifts of the Stokes profiles (given its shape that mimics the derivative of the first eigenvector), is only marginally required, indicating that the parent magnetic regions do not travel much throughout the line profile, i.e., are located at high latitudes.We also find that the mean Stokes profile is antisymmetric with respect to the line center and dominates over the meansubtracted profiles, i.e., that the large-scale field is mainly poloidal and axisymmetric, as expected from the nearly pole-on configuration of TW Hya that renders us almost insensitive to axisymmetric toroidal fields at the surface of the star (nearly perpendicular to the line of sight).The 1 coefficient associated with the first PCA eigenvector (scaled and shifted to match ℓ , as in Lehmann et al. 2024) exhibits a time dependence very similar to that of ℓ (see Fig. 2), with little to no rotational modulation depending on the season (see Fig. D1), typical of a simple poloidal field only slightly tilted to the rotation axis.These preliminary conclusions are confirmed with the full magnetic modeling of TW Hya presented in Sec. 5.  2023a), which incorporates magnetic fields as well as a MCMC process to determine optimal parameters and their error bars.We find that TW Hya hosts a small-scale magnetic field of < > = 3.60±0.04kG, with the relative area of non-magnetic regions being 0 = 7 ± 2%.
Looking at the Zeeman broadening of atomic and molecular lines with ZeeTurbo (Cristofari et al. 2023a,b) applied to our median spectrum of TW Hya, we find that 4 components, associated with small-scale magnetic fields of strengths 0, 2, 4 and 6 kG, and respective filling factors 0 = 7 ± 2%, 2 = 40 ± 3%, 4 = 20 ± 3% and 6 = 33 ± 2% of the visible stellar surface, are needed to obtain a good fit, yielding an overall small-scale field measurement of < > = 3.60 ± 0.04 kG (see Fig. 3).Our estimate of the smallscale field of TW Hya from SPIRou spectra is consistent with those derived in previous studies (Yang et al. 2005;Sokal et al. 2018;Lavail et al. 2019;López-Valdivia et al. 2021), given the expected temporal variability.On the timescale of our observations, we find only marginal variations of the small-scale field (typically 0.1 kG rms on < >) between our 4 observing seasons, and detect no rotational modulation of < > on measurements from individual spectra (which is unsurprising given the nearly pole-on viewing configuration of TW Hya).

ZEEMAN-DOPPLER IMAGING OF TW HYA
In this section, we analyse the Stokes and LSD signatures of TW Hya from each season using Zeeman-Doppler Imaging (ZDI), in order to simultaneously reconstruct the topology of the largescale field and the associated brightness map, as well as their temporal evolution over the four seasons of our observations.We achieve this through an iterative process that progressively adds information at the surface of the star, starting from a small magnetic seed and a featureless brightness map and exploring the parameter space with a variant of the conjugate gradient technique that aims at efficiently minimizing the discrepancy between the observed and synthetic Stokes and LSD profiles (e.g., Skilling & Bryan 1984;Brown et al. 1991;Donati & Brown 1997;Donati et al. 2006;Kochukhov 2016).Since the problem is ill posed, regularization is needed to ensure a unique solution.ZDI uses the principles of maximum entropy image reconstruction, which aims at reaching a given agreement with the data, usually2 r ≃ 1, while minimizing information in the derived maps to ensure that reconstructed features are mandatory to reproduce the data.
In practice, we describe the surface of the star as a grid of 5000 cells and compute synthetic Stokes and profiles at each observation epoch by summing up the spectral contributions of all grid cells, taking into account the main geometrical parameters such as the cell coordinates, ≃ 10 • , sin = 3 km s −1 , and the linear limb darkening coefficient (set to 0.3).We also assume that the surface of TW Hya rotates as a solid body over each season, consistent with the low level of differential rotation observed on TW Hya (see Secs. 4 and 6).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 LSD profiles, and a Doppler width D = 3 km s −1 (including thermal, micro and macrotubulent broadening) for the local profile (as for CI Tau, see, e.g., Donati et al. 2024).
The relative brightness at the surface of the star is described as a series of independent pixels, whereas magnetic field is expressed 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 depend on 3 sets of complex SH coefficients, ℓ, and ℓ, for the poloidal component, and ℓ, for the toroidal component 2 , where ℓ and note the degree and order of the corresponding SH term in the expansion.Given the low sin of TW Hya, we can safely limit the expansion to terms up to ℓ = 5.As for other cTTSs magnetically imaged with ZDI to date, we favour large-scale field configurations that are mostly antisymmetric with respect to the centre of the star, in which accretion funnels linking the inner disc to the star are anchored at high latitudes, which is achieved in practice by penalizing even SH modes with respect to odd ones in the entropy function (as in, e.g., Donati et al. 2011).
Finally, we assume that only a fraction of each grid cell (called filling factor of the large-scale field, equal for all cells) contributes to Stokes profiles, with a magnetic flux over the cell equal to (i.e., 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 large-scale field (with a scaling factor of / ), which ensures at least that the resulting Zeeman broadening from small-scale fields is consistent with the reconstructed large-scale field.As for the cTTS CI Tau (Donati et al. 2024), we set ≃ 0.8 (consistent with our results, see Sec. 4, and with those of Yang et al. 2005) and ≃ 0.4, which yields a satisfactory fit to the observed Stokes and profiles of TW Hya, and reproduces in particular the conspicuous triangular shape of the magnetically broadened Stokes LSD profiles.
Fitting our LSD Stokes and profiles with ZDI (down to scale field, shown in Fig. 5.Although brightness maps (not shown) were reconstructed at the same time as magnetic maps (with ZDI simultaneously fitting Stokes and profiles), we find that no brightness feature at the surface of TW Hya is large enough, or exhibits a big enough nIR contrast with respect to the surrounding photosphere, to generate clear Stokes profile distortions and the corresponding rotational modulation, and thereby to show up in the derived brightness maps.This is in contrast with the brightness images we had derived from optical data showing an obvious dark feature coinciding with the magnetic pole (Donati et al. 2011), but similar to our findings on CI Tau where Stokes profile distortions induced by surface brightness features in the nIR were also barely detectable among the dominant ones induced by magnetic fields (conversely to the optical domain where the opposite behaviour holds).
We find that TW Hya hosts a large-scale magnetic field of average strength ≃1.1 kG over the star and reaches a maximum intensity of 1.5-2.0kG, which translates into average and maximum small-scale fields of 2.2 and 3-4 kG respectively (taking into account the / ≃ 2 ratio between both quantities), consistent with literature values (e.g., Yang et al. 2005;Sokal et al. 2018;López-Valdivia et al. 2021).The large-scale field we reconstruct is almost fully poloidal and axisymmetric, and mainly consists of a 1.0-1.2kG dipole inclined at 20 • with respect to the rotation axis.The octupole component is significantly weaker, with a polar strength ranging from 0.2 to 0.5 kG in the different seasons, and adds up to the polar large-scale field values, generating at times local maxima aside the main one in the radial field map, like in the 2021 and 2022 seasons.The main properties of the reconstructed magnetic topologies, consistent with the preliminary conclusions of our PCA analysis (see Sec. 4), are summarized in Table 3.
We can see in particular that the large-scale magnetic topology is not undergoing global polarity switches over our 4-season timespan, despite the longitudinal field changing sign from 2019 to 2020 and again from 2021 to 2022 (see Fig. 1 and Table B1).It shows that, in a nearly pole-on viewing configuration like that of TW Hya, sign switches in the longitudinal field may also reflect changes in the relative contributions of the radial field near the pole and the meridional field at lower latitudes (see Fig. 5), with the first dominating over the second when the strength of the dipole field is largest (i.e., in 2020 and 2021, see Table 3).
We note that the magnetic maps we derive from our SPIRou data differ from those reconstructed from ESPaDOnS data collected a decade earlier (Donati et al. 2011), most likely as a result of changes in the large-scale field topology and accretion pattern, that we know can occur on relatively short timescales on TW Hya (Herczeg et al. 2023).In particular, the octupole component measured from our SPIRou data is much smaller (by a factor of 5-10) than that derived from previous ESPaDOnS observations, whereas the dipole component is about twice larger (Donati et al. 2011).We Table 3. Properties of the large-scale and small-scale magnetic field of TW Hya for our 4 observing seasons.We list the average reconstructed large-scale field strength < > (column 2), the maximum small-scale field strength (column 3), the polar strength of the dipole component d (column 4), the tilt / phase of the dipole field to the rotation axis (column 5) and the amount of magnetic energy reconstructed in the poloidal component of the field and in the axisymmetric modes of this component (column 6).Error bars are typically equal to 5-10% for field strengths and percentages, and 5-10 • for field inclinations.estimate that most of this evolution is real, although we cannot exclude that some of it relates to the difference in the data sets and in particular in the wavelength domains.Future analyses simultaneously combining optical and nIR spectropolarimetric data, such as that recently carried out for CI Tau (Donati et al. 2024) should allow one to address this point in a more extensive way.Despite such changes, the dipole tilts we measure from our new maps, of order 20 • (see Table 3), are consistent with the off-centring of the mainly polar accretion region taking place at the surface of TW Hya, as derived from previous studies including ours (Donati et al. 2011;Sicilia-Aguilar et al. 2023).

RADIAL VELOCITY MODELING OF TW HYA
Using now the spectra of TW Hya reduced with APERO (Cook et al. 2022), we can derive precise RVs with LBL (Artigau et al. 2022), listed in Table B1, with a median RV precision of 2.2 m s −1 that reflects the relatively sharp lines of this star, and in particular the magnetically insensitive (unbroadened) molecular features.We find that TW Hya is RV stable at a rms level of 32.5 m s −1 over our 4 observing seasons, a dispersion about 15× larger than the median error bar on individual RV points.By carrying out a GPR fit to these RVs, we find that 80% of these variations are caused by rotational modulation, with a period equal to 3.5647 ± 0.0024 d, slightly but definitely smaller than the period derived from our ℓ data (see Sec. 4).This difference argues again for the presence of small latitudinal differential rotation at the surface of TW Hya (at a level of a few mrad d −1 ), with the polar regions contributing most to the ℓ variations (see Fig. 5) rotating more slowly than lower latitudes generating most of the RV modulation.The semi-amplitude of this RV modulation is 25.5 +5.6  −4.6 m s −1 , whereas the additional white noise on RVs (presumably caused by accretion-induced intrinsic distortions of spectral lines) reaches 19.0 ± 1.7 m s −1 , i.e., 8.6× the median error bar of our RV measurements (see Table 4).Note that two of the GPR hyper-parameters ( 3 and 4 ), weakly constrained by the data, were fixed to their optimal value from a preliminary run with all GPR parameters free to vary, making no difference on the filtering of activity.
Once activity is filtered out, we find residual power at a period of about 8.3 d in the RV periodogram.We thus ran a new series of GPR fits to our RV points through a Bayesian Monte Carlo Markov chain experiment, including now the presence of a putative close-in planet in circular orbit around TW Hya at a period of about 8.3 d described with 3 additional parameters.We find that a relatively clear RV signal is present in the data at a 4.3 level, with a semi-amplitude of 11.1 +3.3  −2.6 m s −1 and a period of 8.339 ± 0.008 d, corresponding to an orbital distance of 0.075 ± 0.001 au.The RV residuals are now smaller (15.6 m s −1 instead of 17.1 m s −1 rms) than in the reference (no planet) case, and so is the additional white noise (17.1 m s −1 instead of 19.0 m s −1 ).The corresponding increase in marginal log likelihood log BF = Δ log L reaches 8.7 (see Table 4), suggesting that the detected RV signal is real (i.e., log BF > 5, Jeffreys 1961).The corresponding fit to our RV points is shown in Fig. 6, along with the filtered RVs, the fitted RV signal and the RV residuals.The 1-yr alias of the reported RV signal, located at a period of 8.147 ± 0.010 d, is also a potential solution (which we refer to as case b' in Table 4), albeit with a slightly lower confidence level (log BF = 8.1).The corresponding periodograms are depicted in Fig. E1, along with the stacked periodogram illustrating how the main RV signal and its 1-yr alias get stronger and more dominant as more spectra are included in the analysis (see Fig. E2).The associated corner plot is shown in Fig. E3.We note that the period of this RV signal is slightly smaller than that of the photometric one that dominates the March 2021 and 2023 TESS light curves of TW Hya (at about 9 d, see, e.g., Fig. C1), though it is not clear whether both are related.
If this RV signal is generated by a true planet, the minimum mass of this orbiting body would be 30 +9 −7 M ⊕ (or 28 +9 −6 M ⊕ in case b'); further assuming that the planet orbital plane coincides with the equatorial plane of the star yields a planet mass of = 0.55 +0.17 −0.13 M (0.51 +0.17 −0.12 M in case b').The corresponding phase-folded RV curve is shown in Fig. 7 for our main solution (with case b' yielding a very similar plot).Running again the same experiment assuming now a more general Keplerian orbit, we derive an eccentricity consistent with zero (error bar 0.04), and no obvious improvement in log BF with respect to the circular orbit case.
Due to the presence of a disc around TW Hya, one may wonder whether the RV signal we detect truly comes from the reflex motion of the host star under the gravitational pull of a planet.One could for instance suspect the disc itself to contribute to the spectral lines of TW Hya, in particular the molecular lines that dominate the spectrum, and generate a modulated RV signal that would rather reflect, e.g., a non-axisymmetric structure within the disc rather than a genuine planet.If this were the case, one would expect the molecular lines to be more affected than atomic lines, or lines in the and bands to be less impacted than those in the band, as a result of the different temperatures of the stellar photosphere and disc material.We thus also analysed RVs obtained from Gaussian fits to Stokes LSD profiles of atomic lines only, that presumably come from the stellar photosphere with no contribution from the disc.If the detected RV signal were not present in atomic lines (as for CI Tau, e.g., Donati et al. 2024), this would argue for it being induced by the disc.However, atomic lines are 3× broader (as a consequence of magnetic fields) than molecular lines in TW Hya, and telluric residuals affect LSD profiles more than LBL RV measurements.As a result, the RV precision we obtain from atomic lines is significantly worse than in our main analysis, with an excess white noise from the GPR fit reaching 60 m s −1 (instead of 17 m s −1 when using LBL RVs from all spectral lines, see Table 4).Using LSD profiles from CO lines yields better results than from atomic lines (with an excess noise reduced by a factor of 2, down to 30 m s −1 ), but still not good enough to unambiguously detect the RV signal detected from LBL measurements.We also looked at the LBL RV measurements using lines from the , and bands only, and again find that the excess noise (equal to 50, 32 and 25 m s −1 for the , and bands respectively) is still too large to enable a Table 4. MCMC results for the 3 cases (no planet, planet b and planet b') of our RV analysis of TW Hya.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 model without planet.Two GPR hyper-parameters, weakly constrained by the data, were fixed to their optimal value from a preliminary run with all GPR parameters free to vary.

Parameter
No planet b b' Prior 1 (m s −1 ) 25.5  firm detection of the RV signal in each band (with respective error bars of 10, 6 and 5 m s −1 on the semi-amplitude) and thus to look for potential differences between them.It is therefore not possible at this stage to either confirm nor refute the planetary origin of the RV signal we report here.

EMISSION LINES OF TW HYA
In this penultimate section, we discuss the main emission lines present in the nIR spectra of TW Hya, and in particular the 1083.3nmHe triplet, as well as the 1282.16-nmPa and 2166.12-nmBr lines, known to probe accretion flows as well as outflows, in particular for the He triplet with its conspicuous P Cygni profile featuring a broad and strong blue-shifted absorption component.
The stacked profiles and the associated 2D periodograms over the full data set are shown in Fig. 8 for He and Pa , whereas those of Br are depicted in Fig. F1.We note that the He blue-shifted absorption, whose shape suggests it is formed within the stellar wind rather than from a disc wind (Kwan et al. 2007), is strongly variable with time, sometimes extending blue-wards as far as −300 km s −1 but only down to −150 km s −1 at other epochs.Its median equivalent width (EW) is about 100 km s −1 (0.36 nm, with no scaling from veiling).With a median EW of about 150 km s −1 (0.54 nm), the emission component is also quite variable, sometimes dominating the whole profile and at other times almost non-existent.About half the spectra show red-shifted absorption at velocities of 200 km s −1 , likely tracing accreted material from the disc falling onto the polar regions of TW Hya.
Neither the blue-shifted nor the red-shifted absorption components are modulated with rotation, even in individual seasons.This may sound surprising at first glance, at least for the red-shifted absorption given the conclusion of Sec. 5 that accretion occurs mostly towards the pole on TW Hya in a more or less geometrically stable fashion.However, TW Hya being viewed almost pole-on, accreted material only comes in front of the stellar disc once close to the surface of the star where it ends up being visible all the time, thereby rendering rotational modulation much smaller than intrinsic variability.This is likely the same for the blue-shifted component, especially if formed within a stellar wind; alternatively, it may suggest that this component traces a wind from the inner disc (rather than from the star) for which no rotational modulation is expected, and no more than marginal evidence for longer periods is observed (apart from those attributable to the window function at the synodic period of the Moon, i.e., 29.5 d, and its 1-yr aliases, see Fig. E1).Such incoherent variability of emission lines is similar to what is seen in photometric light curves (e.g., Rucinski et al. 2008;Siwak et al. 2014Siwak et al. , 2018;;Sicilia-Aguilar et al. 2023), including the 2019 and 2021 TESS light curves collected at the same time of our SPIRou observations (see Fig. C1 for their stacked periodograms).
We also note that a clear Stokes Zeeman signature is visible in the weighted average of all He profiles, as well as in those of each individual season.These signatures centred in the stellar rest frame and falling in conjunction with the emission peak (see Fig. 9), demonstrate that at least part of the He emission comes from the footpoints of accretion funnels, i.e., close to where the large-scale field is strongest, and indicate the presence of an axisymmetric magnetic field component of negative polarity that is visible at all times.This agrees well with our reconstructed ZDI maps that indeed show a negative radial field region close to the pole (see Fig. 5) and thereby always visible to the observer given the viewing angle of TW Hya.Assuming that the longitudinal field over the accretion region is similar to that previously probed by the 588 nm He 3 line (i.e., 2.5-3.0 kG, see Donati et al. 2011), it implies that about 20% of the emission flux in the 1083.3-nmHe triplet, i.e., a component of EW ≃30 km s −1 (0.11 nm), is coming from the post-shock accretion region within the chromosphere of TW Hya.
The Pa line of TW Hya shows a simpler profile, with a main emission peak that is slightly blue-shifted (by -2.7 km s −1 in average) and features an extended blue wing.The EW of Pa , measured through a simple Gaussian fit without any scaling for veiling, are listed in Table B1 and vary from 120 to 616 km s −1 (0.51 to 2.63 nm), with a median of 300 km s −1 (1.28 nm).When scaled up for veiling, these EWs translate into logarithmic luminosities (relative to L ⊙ ) in Pa of −3.79 ± 0.18 (with the error bar corresponding to temporal variability), and thus into logarithmic accretion luminosities (again relative to L ⊙ ) of −1.25 ± 0.19 (following the calibrated relations of Alcalá et al. 2017).This implies an average logarithmic mass accretion rate of −8.48 ± 0.19 (in units of M ⊙ yr −1 ), with temporal variations in the range −8.88 to −8.13 (i.e., by a factor of 5.6 peak to peak).
The conspicuous absorption component that shows up in the red wing at a velocity of 84 km s −1 is caused by a photospheric line from Ti (located at 1282.52 nm).The red-shifted absorption visible in the He line (tracing accreting material from the disc about to reach the stellar surface) is apparently also present in Pa , though much shallower.The 2D periodogram of Pa over our full data set shows no clear feature apart from those attributable to the window function (see Fig. E1) and already present in the He periodogram.By running GPR through the EWs of Pa , we further confirm that no clear period emerges from the noise, dominated by accretioninduced intrinsic variability.As for He , we speculate that the nondetection of rotational modulation is due to the viewing angle under which TW Hya is seen from the Earth.We also detect a clear Zeeman signature in conjunction with Pa (see Fig. 9), albeit with a lower amplitude than that in He , and again probing the presence of an axisymmetric magnetic field component of negative polarity at the surface of TW Hya.Assuming this axisymmetric magnetic component is the same as that detected in He , it implies that the corresponding post-shock region in the chromosphere of TW Hya contributes to the emission of Pa at an average EW of ≃20 km s −1 (0.09 nm), i.e., about 7% that of the whole Pa emission.
Similar results are derived from Br (see Fig. F1), with a comparable overall blue-shift (of -2.3 km s −1 ) and EWs (listed in Table B1) that vary from 19 to 127 km s −1 (0.14 to 0.92 nm) with a median of 57 km s −1 (0.41 nm).These EWs translate into logarithmic luminosities in Br of −4.83 ± 0.21 and into logarithmic accretion luminosities of −1.73 ± 0.24 (both relative to L ⊙ ), implying an average logarithmic mass accretion rate of −8.96 ± 0.24 (in units of M ⊙ yr −1 , with temporal variations in the range −9.51 to −8.53, i.e., by a factor of 9.5 peak to peak).A Zeeman signature is again detected in Br with the same characteristics as that of Pa ; assuming once more that it probes the same axisymmetric magnetic component (of negative polarity), we can infer that the post-shock (dashed horizontal lines), apart from peaks close to the synodic period of the Moon (29.5 d) and its 1-yr aliases (also present in the window function, see Fig. E1).
region at the footpoint of accretion funnels contributes to the emission of Br at an average EW of 7 km s −1 (0.05 nm), i.e., about 12% that of the whole Br emission.
The average logarithmic mass-accretion rate that we derive for TW Hya from our 2019 to 2022 SPIRou observations (using both Pa and Br ) is thus equal to −8.72±0.22(in units of M ⊙ yr −1 , with temporal variations in the range −9.19 to −8.34) in good agreement with the results of Herczeg et al. (2023) derived from 2.5 decades of irregular monitoring at optical wavelengths.
We finally note that no power is detected in either lines at the period of the RV signal reported in Sec.6, which is what we expect if the RV signal is not attributable to activity.However, as no power is detected at rot either, where one usually expects activity to show up, the non-detection at does not provide definitive evidence that this period is unrelated to activity.

SUMMARY AND DISCUSSION
We monitored the cTTS TW Hya with the SPIRou high-resolution nIR spectropolarimeter / velocimeter at CFHT over 4 consecutive seasons (from 2019 to 2022), in the framework of the SLS large program and of a PI program.We obtained a total of 82 usable Stokes and spectra of TW Hya on which the LSD and LBL methods were applied to derive Zeeman signatures and precise RVs for each of our observing nights.
The longitudinal field ℓ measured from Stokes and LSD profiles evolved with time, and even switched sign between 2019 and 2020, and again between 2021 and 2022.Rotational modulation of ℓ is also detected, yielding a period equal to rot = 3.606±0.015d, slightly but significantly larger than that derived from optical RVs collected in 2008 (i.e., 3.5683 ± 0.0002 d, Huélamo et al. 2008).We also detect rotational modulation of RVs in our nIR data, with a period of 3.5649 ± 0.0024 d, consistent at 1.4 with the estimate from optical RVs.It demonstrates that latitudinal differential rotation is present at the surface of TW Hya, with the polar regions (mostly probed by ℓ ) rotating more slowly than lower latitudes (to which RVs are mostly sensitive), but at a level of only a few mrad d −1 between the equator and pole, i.e., consistent with previous results on TTSs similar to TW Hya (e.g., Finociety et al. 2023).By modeling the Zeeman broadening of atomic and molecular lines of TW Hya, we also measured a small-scale field of 3.60 ± 0.04 kG (with 7 ± 2% of the visible stellar surface free of such field) and only small season-to-season variations, in rough agreement with previous literature estimates (Yang et al. 2005;Sokal et al. 2018;Lavail et al. 2019;López-Valdivia et al. 2021).
By carrying out a PCA analysis of our Stokes profiles (as advocated by Lehmann & Donati 2022;Lehmann et al. 2024), we find that the large-scale field of TW Hya is mostly poloidal and axisymmetric at all epochs.This conclusion is confirmed with a thorough modeling with ZDI, thanks to which we reconstructed the large-scale magnetic field of TW Hya, as well as the photospheric brightness at nIR wavelengths, for each observing season, from a simultaneous fit to the corresponding set of Stokes and LSD profiles.We find that the large-scale field is fairly stable with time, despite the sign switches in ℓ , with a dominant dipole component evolving from 1.0 kG (in 2019 and 2022) to 1.2 kG (in 2020) and 1.1 kG (in 2021).The sign switches that ℓ exhibits directly reflects this temporal evolution of the large-scale field, with the polar and lower latitude regions both contributing to ℓ through the radial and meridional field respectively, in the mostly pole-on viewing configuration of TW Hya.Besides, we find that the nIR brightness inhomogeneities at the surface of TW Hya only feature a low contrast with respect to the quiet photosphere, generating, along with the large-scale field, a rotational modulation of LBL RVs (from the narrower molecular lines mainly) whose semi-amplitude is only 25 ± 5 m s −1 , i.e., much lower than that from TTSs with high levels of spot coverage (e.g., Finociety et al. 2023).We also find that the small-scale fields derived from the shape and width of LSD profiles of atomic lines of TW Hya are consistent with typical values of the small-scale and large-scale filling factors of cTTSs, i.e., ≃ 0.8 and ≃ 0.4 (e.g., Donati et al. 2024), yielding average and maximum values of 2.2 and 3-4 kG respectively for the small-scale field at the surface of the star, again consistent with previous measurements including ours.
Comparing with our previous large-scale magnetic field maps of TW Hya from optical spectropolarimetric data collected a decade ago (Donati et al. 2011), we find a clear evolution, with the dipole component about twice stronger and the octupole component much weaker (by a factor of 5-10) than it used to be.We believe that most of this evolution is real, but cannot exclude that some of it reflects the difference in wavelength domains between both data sets.Similarly, we note that the semi-amplitude of the RV modulation, equal to 25 ± 5 m s −1 in our nIR data, is an order of magnitude smaller than that reported from optical RVs collected 1.5 decades ago (Huélamo et al. 2008;Donati et al. 2011), which again argues for intrinsic variability of the magnetic activity at the surface of TW Hya, the typical ratio between the nIR and optical RV jitter being usually much smaller than 10 (e.g., Mahmud et al. 2011;Crockett et al. 2012).These results strongly argue in favour of collecting spectropolarimetric and velocimetric data in both optical and nIR domains at the same time so that one can simultaneously use information from both spectral ranges, as recently done in the case of the cTTS CI Tau (Donati et al. 2024).This should be routinely possible in a few months once SPIRou and ESPaDOnS are merged into a single instrument (called VISION) allowing one to simultaneously observe the same star in both wavelength domains.
Given the reported mass accretion rate at the surface of TW Hya over the previous 2 decades (with log ≃ −8.65 M ⊙ yr −1 , ranging from -9.2 to -8.2, Herczeg et al. 2023, and consistent with the one we derive from our SPIRou observations, see Sec. 7), we find that the magnetic truncation radius mag (as defined by Bessolaz et al. 2008, and up to which the large-scale field is able to disrupt the Keplerian disc) is equal to mag = 4.5 +2.0 −1.1 ★ , i.e., 0.57 +0.25  −0.15 cor , the error bar reflecting mostly the reported variation in mass accretion rate rather than that in the dipole component of the large-scale field.Our result is consistent with the recent interferometric measurement of the magnetospheric radius of TW Hya, equal to 4.50 ± 0.26 ★ (Gravity Collaboration et al. 2020) at the time of their observations.The average value of mag / cor ≃ 0.57 means in particular that the dipole field of TW Hya is not strong enough to spin the star down given the average mass accretion rate (Zanni & Ferreira 2013;Pantolmos et al. 2020), consistent with the rotation period of TW Hya being shorter than that of most prototypical cTTSs (e.g., CI Tau, Donati et al. 2020a).However, the dipole field is nonetheless sufficiently intense to ensure that the accretion pattern is geometrically stable, with magnetic funnels linking the inner accretion disc to the polar regions at the surface of the star (e.g., Sicilia-Aguilar et al. 2023), rather than to lower stellar latitudes through chaotic accretion tongues (Blinova et al. 2016), at least when the accretion rate is not too large.We suspect that previous reports of equator-ward accretion at the surface of TW Hya (e.g., Argiroffi et al. 2017) correspond to epochs where the accretion rate was close to its maximum.We note that no rotational modulation is detected in the 1083.3-nmHe , Pa and Br accretion / ejection proxies despite our conclusion that the accretion pattern is geometrically stable.This likely reflects the specific viewing angle of TW Hya, causing rotational modulation to be minimal and thus easily hidden behind accretion-induced intrinsic temporal variability.Besides, we report the detection of Zeeman signatures in the He , Pa and Br lines of TW Hya, suggesting that 7-20% of the line fluxes come from the hot chromospheric region at the footpoints of accretion funnels.
Last but not least, we report that RVs of TW Hya are also modulated with a period of 8.339 ± 0.008 d (or its 1-yr alias 8.147 ± 0.010 d), with a semi-amplitude of 11.1 +3.3  −2.6 m s −1 .This modulation may reflect the presence of a planet in a circular orbit around TW Hya at a distance of 0.075 ± 0.001 au, i.e., beyond both the magnetospheric and corotation radii, and that would be detected with a confidence level of 4.3 (log BF = 8.7).If the orbit of this candidate planet is coplanar with the rotation of the star, this would imply a planet mass of 0.55 +0.17 −0.13 M .An alternative option is that the RV signal we detect is caused by a non-axisymmetric density structure in the inner disc of TW Hya (possibly induced by a lower mass planet), generating a spectral contribution to some of the spectral lines (e.g., the molecular lines) and thereby inducing a small amplitude modulation in the measured RVs.We note that the period of this RV signal is slightly smaller than the photometric one that dominates the 2021 and 2023 TESS light curves; it is too early to speculate whether both are physically related (e.g., with a planet or disc structure regularly triggering enhanced accretion), or rather simply coincide by chance.Additional SPIRou observations are needed to further investigate the spectral properties of this RV signal and unambiguously diagnose its origin, before claiming the detection of a close-in massive planet orbiting TW Hya.If confirmed, this detection would demonstrate that planet formation and migration is actively going on within the protoplanetary disc of TW Hya and likely participates in generating the reported disc structures as previously suspected (e.g., Tsukagoshi et al. 2019).In particular, the candidate inner planet we report here may possibly be at the origin of the innermost gap at 1 au (Andrews et al. 2016) in the disc of TW Hya and / or contribute to the misalignment of the inner disc rings within 7 au (Debes et al. 2023).It is however unlikely to have caused the more distant multiple gaps (e.g., at about 25 and 40 au, van Boekel et al. 2017;Huang et al. 2018) in the outer disc, that may probe the presence of additional embedded massive planets within the protoplanetary disc of TW Hya.     4) get stronger and increasingly dominant as more spectra are added to the analysis, especially in the filtered RVs but also in the raw RVs (though at a lower significance level as expected).The horizontal dashed lines illustrate the transition between seasons 2020 and 2021, and between 2021 and 2022.Note the difference in color scale between both panels.

Figure 1 .
Figure 1.Longitudinal magnetic field ℓ of TW Hya (red open circles) 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 panel.The rms of the residuals is about 15 G ( 2 r = 1.13), consistent with our median error bar (of 15 G as well).The 2 r with respect to the ℓ = 0 G line is equal to 32.3.

Figure 2 .
Figure 2. Same as Fig. 1 for the 1 coefficient of the first PCA eigenvector (scaled and shifted to match ℓ ).The GPR fit and associated error bars are now shown in green.The rms of the residuals is about 9 G.

Figure 3 .
Figure 3. Magnetic parameters of TW Hya, derived by fitting our median SPIRou spectrum using the atmospheric modeling approach of Cristofari et al. (2023a), which incorporates magnetic fields as well as a MCMC process to determine optimal parameters and their error bars.We find that TW Hya hosts a small-scale magnetic field of < > = 3.60±0.04kG, with the relative area of non-magnetic regions being 0 = 7 ± 2%.

Figure 4 .
Figure 4. Observed (thick black line) and modelled (thin red line) LSD Stokes (top row) and (bottom row) profiles of TW Hya for seasons 2019 (first colum), 2020 (second) column, 2021 (third colum) and 2022 (fourth column).Observed profiles were derived by applying LSD to our SPIRou spectra, using the atomic line mask outlined in Sec. 2. Rotation cycles (counting from 0, 82, 188 and 294 for the 2019, 2020, 2021 and 2022 seasons respectively, see TableB1) are indicated to the right of all LSD profiles, while ±1 error bars are added to the left of Stokes signatures.

Figure 5 .
Figure 5. Reconstructed maps of the large-scale field of TW Hya (left, middle and right columns for the radial, azimuthal and meridional components in spherical coordinates, in G), for seasons 2019, 2020, 2021 and 2022 (top to bottom rows respectively), derived with ZDI from the Stokes and LSD profiles of Fig. 4. The maps are shown in a flattened polar projection down to latitude −10 • , 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.

Figure 6 .
Figure 6.Raw (top), filtered (middle) and residual (bottom) RVs of TW Hya (red open circles).The top plot shows the MCMC fit to the RV data, including a QP GPR modeling of the activity and the RV signature of a putative close-in planet of orbital period 8.339 ± 0.008 d (cyan), whereas the middle plot shows the planet RV signature alone once activity is filtered out.The rms of the RV residuals is 15.6 m s −1 .

Figure 7 .
Figure 7. Filtered (top plot) and residual (bottom plot) RVs of TW Hya phase-folded on the 8.34-d period.The red open circles 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. 6, the dispersion of RV residuals is 15.6 m s −1 .

Figure 8 .
Figure 8. Stacked Stokes profiles and 2D periodograms of the 1083.3-nmHe triplet (left panels) and of the 1282.16-nmPa line (right) in the stellar rest frame, for our complete data set.The color scale depicts the logarithmic power of the periodogram.Both lines are strongly variable with time, but with no obvious periodicity showing up, in particular at rot and(dashed horizontal lines), apart from peaks close to the synodic period of the Moon (29.5 d) and its 1-yr aliases (also present in the window function, see Fig.E1).

Figure 9 .
Figure 9. Weighted-average Stokes (blue) and (red) LSD profiles of the 1083.3-nmHe (left) and Pa (right) lines of TW Hya over our full data set.Zeeman signatures are clearly detected in both cases in conjunction with the emission peak.The Stokes LSD profiles are both expanded by a factor of 20× and shifted upwards (by 2 and 3 respectively) for graphics purposes.

Figure C1 .Figure D1 .
Figure C1.Stacked periodograms of the TESS light curves of TW Hya collected in March 2019 (top) and 2021 (bottom), with the TESS data binned by groups of 10 (in 2019) and 20 (in 2021) adjacent points.In these plots, each horizontal line corresponds to a color-coded periodogram, computed on an increasing number of binned points (starting from the first binned point).The color scale depicts the logarithmic power of the periodogram.At both epochs, a signal at a period of about 4 d shows up, but only dominates for a limited time before splitting itself into 2 weaker signals (one of which at the rotation period, in 2019) or vanishing (in 2021).In 2021 (and 2023, not shown), power at a period of about 9 d dominates the overall light curve, but not in 2019 where power at about 10 d quickly weakens after showing up.The dashed line depicts the rotation period.

Figure E1 .
Figure E1.Periodogram of the raw (top), filtered (middle) and residual (bottom) RVs when including a fit to the RV signal at a period of = 8.34 d in the MCMC modeling.The cyan vertical dashed lines respectively trace rot and, whereas the horizontal dashed lines indicate the 10 and 0.1% FAP levels in the periodogram of our RV data.The orange curve depicts the window function, whereas the orange vertical dashed and dotted line outline the 1-yr period, the synodic period of the Moon (at 29.5 d) and its 1-yr aliases.

Figure E2 .
Figure E2.Stacked periodograms of the raw (top) and filtered (bottom) RVs, as a function of the number of successive spectra taken into account in the Fourier analysis (starting from the first collected spectrum).The main RV signal and its main 1-yr alias (outlined with vertical dashed lines at 8.34 d and 8.15 d, see Table4) get stronger and increasingly dominant as more spectra are added to the analysis, especially in the filtered RVs but also in the raw RVs (though at a lower significance level as expected).The horizontal dashed lines illustrate the transition between seasons 2020 and 2021, and between 2021 and 2022.Note the difference in color scale between both panels.

Figure E3 .
Figure E3.Corner plot of the MCMC fit to our RV data.The yellow, red and blue regions depict the 1, 2 and 3 confidence regions, whereas the cyan dashed and dotted lines denote the optimal values and corresponding error bars.

Table 1 .
Parameters of TW Hya used in / derived from our study

Table 2 .
Results of our MCMC modeling of the ℓ curve of TW Hya 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.

Table A1 .
Abbreviations used in the paper, in alphabetical order

Table B1 .
Observing log of our SPIRou observations of TW Hya in seasons 2019, 2020, 2021 and 2022.All exposures consist of 4 sub-exposures of equal length.For each visit, we list the barycentric Julian date BJD, the UT date, the rotation cycle c and phase (computed as indicated in Sec. 2), the total observing time t exp , the peak SNR in the spectrum (in the band) per 2.3 km s −1 pixel, the noise level in the LSD Stokes profile, the estimated ℓ with error bars, the nightly averaged RVs and corresponding error bars derived by APERO and LBL, and finally the EWs of the Pa and Br emission lines with error bars (no scaling up from veiling).