A multi-kiloGauss magnetic field driving the magnetospheric accretion process in EX Lupi
Kim Pouilly Marc Audard Alexis Lavail Ágnes Kóspál
 


   


   

   Department of Astronomy, University of Geneva, Chemin Pegasi 51, CH-1290 Versoix, Switzerland

              Kim.Pouilly@unige.ch
        Konkoly Observatory, HUN-REN Research Centre for Astronomy and Earth Sciences, MTA Centre of Excellence, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary
        Institute of Physics and Astronomy, ELTE Eötvös Loránd University, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary
        Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, IRAP/UMR 5277, 14 avenue Edouard Belin, 31400 Toulouse, France
             

   



 
  
  
   EX Lupi is the prototype of EX Lup-type stars, meaning classical T Tauri stars (cTTSs) showing luminosity bursts and outbursts of 1 to 5 magnitudes lasting for a few months to a few years. These events are ascribed to an episodic accretion that can occur repeatedly but whose physical mechanism is still debated.
  
   In this work, we aim to investigate the magnetically-driven accretion of EX Lup in quiescence, including for the first time a study of the small and large-scale magnetic field. This allows us to provide a complete characterisation of the magnetospheric accretion process of the system.
  
   We use spectropolarimetric times series acquired in 2016 and 2019 with the Echelle SpectroPolarimetric Device for the Observation of Stars and in 2019 with the SpectroPolarimètre InfraRouge at the Canada-France-Hawaii telescope, during a quiescence phase of EX Lup. We were thus able to perform a variability analysis of the radial velocity, the emission lines and surface averaged longitudinal magnetic field along different epochs and wavelength domains. We also provide a small-scale magnetic field analysis using Zeeman intensification of photospheric lines and large-scale magnetic topology reconstruction using Zeeman-Doppler Imaging.
  
   Our study reveals a typical magnetospheric accretion ongoing on EX Lup, with a main accretion funnel flow connecting the inner disc to the star in a stable fashion and producing an accretion shock on the stellar surface close to the pole of the magnetic dipole component. We also measure one of the strongest fields ever observed on cTTSs. Such a strong field indicates that the disc is truncated by the magnetic field close but beyond the corotation radius, where the angular velocity of the disc equals the angular velocity of the star. Such a configuration is suitable for a magnetically-induced disc instability yielding episodic accretion onto the star.
  
   

   

   A multi-kiloGauss magnetic field driving the magnetospheric accretion process in EX Lupi
    K. Pouilly1
          M. Audard1
          Á. Kóspál23
          A. Lavail4
          
    Received 16 July 2024; accepted 04 September 2024
==========================================================================================






§ INTRODUCTION

   
    EX Lup-type objects (EXors) are classical T Tauri stars (cTTSs), meaning low-mass pre-main sequence stars surrounded by an accretion disc, that show burst and outburst events ascribed to episodic accretion. 
    During these phases, they can increase their optical luminosity from 1 to 5 magnitudes, lasting typically for a few months to a few years, and can occur repeatedly <cit.>. 
    These events are therefore more moderate, both in duration and luminosity increase, than the FU Orionis-type stars (FUors).
    
    While the magnetospheric accretion of cTTSs, where the strong stellar magnetic field truncates the disc, forcing the accreted material to follow the magnetic field lines <cit.>, seems to be ongoing on EXors as well, the origin of this episodic accretion is still highly debated. 
    The different hypotheses can be gathered in three groups <cit.>: (i) the magnetospheric accretion itself, being inherently episodic when a strong magnetic field truncates the disc close to the corotation radius <cit.>. (ii) The disc, showing viscous-thermal <cit.>, or gravitational and magneto-rotational <cit.> instabilities, or accretion clumps in a gravitationally unstable environment <cit.>. (iii) The presence of a companion, perturbing the accretion trough tidal effect <cit.>,  or thermal instabilities <cit.>.

    This work aims to investigate the magnetospheric accretion process of the prototypical EXors, EX Lup, using high-resolution spectropolarimetric time series, as it was done for cTTSs not displaying episodic accretion <cit.>.
    This object is a young M0.5-type star <cit.>, known to have both moderate and short-timescale variability and rare extreme episodic outbursts.
    It is located at 154.7±0.4 pc <cit.> and has a rotation period of 7.417 days, from its radial velocity modulation that was first ascribed to a low-mass companion <cit.>, before being ascribed to stellar activity <cit.>.
    The system has a moderate-to-low inclination of its rotation axis (between 20^∘ and 45^∘) according to modelling of the spectral energy distribution <cit.> and emission lines analysis <cit.>, and a projected rotational velocity vsin i=4.4±2.0  <cit.>.
    
    This object is one of the most studied EXors using spectroscopy, both in quiescence <cit.> and in outburst <cit.>, but the present study is the first one including spectropolarimetry, giving access to information on the magnetic field together with the accretion diagnostics.
    Its spectrum contains the typical accretion-related emission lines observed in cTTSs, in addition to numerous neutral metallic emission lines superimposed to photospheric absorption in quiescence, and overwhelming any absorption feature in outburst <cit.>.
    The present dataset was acquired during quiescence, meaning that we will characterise the "stable" accretion of the system, even if <cit.> have shown that the accretion pattern seems stable in quiescence and outburst, only the amount of accreted material is affected.
    This accretion pattern is consistent with the typical cTTS magnetospheric accretion, through accretion funnel flows connecting the disc to the stellar surface, except that "clumps" of material are also accreted through these funnel flows, detected thanks to the day-to-day variation of the emission lines' broad component (BC) <cit.>.

    This article is organised as follows: we describe the observations in Sect. <ref>, the analysis and results are presented in Sect. <ref> and discussed in Sect. <ref>.
    We conclude this work in Sect. <ref>.
    




§ OBSERVATIONS



    The spectropolarimetric time series used in this work were acquired at the Canada-France-Hawaii Telescope at two different epochs (2016 and 2019, proposals 16AF03 and 19AF50, respectively).
    The second data set is composed of two subsets using two different instruments: the Echelle SpectroPolarimetric Device for the Observation of Stars <cit.> in the optical and the SpectroPolarimètre InfraRouge <cit.> in the near-infrared, both used in polarimetric mode, while the first one only used ESPaDOnS.
    This means that each observation is composed of four sub-exposures taken in different polarimeter configurations, which are then combined to obtain the intensity (Stokes I), the circularly polarised (Stokes V), and the null polarisation spectra.
    A complete journal of observation is provided in Table <ref>.

    

 §.§ ESPaDOnS

    The ESPaDOnS observations, which cover the 370 to 1050 nm wavelength range and reach a resolving power of 68 000, consist of 11 nights between 2016 June 09 and 2016 June 24, with an approximately nightly cadence, and 6 nights between 2019 May 31 and 2019 June 12, the 5 latter respecting a 1-day sampling.
    The signal-to-noise ratio (S/N) of the 2016 (2019) observations ranges between 69 and 142 (111 and 140) for the Stokes I at 731 nm. 
    Each observation was reduced using the  package <cit.>.

    

 §.§ SPIRou

    The SPIRou observations are covering the 960 to 2350 nm wavelength range with R∼75 000.
    They were acquired during 8 consecutive nights, between 2019 June 14 and 2019 June 21, and the S/N of the unpolarised spectra in the H-band range between 107 and 175.
    The observations were reduced using the  pipeline <cit.>.

    
    

    
    




§ RESULTS


    
    In this section, we present the results obtained from the analysis of the observations described in Sect. <ref>. They consist of the analysis of the radial velocity, emission lines, and stellar magnetic field.
    
    

 §.§ Radial Velocity

    
        To determine the radial velocity of EX Lup, we cross-correlated each spectrum with a synthetic spectrum computed using the  code <cit.>, with  atmospheric models <cit.> and  <cit.> line lists adapted to EX Lup stellar parameters for ESPaDOnS and SPIRou wavelengths.
        We computed the cross-correlation function (CCF) over 27 (14) wavelength windows of about 10 nm ranging from 441 to 890 nm (1150 to 2290 nm) for ESPaDOnS (SPIRou) observations.
        Then we performed a sigma-clipping across all the CCFs for each observation and used the mean and standard deviation of the remaining values as measurement of the radial velocity and its uncertainty.
        The results are plotted in Fig. <ref>.
        A quick sinusoidal fit of the values for each data set allowed us roughly measure the periodicity and mean value of the radial velocity and yielded P=7.55±0.24 d and ⟨⟩=-0.48±0.07  (ESPaDOnS 2016), P=7.63±0.18 d and ⟨⟩=-0.59±0.14  (ESPaDOnS 2019), and P=8.36±0.97 d and ⟨⟩=-0.58±0.14  (SPIRou).
        These measures are consistent within the uncertainties with previous results obtained by <cit.>, =7.417±0.001 d and ⟨⟩=-0.52±0.07 , we will thus adopt the latter values for our analysis.
        Finally, we folded all the measurements in phase using =7.417 d and an arbitrary T0, and we fitted this curve with a sinus to estimate the T0 required to set ϕ = 0.5 at the mean velocity between the maximum and minimum of the modulation <cit.>.
        The resulting T0 is HJD 2 457 544.40981.
        We will thus use the following ephemeris for the rest of this work:
        
    HJD (d) = 2 457 544.40981 + 7.417 E,

        where E is the rotation cycle.
        All our radial velocity measurements are in phase (Fig. <ref>), but the 2019 measurements show a larger amplitude of its modulation.
        This indicates an evolving feature modulating the radial velocity but located at the same longitude in 2016 and 2019.
        The values are summarised in Table <ref>.

        

        

       
    

 §.§ Emission line variability

    
        In this section, we present the analysis of emission lines tracing either the accretion funnel flow, here the Balmer lines <cit.>, or the accretion shock such as the CaII infrared triplet (IRT), the HeI D3 <cit.>, or the HeI at 1083 nm lines.
        For each line, we analysed the profile variability, their periodicity (except for the 2019 ESPaDOnS data set which does not cover a sufficient time span), and the correlations of these variabilities.
        Most of the analyses of this section were performed using [<https://github.com/pouillyk/PySTELLA>], a Python tool for SpecTral Emission Lines (variabiLity) Analysis.
        
        

  §.§.§ Balmer lines

            The Balmer lines are partly formed in the accretion funnel flow and are thus tracing the magnetospheric accretion process.
            In this work, we focussed on Hα, Hβ, and Hγ, which we corrected from the photospheric contribution using the moderately active M-dwarf HD 42581 as template <cit.>, broadened to the vsin i of EX Lup.
            The 2016 and 2019 profiles are shown in Fig. <ref> and Fig. <ref>, respectively.
            On both data sets, the profiles are composed of a broad and a narrow component, both highly variable. 
            Furthermore we can notice a flux depletion, going below the continuum, around +200 , and extending up to +300 .
            This behaviour is characteristic of the so-called inverse P Cygni (IPC) profiles, the red-shifted absorption forming due to infalling material.
            
            The 2D-periodograms, consisting of a Lomb-Scargle periodogram computed in each velocity channel, of the 2016 lines are presented in Fig. <ref>, and show a periodic signal along the whole Hβ and Hγ lines consistent with the rotation period of the star, with a false alarm probability <cit.> of 0.04 and 0.01, respectively. 
            The Hα line is showing this signal in a less continuous way with a much higher FAP (0.25).
            The symmetric signal around 0.9 d^ -1 is the 1-day alias, a spectral leakage of the Fourrier transform reconstructing (with the real period) the observation sampling, meaning approximately one observation per day.
            
            

            To separate the various parts of the line profile undergoing different variability patterns, we computed the auto-correlation matrices of each line.
            This tool consists of computing a linear correlation coefficient (here a Pearson coefficient) between the velocity channels of the line.
            A correlated region (close to 1) indicates a variability dominated by a given physical process. 
            An anti-correlated region (close to -1) indicates a variability dominated by a given physical process or, at least, linked physical processes.
            The auto-correlation matrices of 2019 Hα, Hβ, and Hγ lines are shown in Fig <ref>.
            Here again, Hα is showing a different behaviour than Hβ and Hγ.
            
            
            
            Hα shows three main correlated regions between -100  and +100 , corresponding to the core of the line, around +150 , corresponding to a slight emission excess in the IPC profile (occurring around HJD 2 458 8634.95 and 2 458 641.95), and around +250 , corresponding to the IPC profile itself. 
            The two latter regions are anti-correlated between them, but the ∼+150  region is also slightly anti-correlated with the line center.
            This means that this region might be a broadening of the BC invoked by <cit.>, causing a global decrease of the line (and so a anti-correlation with the whole line profile).
            
            
            Hβ and Hγ are showing the same correlation around the line centre and +350 , without significant anti-correlation around +150 , and a moderate anti-correlation of the profile with a region around +350  which seems to be an artefact as it is located in the continuum.
            
            

        

  §.§.§ HeI D3 587.6 nm

        
            The HeI D3 lines of 2016 and 2019 data sets are presented in Fig. <ref> and are composed of a narrow component (NC) only, extending from approximately -35 to  +50 , without significant BC.
            The NC is formed in the post-shock region of the accretion spot <cit.>, and can thus trace the accretion close to the stellar surface. 
            One can note the strong variability of this NC with a maximum reached around ϕ=0.6 [HJD 2 457 548.91, 2 458 645.88, and 2 458 646.87].
            
            The auto-correlation matrices shown in Fig. <ref> confirm that this region is formed by one physical process, and the periodicity of the 2016 NC's variation, consistent with the stellar rotation period (see Fig. <ref>, FAP=0.06), indicates that this component is tracing an accretion shock at the stellar surface.

            

            

            We thus performed a fit of the HeI D3 NC's radial velocity following the method described in <cit.> to recover the emitting region's location.
            The results are shown in Fig. <ref> and summarised below:
            
                
  * V_ flow = 7.091_-1.03^+2.00 ,
                
  * V_ rot = 3.85 ± 5.0 ,
                
  * dϕ = 0.10_-0.11^+0.15,
                
  * θ = 12.05_-6.7^+48.00 ^∘,
                
  * α = 55.40_-16.36^+20.00 ^∘,
            
            where V_ flow is the velocity of the material in the post-shock region, V_ rot the equatorial velocity, 0.5+dϕ the phase where the emitting region is facing the observer, θ the colatitude of the spot, and α = 90^∘-i, i being the inclination of the rotation axis.
            This means that the emitting region is located at ϕ = 0.6,  70^∘ latitude.
            
            These results are consistent within uncertainties with the measurements by <cit.>, who give for the HeI emitting region[Values estimated from their Fig. 7.] a latitude of 60±25^∘, longitude 40±5^∘, meaning  ϕ = 0.1±0.01.
            However, the authors are using the first date of observation as T0, translated to our ephemeris (Eq. <ref>), this yields ϕ = 0.7.
            
            
        

  §.§.§ CaII infrared triplet

            As the HeI D3 NC, the CaII IRT NC is formed in the post-shock region, we thus studied these lines as well. 
            The three components of this triplet show identical shape and variability, we thus focussed on one of them located at 854.209 nm. 
            The 2016 and 2019 residual line profiles are shown in Fig. <ref> and Fig. <ref>, respectively.
            The two sets of lines show an IPC profile around 200 , which is, for the 2016 line, periodic with the stellar rotation period (see the 2D-periodogram in Fig. <ref>, FAP=0.05).

            The 2016's auto-correlation matrix (Fig. <ref>) exhibits several correlated regions: from -130 to -60, -50 to -10, +50 to +100, +110 to +160, and +170 to +200 .
            However, these regions are less numerous on the 2019 matrix (Fig. <ref>), with only three correlated regions from -40 to +40, +80 to +130, and +150 to +200 , but we retrieve the correlated region around the IPC profile in both matrices, which is anti-correlated with the NC in 2019.
            
            This goes with the smaller (larger) variability of the NC (BC) observed in 2016 compared to 2019, showing the different origins of the two components and a small change in the accretion pattern between the two epochs.
            
            
            
            
        

  §.§.§ HeI 1083 nm

            The only accretion tracer in emission in EX Lup's SPIRou observation is the HeI line at 1083 nm. 
            The profiles, the 2D-periodogram and the auto-correlation matrix are shown in Fig. <ref>.
            The profiles seem composed of two peaks blue- and red-shifted around -50 and +100, and two absorptions, blue-and red-shifted at higher velocity (-150 and +200).
            
            The red-shifted peak and the two absorptions display significant variability and seem modulated on the stellar rotation period with FAPs reaching 0.001, 0.01, and 0.02 for the red-shifted peak, the blue- and the red-shifted absorption, respectively.
            
            The auto-correlation matrix revealed a more complex decomposition.
            Indeed, if the four substructures seen in the profiles are represented, it seems that the two absorptions can both be separated in two regions, from -250 to -160 and -160 to -110  for the blue-shifted absorption, and from +130 to +170  and +210 to +270  for the red-shifted absorption. 
            Furthermore, the main peak at ∼+100  is anti-correlated with the most blue- and red-shifted regions only.
            
            This can be interpreted as follows: the blueshifted absorption, probably a P-Cygni profile traditionally ascribed to a wind, is also associated with a redshifted emission excess, producing the two substructures seen in the redshifted absorption. 
            The IPC profile, as the opposite physical phenomenon, is also associated with a blueshifted emission excess, responsible for the two substructures in the blueshifted absorption.

        

  §.§.§ Correlation matrices ESPaDOnS

        
            As the several lines studied are tracing different regions of the accretion, we can compute correlation matrices between two different lines to analyse the link between the different regions identified from the auto-correlation matrices.
            The correlation matrices of the 2016 and 2019 lines are presented in
            
            Appendix <ref>.

            In 2016, the Hα line centre is correlated with the HeI D3 and the CaII IRT NCs, and slightly anti-correlated with the region of CaII IRT IPC profile.
            The latter is also strongly anti-correlated with the HeI D3 NC.
            The region of the Hα's IPC profile is also anti-correlated with the HeI D3 NC, and correlated with the region of the CaII IRT IPC profile.
            In 2019 matrices, we observe the same behaviour between the NCs and IPC regions of the various lines, but both the correlation and anti-correlation coefficients are stronger.
            
            
            
            
            
            
            
            
            
            
            
            
            
        
    

 §.§ Magnetic field

        In this section, we present the magnetic analysis of EX Lup.
        This was done at two scales: the large scale using the Zeeman-Doppler Imaging technique <cit.>, and the small scale using the Zeeman intensification of photospheric lines.
        
        

  §.§.§ Large-scale

        
            In this section, we used the Least-Squares Deconvolution method <cit.> to study the large-scale magnetic field.
            This method allows us to increase the S/N of the Stokes I (unpolarised) and Stokes V (circularly polarised) profiles, by using as many photospheric lines as possible.
            To compute the LSD profiles we used the [<https://github.com/folsomcp/LSDpy>] Python implementation.
            We normalised our LSD weights using an intrinsic line depth, a mean Landé factor and a central wavelength of 0.2, 1.2, and 500 nm (respectively) for ESPaDOns, and 0.1, 1.2, and 1700 nm for SPIRou observations. 
            The photospheric lines were selected from a mask produced using the same  line list and  atmospheric models as in Sect. <ref>, and removing the emission lines, the telluric and the heavily blended lines regions using the [<https://github.com/folsomcp/specpolFlow>] Python package. 
            About 12 000 lines were used for ESPaDOnS, and 1600 for SPIRou observations.
            The LSD profiles of ESPaDOnS 2016, 2019, and SPIRou observations are presented in Figs. <ref>, <ref>, and <ref>, respectively, and the S/N of the profiles are available in Table <ref>. 
            Please note that the observation at HJD 2 457 551.87135 was remove from this analysis because of its low S/N.
            
            A clear Stokes V signature is detected for all ESPaDOnS observations and most of the SPIRou observations.
            Furthermore, one can note that the SPIRou Stokes V signatures are much weaker than ESPaDOnS and the signal at ϕ∼0.6–0.7 almost vanishes on SPIRou when maximal on ESPaDOnS.
            
            
            

            
            
            

            Directly from the LSD profiles, one can estimate the surface averaged longitudinal magnetic field <cit.>:
            
    B_ℓ = -2.14 × 10^11×∫ vV(v)dv/λ g c ∫ (1-I(v))dv,

            where B_ℓ is in Gauss, v the velocity relative to the line centre, and λ and g  the central wavelength and the mean Landé factor used for the LSD computation.
            The integration was performed on a ±25  (±35 ) velocity range around the stellar rest frame to minimise the uncertainties without losing any magnetic information on ESPaDOnS (SPIRou) observations.
            The B_ℓ curves are shown in Fig. <ref>.
            The three curves are modulated with the stellar rotation period and show a maximum around ϕ=0.6.
            In the optical frame, the modulation's amplitude is slightly larger in 2019 than in 2016, which is reminiscent of the radial velocity behaviour (see Sect. <ref>).
            Finally, as expected from the Stokes V signatures, the SPIRou measurements are far weakest, and mostly negatives.

            The same analysis can be performed on the NC of the HeI D3 line (on a ±50  velocity range), formed close to the accretion shock and thus gives access to the magnetic field strength at the foot of the accretion funnel flow.
            The B_ℓ obtained are shown in Fig. <ref> and range between -1.4 and -3.7 kG in 2016 and between -0.5 and -4.5 kG in 2019.
            As for the radial velocity (see Sect. <ref>), the amplitude of the modulation is larger in 2019, but both curves are in phase, with a minimum reached at ϕ=0.6, which is also consistent with the emitting region's position obtained from the radial velocity modulation of the HeI D3 NC (see Sect. <ref>).
            These measurements are also opposed in phase and sign with the LSD B_ℓ. This is a common behaviour on cTTSs <cit.> reflecting that LSD and HeI D3 NC are two diagnostics probing different regions of different polarities.
            
            
                        
            
            
            Finally, we performed a complete ZDI analysis of the three data sets using [<https://github.com/folsomcp/ZDIpy>]<cit.> with a recent implementation of Unno-Rachkovsky’s solutions to polarised radiative transfer equations in a Milne-Eddington atmosphere <cit.> presented in <cit.> to used a more general description than the weak-field approximation used initially in . 
            The magnetic topology is reconstructed in two steps: (i) the building of a Doppler image (DI), starting from a uniformly bright stellar disk and iteratively adding dark and bright features to fit the observed LSD Stokes I profiles, and (ii) fitting the LSD Stokes V profiles to derive the magnetic topology by adjusting its spherical harmonic components <cit.>, here with a maximum degree of harmonic expansion ℓ_ max=15.
            We used as input parameter P_ rot=7.417 d, vsin i=4.4 .
            Then we ran a grid of ZDI Stokes V reconstruction over the inclination range derived in the literature (20-45^∘) and used the minimal χ^2 obtained to set our input inclination value (i=30^∘).
            The resulting maps are presented in Fig. <ref> and <ref> and the corresponding LSD profile fits are shown in Appendix <ref>.

            
            
            
            
            The ESPaDOnS 2016's DI map shows a main dark spot around ϕ=0.55 and extending between 70 and 20^∘ latitude, which is consistent with HeI D3 emitting region, and a bright structure from ϕ=0.95 to 0.20 around the equator which is certainly a plage as even the accretion shock is dark at the photospheric level.
            The magnetic topology, mostly toroidal (61%) is dominated by the dipolar (60%) and the quadrupolar (17%) components, with a mean magnetic field strength of 1.00 kG (B_ max = 3.12 kG), and a magnetic dipolar positive pole of 0.728 kG located at about 30^∘ latitude and 166^∘ longitude (ϕ=0.46).

            
            
            
            
            
            
            
            
            
            The brightness map obtained from SPIRou revealed a main dark feature extended from ϕ = 0.0 to 0.25, and a bright plage around ϕ = 0.6, which is perfectly consistent with the 2016 map.
            However, the magnetic topology seems less complex than in 2016, almost fully poloidal (99%) and more dominated by the dipole component (77%).
            As expected, the recovered field strength is much smaller (⟨B⟩ = 0.131 kG, B_ max = 0.389 kG).
            The dipole negative pole is located at about 23^∘ latitude, 329^∘ longitude (ϕ=0.91), with a -0.231 kG-strength.
            
            
            

            
            
            
            
            
            
            
            
            
            
            
            Given the poor rotational phase coverage of the ESPaDOnS 2019's data set (ϕ∈[0.03,0.64]), we needed to guide the reconstruction instead of starting from a uniform map to avoid a too strong extrapolation on the missing phases.
            We do not expect a similar brightness contrast between ESPaDOnS and SPIRou, but given the very similar spectroscopic behaviour between the 2016 and 2019 ESPaDOnS data sets (see Sect. <ref> and <ref>), we could expect similar features, we thus used the ESPaDOnS 2016 brightness map as input to guide its reconstruction.
            To check this assumption, we reproduce the Stokes I profiles resulting from the 2016 reconstruction on the 2019 phases, yielding a consistent behaviour (reduced χ^2=1.1).
            The obtained final brightness reconstruction shows a main dark feature, slightly shifted in phase compared to 2016 (ϕ≈0.5) and less extended in latitude.
            A polar bright feature is also located around ϕ=1, which is reminiscent of the SPIRou's maps.
            Concerning the magnetic reconstruction, we expect a similar topology of the magnetic field between SPIRou and ESPaDOnS 2019, with different magnetic strengths as pointed out by the B_ℓ analysis. 
            We thus used the SPIRou magnetic maps as input to guide the reconstruction.
            The resulting topology is less poloidal-dominated than SPIRou (91%), and the dipolar component occupies a smaller fraction of this poloidal field (65%) with a similar contribution of the quadrupolar and octupolar components (about 15%).
            Surprisingly the mean magnetic field strength is similar to 2016 (1.08 kG) but the maximum field strength is much higher (4.79 kG), as well as the strength of the dipolar pole B_ dip=1.87 kG which is located at the same position (ϕ =0.43 and 30^∘ latitude).
                
            
            
            
            
            
            
            
            
            

            
              
            
            
            
        

  §.§.§ Small-scale

            Even if the ZDI analysis gives access to the magnetic topology, it neglects the small-scale magnetic field, which contains a major part of cool stars’ magnetic energy.
            This is analysed by studying the change it induces in the shape of magnetically sensitive lines <cit.>, a technique called Zeeman intensification.
            To perform this analysis, we used the algorithm of <cit.>, performing a Markov Chain Monte Carlo (MCMC) sampling, using the  library <cit.> on a grid of synthetic spectra produced by the  code, a polarised radiative transfer code described by <cit.>.
            This grid is computed from the  line lists used in Sect. <ref>, and  atmospheric models <cit.>.
            We parametrized a uniform radial magnetic field as a sum of magnetic field strength ranging from 0 to 6 kG, with a 2 kG step, weighted by filling factors representing the amount of stellar surface covered by this magnetic field.
            For ESPaDOnS observations, as in previous studies <cit.>, we used the 963.5–981.2 nm region.

            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
            
                 
            
            
            
            
            

            
            
            This region contains a group of TiI lines with different magnetic sensitivity (g_ eff summarised in Table <ref>), allowing us to disentangle the effect of the magnetic field on the equivalent widths from the effect of any other parameters, such as the TiI abundance. 
            However, this region also contains many telluric lines which are superimposed on the TiI lines from the stellar spectra and which need to be removed from the observed spectrum to perform the magnetic analysis. 
            To do so, we used the molecfit package <cit.>, developed to model and remove telluric lines from spectra obtained with instruments at the European Southern Observatory, and which can be used on spectra from any instrument.

            Finally, as EX Lup is an accreting star, showing a signature of an accretion shock, the veiling might disturb the inference results, in particular the abundance, the vsin i and the radial tangential macroturbulent velocity, v_ mac.
            We thus estimated the veiling using the magnetic null line of the TiI multiplet at 974.36 nm by performing a χ^2 minimisation using  synthetic spectra letting the abundance, the vsin i, and the v_ mac as free parameters and adding a fractional veiling defined as:
            
    I_veil = I+r/1+r,

            where I_ veil is the veiled spectrum, I the spectrum without veiling, and r the fractional veiling.
            We stress to the reader that deriving a precise value of the veiling is outside the scope of the present study
            
            , our aim was only to get an estimation of the EX Lup mean spectrum at a wavelength close to the TiI lines in order to minimise its effect on the other inferred parameters.
            For the 2016's ESPaDOnS mean spectrum, the minimal χ^2 is reached for r=0.5, vsin i=4.9 , v_ mac=1.9 , and a TiI abundance of -7.2.
            For the mean spectrum of the 2019's ESPaDOnS observations, we obtained r=0.48, vsin i=5.0 , v_ mac=2.5, and a TiI abundance of -7.2.
            We will thus assume the values obtained for r and v_ mac, and use the others as initial guesses for the MCMC sampling.

            We assumed a multicomponent model given by:
            
    S=∑f_i S_i,

            where f_i are the filling factors, meaning the fraction of the stellar surface covered by a field strength B_i, and S_i the synthetic spectra of the corresponding magnetic field strength.
            The averaged magnetic field is thus given by:
            
    ⟨ B ⟩ = ∑f_i B_i.

            To set the number of filling factor to use, we iteratively added filling factors, with a 2-kG step in the corresponding magnetic field strength, and used the Bayesian information criterion <cit.> to only include the filling factors that significantly improve the fit.
            The suitable solutions are components of 0, 2, 4 kG for ESPaDOnS 2016 and SPIRou observations, and 0, 2, 4, 6, 8 kG for ESPaDOnS 2019.
            The free parameters of the analysis are thus the following: f_i, vsin i, v_r, and the TiI abundance, for which uniform prior were adopted.
            Finally, we used an effective sample size of 1000 <cit.>.
            
            The resulting line fit and magnetic field strength posterior distributions are presented in Fig. <ref>.
            The 2016's magnetic field strength (3.08 ± 0.04 kG) is consistent with 2019 within uncertainties (3.16 ± 0.05 kG). 
            
            The inferences of all parameters are summarised in Table <ref>.
            
    

            


            For SPIRou observations, we used here again a set of TiI lines around 2200 nm (see Table <ref>). 
            Unfortunately, this wavelength region does not contain any magnetically null line deep enough at our S/N to perform the veiling study we have done on ESPaDOnS observations.
            Only a visual inspection of the line at 974.4 nm is possible, indicating that at this wavelength the parameters obtained for ESPaDOnS are consistent with the SPIRou observations.
            We thus used the veiling values obtained for the 2019 ESPaDOnS data set, fixed the vsin i at the literature value and let the inference compensate for the eventual error with the non-magnetic parameters.
            The inferred parameters are summarised in Table <ref>, and the line fit and inferred magnetic field strength are shown in Fig. <ref>. 
            
            The magnetic field strength recovered (2.00±0.03 kG) is significantly lower than the values found for ESPaDOnS observations.
            As highlighted by <cit.>, the small-scale field might be overestimated when using optical wavelength. 
            In our case, a second explanation might come from the high v_ mac obtained, probably needed to compensate for an underestimated veiling, which lowers the effect of the magnetic field.
            
            
            
            
            
            



§ DISCUSSION


 
    This work aimed at characterising the accretion process of the prototypical EXor, EX Lup, in quiescence, together with its magnetic field at small and large scales.
    This study confirms that the typical magnetospheric accretion process of cTTSs is ongoing on this system, which seems stable between the two epochs studied (2016 and 2019), with a main accretion funnel flow connecting the disc to the stellar surface.
    This produces the IPC profile observed in the H lines studied and their modulation with the stellar rotation period.
    The accretion shock at the stellar surface produces the emission of the NC of the HeI D3 and CaII IRT lines, which are modulated with the stellar period and anti-correlated with the IPC profiles as expected in the magnetospheric accretion scheme.
    This is consistent with the maximum IPC profile occurring at the same phase as the HeI D3 emitting region (ϕ≈ 0.6), indicating a funnel flow aligned with the accretion shock and inherently meaning that the truncation radius is located at the stellar corotation radius.
    As expected, this phase is also associated with an extremum of the B_ℓ, for both LSD and HeI D3 measurements, and with the ESPaDOnS ZDI brightness and magnetic topology reconstructions,  showing the connection between accretion and magnetic field.
    To investigate the truncation radius r_ mag, we used the expression given by <cit.>:
    
    r_ mag/R_⋆  = 2 m_ s^2/7 B_⋆^4/7Ṁ_ acc^-2/7 M_⋆^-1/7 R_⋆^5/7,

    where the Mach number m_s≈1, B_⋆ is the equatorial magnetic field strength <cit.> in units of 140 G,  is the mass accretion rate in units of 10^-8 , M_⋆ the stellar mass in units of 0.8 M_⊙, and R_⋆ the stellar radius in units of 2 R_⊙.
    As  in quiescence, we used the value before the 2022 outburst given by 
    
    <cit.>, 1.8×10^-9 .
    
    
    The stellar mass and radius are given by <cit.> (0.6 M_ and 1.6 R_), and magnetic obliquity is obtained from ZDI analysis (30^∘ for ESPaDOnS, 23^∘ for SPIRou, see Sect. <ref>).
    The dipole strengths needed to obtain r_ mag = r_ corot = 8.5±0.5 R_⋆ are thus 
    
    2.10±0.46 kG and 1.98±0.43 kG for ESPaDOnS and SPIRou, respectively.
    
    Only an estimate of the truncation radius can be obtained because we do not have a precise value of the dipole strength.
    Indeed, the B_ℓ measurements in the HeI D3 line and Zeeman intensification values also contain the higher-order components of the magnetic field, and the ZDI spherical harmonic decomposition from LSD might include flux cancellation lowering the large-scale strength obtained.
    Even if the ZDI results on the 2019 datasets indicate a topology largely dominated by the dipole component, using these values should yield an overestimation of r_ mag.
    The results using our magnetic field measurements are summarised in Table <ref>.
    
    The values obtained from the B_ℓ and the small-scale field are not consistent with the corotation radius, as expected, except for the SPIRou measurement which can be explained by the smaller recovered field in the infrared domain. 
    However, the optical ZDI values yield truncation radii consistent with the corotation radius, due to the dipole-dominant magnetic topology of the system allowing us to recover a good estimate of the dipole pole strength using ZDI.
    Here again, the SPIRou value yields a lower truncation radius due to the lower magnetic field strength obtained.

    The difference in the results of the various magnetic analyses between the optical and infrared frames in 2019 needs to be discussed. 
    From the B_ℓ computed from the LSD profiles, the maximum values obtained are 911±26 and 64±20 G, for ESPaDOnS and SPIRou observations, respectively. 
    These maxima both occurred around ϕ=0.6, where the Stokes V signature is maximal with ESPaDOnS but almost vanishes with SPIRou, it is thus not surprising to have a large discrepancy here.
    Concerning the minimum B_ℓ values, they reach -155±11 and -178±13 G, for ESPaDOnS and SPIRou observations, respectively, both around ϕ=0.1 and are thus consistent.
    This behaviour is also visible on the ZDI reconstruction, where the strong positive radial magnetic field region at phase 0.6 on ESPaDOnS map completely vanishes on SPIRou map.
    The two qualitative explanations we can provide are the following: (i) this maximum value is located in an optically dark region of the photosphere, which is consistent with simultaneous optical photometry (Kóspál et al., in prep.), and thus less contrasted in the SPIRou domain.
    The magnetic contribution of this region to the Stokes V signal might be thus lowered in the infrared domain if the smaller-scale negative field is obscured in the optical domain.
    And (ii) the two wavelength domains allowing to trace different heights in the photosphere, that such a difference might be an indication of a vertical structure of the magnetic field.

    
    
    

    Concerning the two ZDI reconstructions obtained from the two ESPaDOnS epochs, they seem to point a similar brightness and magnetic radial maps even if the overall topology has evolved from a toroidal- to a poloidal-dominated state, as seen on previous objects <cit.>. 
    The strong positive radial field spot, associated with the dark spot of the stellar brightness and with the dipole pole, is also consistent with the modulation of HeI D3 NC radial velocity and B_ℓ.
    However, the latter is pointing to an accretion shock associated with a region with a strong negative field.
    Given the large amplitude of the HeI D3 NC radial velocity variation compared to EX Lup vsin i, the emitting region has probably a very small extent, meaning that its magnetic field information is lost at the photospheric level, explaining this disparity <cit.>.  

    Despite this stable pattern between 2016 and 2019, we observed some disparities in some parameters.
    Even if the radial velocity and the B_ℓ (from LSD and HeI D3) modulations are well in phase, their amplitudes are slightly larger in 2019.
    If the stronger magnetic field obtained in the optical frame, at small- and large-scale, explains the larger amplitude of the B_ℓ modulation, the accompanying effect on the radial velocity points to a modulation by the hotspot.
    This is not a surprising behaviour, but the consistency between the ESPaDOnS and SPIRou radial velocity measurements is.
    Indeed, the stellar activity effect on the photospheric lines, producing the apparent radial velocity modulation, is a wavelength-dependent phenomenon, an amplitude consistency between the radial velocity measurements in the optical and infrared frames is thus not expected.
    
    To investigate in our data sets the relation between the stellar activity and the radial velocity modulation, we compared the latter to an activity indicator, the bisector inverse slope <cit.>.
    The BIS was computed from the LSD profiles presented in Fig. <ref>, <ref>, and <ref>, and is defined as the difference between the mean velocity of the bisector at the top and at the bottom of the line.
    On ESPaDOnS profiles, the first 15%, which contains the continuum and the wings, were ignored, as well as the last 15%, where the noise or an activity signature splitting the profile into two parts can affect the computation.
    The top and bottom regions used to compute the BIS are the 25% top and bottom parts of the remaining profile.
    For SPIRou profiles, the same conditions were used except that we had to ignore the first 25% of the profiles as it was affected by stronger wings.
    In the case where the radial velocity modulation is only induced by the stellar activity, the line deformation, indicated by the BIS, is completely responsible for this modulation.
    This means that a strong linear correlation should appear between the BIS and the radial velocity with, in a perfect situation, a -1-slope and a BIS = 0  at the mean velocity.
    
    The BIS versus radial velocity plots are presented in Appendix <ref>.
    For each data set, the Pearson correlation coefficient indicates a strong anti-correlation (2016: r=-0.81, p-value=0.005; ESPaDOnS 2019: r=-0.95, p-value=0.004; SPIRou: r=-0.72, p-value=0.05).
    However, the optical measurements show much shallower slopes (-0.51±0.13 and -0.43±0.07 in 2016 and 2019, respectively), and an intersect far from the mean velocity expected (-0.34±0.10 and -0.08±0.06). 
    Only the infrared measurements are consistent (slope:-0.68±0.27, intersect: -1.37±0.25), but only thanks to the larger uncertainties, probably due to the lower correlation.
    Furthermore, the much lower BIS obtained (exclusively negative) in the infrared compared to the optical frame, indicates chromatic effects that are not observed in the radial velocity modulation.
    Stellar activity is thus probably dominating the radial velocity modulation, but another effect, such as Doppler shift induced by a companion, can not be completely excluded and needs further investigations that are unfortunately beyond the scope of the present work.

    
    
    
    

    
    
    
    
    
    
    
    

    Finally, we would like to address the strong magnetic field recovered that drives the accretion process on EX Lup.
    With the B_ℓ reaching 4.2 kG in the accretion shock, a dipole strength of 1.9 kG within a large-scale field of 1 kG reaching 4.8 kG locally, and a small-scale field exceeding 3 kG in the optical domain, EX Lup has one of the strongest magnetic field among cTTSs, and the strongest among those with dipole-dominated topology.
    Such a configuration can set the suitable condition invoked by <cit.> for their hypothesis of an episodic accretion due to the magnetospheric accretion process itself.
    Indeed, with such a dipolar strength, the magnetospheric radius is outside but close to the corotation radius (see Table <ref>).
    In such a situation, instability arises due to the fact that angular momentum is transferred from the star to the disc (the so-called propeller regime), without being enough to drive an outflow.
    This magnetic interaction only prevents accretion, pilling up the gas in the inner disc and increasing its pressure, thus forcing the inner edge of the disc to move inward and cross the corotation radius allowing the accretion to occur.
    Once the gas reservoir has been accreted, the inner edge of the disc moves outward and another cycle starts <cit.>.
    This phenomenon is different from the magnetospheric inflation reported for other cTTSs <cit.>, where a significant difference between the truncation and corotation radii induces a torsion of the magnetic field lines producing a toroidal field that inflates the whole magnetosphere. 
    This inflation goes up to an opening of the magnetic field lines that produces a magnetospheric ejection <cit.> before reconnecting to the disc.
    Such a phenomenon is recurrent and produces an accretion variability as well.
    Still, the time scale of such a cycle is of the order of the stellar rotation period, far shorter than the gas pilling up invoked by <cit.>, and no signature of magnetospheric inflation, nor ejection, were detected in EX Lup.
    However, we would like to stress to the reader that no accumulation of matter at the magnetospheric radius was detected either, EX Lup seems only to be in the same initial conditions as the latter authors' theory. 
    In addition, the BIS suggest another source for radial velocity variation, such as a companion. 
    Tidal interactions <cit.> or thermal instabilities in the disc <cit.> induced by a companion are also existing hypotheses for EXor behaviour. 
    Moreover, recent works by <cit.> favour the latter instability as the origin of the episodic accretion of FUor objects.
    Finally, the disc itself is not investigated in this work, instabilities in the disc are thus hypotheses that cannot be excluded <cit.>.
    
    





§ CONCLUSIONS


    
    EX Lup is the prototypical EXor-type object whose recurrent bursts and outbursts were previously studied in detail using spectroscopy and photometry.
    However, until this work, no information about its magnetic field was derived, despite its key role in the accretion process of cTTS and its possible origin for episodic accretion.
    Here we provide the first spectropolarimetric time series, over two epochs (2016 and 2019) and two wavelength domains (optical and infrared), study of EX Lup, the first EXor whose magnetic field is studied.

    We confirmed an ongoing magnetospheric accretion process as seen on many cTTSs.
    It is represented by an accretion funnel flow and an accretion shock corotating with the stellar surface and driven by a kG dipolar magnetic field.
    The funnel flow seems aligned with the accretion shock, itself located near the magnetic dipole pole, which is consistent with the stable pattern observed between the epochs studied.

    The magnetic field of EX Lup has shown some disparities between wavelength domains, much weaker in the infrared.
    This can be understood as a wavelength dependency of the parameters studied, but can also point to a vertical structure of the magnetic field, as different wavelengths are tracing different heights in the photosphere.
    An expected small-to-large scale effect is also observed, by the different field strengths recovered using ZDI and Zeeman intensification, but also by the opposite polarity of the field associated with the accretion shock and the one recovered using LSD. 
    This indicates a very small accretion shock, as expected from the low mass accretion rate in quiescence and the large radial velocity variation of the emitting region.

    Finally, the multi-kG field recovered for EX Lup is pointing to a magnetospheric radius being close but outside the corotation radius. 
    This configuration is suitable for disc instabilities induced by the magnetic field that yield accretion cycles.
    These cycles might explain the accretion bursts observed on EX Lup, suggesting an inherently episodic magnetospheric accretion process.
    However, a definite identification of the origin of EXor behaviour is beyond the scope of this paper, and other theories implying the disc itself or a companion cannot be excluded.

    



    We would like to warmly thanks Oleg Kochukhov for useful discussions about the magnetic field of EX Lup, as well as Colin P. Folsom for his help in using the new version of his  package.
    
    The  package is available at <https://github.com/folsomcp/specpolFlow>.

    The  package is available at <https://github.com/pouillyk/PySTELLA>
    
    This research was funded in whole or in part by the Swiss National Science Foundation (SNSF), grant number 217195 (SIMBA). For the purpose of Open Access, a CC BY public copyright licence is applied to any Author Accepted Manuscript (AAM) version arising from this submission.
    
    Based on observations obtained at the Canada–France– Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada–France–Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site.

    This work was also supported by the NKFIH excellence grant TKP2021-NKTA-64.



aa






    

§ CROSS-CORRELATION MATRICES OF OPTICAL EMISSION LINES

    
        Here we present the cross-correlation matrices of the ESPaDOnS emission lines studied in this work.
        These matrices are discussed in Sect. <ref>.
        
        

        
    

    

§ STOKES I AND V FIT FROM ZDI RECONSTRUCTION

    
        In this section, we show the fit of the LSD profiles by the ZDI reconstruction for the three data sets studied.
        The Stokes I profiles are shown in Fig. <ref> and the Stokes V profiles in Fig. <ref>.
        
    
        

    

§ BIS VS. RADIAL VELOCITY

    
        This appendix presents the analysis of the BIS and radial velocity correlation. 
        These results are discussed in Sect. <ref>.
        




















    






        






        



