SPIRou reveals unusually strong magnetic fields of slowly rotating M dwarfs

In this paper, we study six slowly rotating mid-to-late M~dwarfs (rotation period $P_{\mathrm{rot}} \approx 40-190\,\mathrm{dy}$) by analysing spectropolarimetric data collected with SPIRou at the Canada-France-Hawaii Telescope as part of the SPIRou Legacy Survey from 2019 to 2022. From $\approx$100--200 Least-Squares-Deconvolved (LSD) profiles of circularly polarised spectra of each star, we confirm the stellar rotation periods of the six M~dwarfs and explore their large-scale magnetic field topology and its evolution with time using both the method based on Principal Component Analysis (PCA) proposed recently and Zeeman-Doppler Imaging. All M~dwarfs show large-scale field variations on the time-scale of their rotation periods, directly seen from the circularly polarised LSD profiles using the PCA method. We detect a magnetic polarity reversal for the fully-convective M~dwarf GJ~1151, and a possible inversion in progress for Gl~905. The four fully-convective M~dwarfs of our small sample (Gl~905, GJ~1289, GJ~1151, GJ~1286) show a larger amount of temporal variations (mainly in field strength and axisymmetry) than the two partly-convective ones (Gl~617B, Gl~408). Surprisingly, the six M~dwarfs show large-scale field strengths in the range between 20 to 200\,G similar to those of M~dwarfs rotating significantly faster. Our findings imply that the large-scale fields of very slowly rotating M~dwarfs are likely generated through dynamo processes operating in a different regime than those of the faster rotators that have been magnetically characterized so far.


INTRODUCTION
M dwarfs are known to host strong magnetic fields with large-and small-scale field strengths that may exceed 1 kG (Morin et al. 2010;Kochukhov 2021).Zeeman-Doppler-Imaging (ZDI, Donati & Brown 1997;Donati et al. 2006) revealed different types of large-scale field topologies for M dwarfs: the partly convective early M dwarfs usually showing more complex, relatively weaker fields with nonaxisymmetric poloidal fields and significant toroidal components (Donati et al. 2008).Mid M dwarfs often display simpler and stronger, mainly poloidal and axisymmetric large-scale fields (Morin et al. 2008) whereas the fully-convective late M dwarfs of lowest masses may end up showing large-scale fields in either configuration, (Morin et al. 2010, see also the review by Donati & Landstreet 2009).
Besides, magnetic activity (diagnosed by various proxies) increases for shorter rotation periods until it saturates, i.e. no longer increases with decreasing rotation periods (see e.g.Saar 1996;Wright et al. 2011).In the unsaturated regime, both large-and small-scale fields, diagnosed by polarised and unpolarised Zeeman signatures on ★ E-mail: lisa.lehmann@irap.omp.euline profiles, increase with decreasing rotation periods (Vidotto et al. 2014;See et al. 2015;Reiners et al. 2022).
The main parameter that describes magnetic fields and activity of most M dwarfs is found to be the Rossby number , equal to the rotation period divided by the convective turnover time (e.g.Noyes et al. 1984;Wright et al. 2018 and references therein), with magnetic fields and activity increasing with decreasing  until saturation occurs at  ∼ 0.1 and below.Whereas large-scale fields of M dwarfs featuring  < 1 have been extensively studied, very little is known about those of very slow rotators with  ∼ 1 or larger.
In this paper, we explore large-scale fields of a small sample of very slowly rotating M dwarfs, whose fields and rotation periods were inaccessible to optical instruments.The six M dwarfs were observed with the SpectroPolarimetre InfraRouge (SPIRou), the nearinfrared spectropolarimeter and velocimeter recently mounted on the Canada-France-Hawaii Telescope (CFHT), in the framework of the SPIRou Legacy Survey (SLS, Donati et al. 2020).The SLS is a Large Programme carried out with SPIRou at CFHT from early 2019 to mid 2022, with a 310-night time allocation spread over this period.The two main goals of the SLS are (i) the search for habitable Earth-like planets around very-low-mass stars and (ii) the study of low-mass star and planet formation in the presence of magnetic fields.Its long timeframe (of 7 semesters) enables us to investigate the temporal variability in the time series of the monitored M dwarfs, and in particular to independently estimate rotation periods of up to a few hundred days and to study the temporal evolution of their large-scale magnetic fields, (Fouqué et al. 2023;Bellotti et al. 2023;Donati et al. 2023, hereinafter D23).
In Section 2 we will present the details of the observations and targets.To analyse the magnetic field properties and to determine the stellar rotation period, we use different methods explained in Section 3, before we present our results for each M dwarf individually in Sec.4-9.We conclude and discuss our results in Section 10.

SPIROU OBSERVATIONS
We analyse here a total of 986 circularly polarised spectra collected with SPIRou.The spectra span a wavelength range of 0.95-2.5m in the near-infrared with a resolving power of  = 70 000.Further details about SPIRou and its spectropolarimetric capabilities can be found in Donati et al. (2020).To process the data, we used the new version of Libre ESpRIT, i.e., the nominal reduction pipeline of ES-PaDOnS at CFHT optimised for spectropolarimetry and specifically adapted for SPIRou (Donati et al. 2020).
We applied Least-Squares Deconvolution (LSD, Donati et al. 1997) to all reduced unpolarised (Stokes ) and circularly-polarised (Stokes ) spectra using a M3 mask constructed from outputs of the VALD-3 database (Ryabchikova et al. 2015) assuming a temperature  eff = 3500 K, a logarithmic surface gravity log  = 5 and a solar metallicity [M/H].Although our 6 stars do not have the exact same atmospheric properties (see Table 1), we nonetheless used a single mask, whose impact on the LSD results is only marginal, especially in the near-infrared domain where the synthetic spectra only provide a rough match to observed ones (e.g., Cristofari et al. 2022).Besides, the mask we chose corresponds to the coolest atmospheric model available by default on the VALD-3 data base.We have selected the atomic lines with a relative depth greater than 10% and resulting in 575 lines for the mask.For further details see D23, Sec. 2. The ephemeris used to calculate the phase and rotation cycle in this paper are summarised in Tab.A1 for all targets of our sample.
The six M dwarfs studied in this paper are part of the 43 star sample analysed by Fouqué et al. (2023) and D23.These two papers aimed at determining, whenever possible, the rotation periods of the sample targets, by applying quasi-periodic Gaussain-Process-Regression (GPR) to times series of their longitudinal fields  ℓ , i.e., the line-of-sight-projected component of the vector magnetic field averaged over the visible stellar hemisphere.All six stars of our sample have well identified rotation periods according to D23, whereas Fouqué et al. (2023), using data reduced with the nominal SPIRou pipeline APERO (optimised for RV precision, Cook et al. 2022) was able to derive rotation periods for four of them (consistent with those of D23).
The key stellar parameters of our sample are presented in Table 1 and are mostly extracted from Cristofari et al. (2022) who studied the fundamental parameters of these stars by analysing their SPIRou spectra.

MODEL DESCRIPTION
To analyse the magnetic properties of our six slowly rotating M dwarfs, we use both the method based on Principal Component Analysis (PCA) recently proposed by Lehmann & Donati (2022), as well as Zeeman-Doppler Imaging (ZDI, Donati & Brown 1997;Donati et al. 2006), applied to our set of LSD Stokes  profiles.Lehmann & Donati (2022) proposed a method to retrieve key information about the large-scale stellar magnetic field directly from time series of Stokes  LSD profiles without the need of an elaborate model of the field topology or several stellar parameters (e.g., the projected equatorial velocity   sin  or the inclination of the stellar rotation axis ).The method provides information about the axisymmetry, the poloidal / toroidal fraction of the axisymmetric component, the field complexity and their evolution with time.

PCA analysis of the LSD Stokes 𝑉 profiles
One first determines the mean profile of the whole Stokes  time series, which stores information about the axisymmetric component of the large-scale field.In contrast to Lehmann & Donati (2022), we use the weighted mean profile providing better results for time series such as ours, where all LSD profiles do not have the same SNR.The averaged Stokes  LSD profile can be decomposed into its antisymmetric component (with respect to the line centre), which is related to the poloidal component of the axisymmetric large-scale field, and its symmetric component (with respect to the line centre), which is probing the toroidal component of the axisymmetric largescale field (see e.g.Fig. 1a).
To evaluate the non-axisymmetric field, we subtract the weighted mean profile (taken over all seasons) from the Stokes  time series removing the signal of axisymmetric field.The resulting meansubtracted Stokes  profiles store now the information about the non-axisymmetric component of the large-scale field and are analysed using a weighted PCA (Delchambre 2015) returning eigenvectors and coefficients.In the weighted PCA, the Stokes  profiles are weighted by the squares of their SNRs, taking into account the different noise levels.Thus, for the long time series analysed in this paper, with uneven SNR over the 7 semesters, the weighted PCA gives better results than a classical PCA where all profiles are treated equally.The PCA coefficients, and in particular their fluctuation with time, can reveal the complexity of the large-scale field and its long-term temporal evolution.
Given the long time range of the SLS data, we further split the Stokes  time series at successive observing seasons into 2-3 seasons per star.To evaluate the evolution of the axisymmetric field from season to season, we determine the weighted mean profiles per season and compare them to one another (see e.g.Fig. 1c left column).To study the evolution of the non-axisymmetric field, we compare the coefficients of the different seasons (see e.g.Fig. 1c middle and right column).We caution that the coefficients are derived from the weighted PCA of the mean-subtracted Stokes  time series using the weighted mean profile computed across all seasons (e.g.Fig. 1a) and not the weighted mean profile of each individual season (e.g.Fig. 1c left column).The usage of the weighted mean Stokes  profiles of each individual season would prevent a direct comparison of the different seasons.For example, it would centre the coefficients for each season, so that we will lose the information if the nonaxisymmetric field becomes more or less positive / negative from one season to another, and also the amplitudes of the coefficients are no longer comparable.Further information about the PCA method can be found in Lehmann & Donati (2022).
In addition, Lehmann & Donati (2022) showed that the sensitivity of the PCA method for toroidal fields decreases for low   sin .As all our stars have   sin  ≤ 0.5 km s −1 , we are likely to miss large-scale toroidal fields.We provide a typical 1 error bar for the axisymmetric toroidal field for each star and each observing season of our sample.Cristofari et al. (2022) including spectral type, effective temperature  eff , stellar mass, stellar radius  ★ , metallicity [M/H], log .The rotation period  rot is copied from D23.For the Rossby number  =  rot /, we use  rot (8th column) and determine the convective turnover time  via Wright et al. (2018, Eq. 6).In the last column, we give the projected equatorial velocity   sin  = 2  ★ rot sin , determined from  ★ (5th column),  rot (8th column) and an assumed inclination of  = 60 • between the stellar rotation axis and the line-of-sight.

Gaussian Process modelling of the time series
We analyse the temporal evolution of the M dwarf's topology with the help of the PCA determined coefficients of the mean-substracted Stokes  time series.For our slowly rotating stars, most of the time only the first eigenvector and therefore only the first coefficient shows a signal.To directly compare the temporal evolution of the coefficients with the result from the longitudinal field  ℓ (presented by D23) for the individual stars, we scale and translate the first coefficient, which we call  1 , using a linear model (scaling factor and offset), that minimises the distance between the first coefficients and  ℓ taking into account the measurement errors on  ℓ .We can re-determine the stellar rotation period  rot of our six M dwarfs using a quasi-periodic (QP) GPR fit to  1 allowing us a direct comparison with the QP GPR results of  ℓ presented by D23.
In contrast to D23, we use the python model presented by Martioli et al. (2022) based on george (Ambikasaran et al. 2015).Our adapted covariance function (or kernel) is given by where    =   −   is the time difference between the observations  and ,  is the amplitude of the Gaussian Process (GP),  is the decay time describing the typical time-scale on which the modulation pattern evolves,  is the smoothing factor indicating the harmonic complexity of the QP modulation (lower values indicating higher harmonic complexity) and  rot is our new estimate of the stellar rotation period.The GP model parameters are fitted by maximising the following likelihood function L using the python package scipy.optimize: where K is the QP kernel covariance matrix,  the diagonal variance matrix of  1 , S the diagonal matrix  2 I (with  an added amount of uncorrelated white noise (Angus et al. 2018) and I the identity matrix),  the number of observations and  corresponds to  1 .The posterior distribution of the free parameters is sampled using a Bayesian Markov Chain Monte Carlo (MCMC) framework applying the package emcee (Foreman-Mackey et al. 2013).For the MCMC, we use 50 walkers, 200 burn-in samples and 1000 samples.Tab. 2 provides a summary of the results for all  1 GPR fits in this study.
For three stars, the decay time  was fixed as in D23 (see Tab. 2).Information about the assumed prior distributions and the posterior distributions for each parameter and each GPR fit can be found in appendix A. Furthermore, we applied the above GP model to the  ℓ values of D23 (see appendix B).This allows a direct comparison of the GP results for the same  ℓ dataset with our GP routine and the GP routine used by D23, as well as a comparison of the GP results for  1 and  ℓ obtained by the same GP routine.In general, we find that  1 has lower RMS, often shows lower  values and provides smaller errors for  rot when the topology is not highly axisymmetric.

Zeeman-Doppler-Imaging
We determined the large-scale vector magnetic field at the surface of the six M dwarfs, for each season, using ZDI.ZDI iteratively builds up the large-scale magnetic field and compares the synthetic Stokes profiles corresponding to the current magnetic map, assuming solid body rotation, with the observed Stokes profiles until it converges on the requested reduced chi-square value  2  between the observed and synthetic data.The problem being ill-posed, i.e., with an infinite number of solutions featuring the requested agreement to the data, ZDI chooses the one with maximum entropy, (i.e., minimum information in our case, Skilling & Bryan 1984).The surface magnetic field is described with a spherical harmonics expansion given in Donati et al. (2006), where the  ℓ, and  ℓ, coefficients of the poloidal component are modified as indicated in Lehmann & Donati (2022, Eq. B1).To compute synthetic Stokes profiles, the stellar surface is decomposed into a grid of 1000 cells.For each cell the local Stokes  and  profiles are determined using Unno-Rachkovsky's analytical solution to the equations of the polarised radiative transfer in a plane-parallel Milne-Eddington atmosphere (Landi Degl 'Innocenti & Landolfi 2004).The Stokes profiles are integrated over the visible hemisphere for each observing phase applying a mean wavelength of 1700 nm and a Landé factor of 1.2.
For the slowly rotating M dwarfs of our sample, we see no obvious variations in the Stokes  LSD profiles beyond those attributable to radial velocity variations, so that we only use the Stokes  profiles for determining the magnetic field map via ZDI.Nevertheless, we make sure that the synthetic Stokes  profiles computed with ZDI agree well with the averaged observed Stokes  profile, especially in terms of width and depth.We assumed a fraction   of each grid cell, which actually contributes to the Stokes  profile.This fraction   is called filling factor of the large-scale field and is set equally to all cells, see also Morin et al. (2010); Kochukhov (2021).For each star, we set   = 0.1 motivated by the results of Klein et al. (2021) for the slow rotator Proxima Centauri and the results of Moutou et al. (2017) for the SPIRou sample.We confirm the choice of   = 0.1 by finding lower  2  values with   = 0.1 compared to   = 1 for each season of the different M dwarfs.The filling factor for the Stokes  profiles is set to   = 1.0 in consistency with the literature, (Morin et al. 2010;Kochukhov 2021).The   sin  of our sample is ≤ 0.5 km s −1 (see Tab. 1) and prevents us from reliably determining Table 2. Summary of the best-fitting parameters of the QP GPR fits applied to  1 for the six M dwarfs in our sample, where rms is the root-mean-square of the residuals and  2  the reduced chi-square value of the GPR fit.Fixed parameters are shown in italics.A comparison with the results of the GPR fit to the  ℓ data, and with those of D23, is given in Table B1.star rotation period decay time smoothing factor amplitude white noise rms Gl 905 111.7 the inclination of the stellar rotation axis for the M dwarfs, so that we set the inclination to 60 • for all M dwarfs.This is motivated by the steep modulation patterns seen for  ℓ and  1 for most targets, that can not be obtained for pole-on viewed stars.Another reason is that higher values of  are intrinsically more likely than smaller ones.We restrict the spherical harmonics of the ZDI reconstructions to ℓ = 7, as we see little magnetic energy stored in ℓ ∼ 5 − 7.
For our analysis we use 219 Stokes  LSD profiles observed with SPIRou between 2019 Apr and 2022 June and split the data in three seasons (2019 Apr -Dec, 2020 May -2021 Jan, 2021 June -2022 Jan) for the per-season analysis.The 15 profiles collected in 2022 June at the beginning of a new season, only covering 9% of a rotation cycle, were left out of the per-season analysis.

PCA analysis of Gl 905
First, we investigate the large-scale field topology using PCA (Lehmann & Donati 2022).The weighted mean profile of all Stokes  profiles is antisymmetric with respect to the line centre indicating a poloidal axisymmetric large-scale field (see Fig. 1a).The symmetric component of Gl 905's mean profile exceeds the noise level ( 2  = 1.4) but is likely due to an uneven phase coverage in the season 2020/21 (see Sec. 4.2).In Fig. 1b, we show the first two eigenvectors of the mean-subtracted Stokes  profiles allowing the analysis of the non-axisymmetric field.Only the first eigenvector shows an antisymmetric signal with respect to the line centre.All other eigenvectors display noise.
In Fig. 2 top, we plot the temporal evolution of  1 , which appears very similar to the temporal evolution of  ℓ determined with our GP model (see Fig. 2 bottom) and to D23's results (see Fig. A12 middle in D23).When only one eigenvector is significant, as it is the case here for Gl 905, the  1 curve mimics that of  ℓ (Lehmann & Donati 2022).
Our QP GPR model of  1 finds a rotation period of  rot = 111.In Fig. 1c, we plot the mean profiles (left column) and the phasefolded coefficient curves colour-coded by rotation phase (middle and right columns) for the three seasons (one season per row).They exhibit large changes in the large-scale field topology from season to season, allowing us to draw first conclusions about the field evolution of Gl 905.We recall that the coefficients for all three seasons are computed from the mean-subtracted Stokes  profiles, using the weighted mean derived from the full data set, and not from the profiles of each season.The same applies to the five M dwarfs discussed in Secs.5-9.
The mean profiles of the first two seasons (2019 and 2020/21) are antisymmetric with respect to the line centre and indicate a mostly poloidal axisymmetric field although the symmetric component is larger for 2020/21 (see Fig. 1c).This may reflect an increasing toroidal field but is more likely due to the uneven phase coverage of this season (with more than 75% of the observations concentrating between phase 0.3 and 0.75).
For the first season 2019,  1 features a roughly sinusoidal behaviour indicating a mainly dipolar configuration.For 2020/21,  1 appears more complex than for 2019 implying that the field becomes more complex, too.For the last season 2021/22, the topology changes more drastically: the mean profile is close to zero indicating a much lower axisymmetric component than before.The phase at which  1 reaches its maximum is shifted, with  1 being more positive now, while it is mainly negative before.Considering the sign of the mean profile and the eigenvector, this suggest that the main polarity of the large-scale field is evolving from a predominantly negative polarity to a positive polarity.We can conclude from the PCA analysis, that the large-scale field topology becomes more complex from 2019 to 2020/21 before it becomes mostly non-axisymmetric and possibly initiates a polarity reversal.

ZDI reconstructions of Gl 905
We conclude our analysis by deriving vector magnetic field maps for Gl 905 using ZDI, for each of the three main observing seasons.The maps are shown in Fig. 3 and their magnetic properties are summarised in Table 3.We were able to fit all three ZDI maps down to  2  ≈ 1.0 assuming The ZDI maps confirm the conclusions we derived from the PCA analysis.The topology gets indeed more complex from 2019 to 2020/21 and the degree of axisymmetry decreases from around 70% to 4% for the last season 2021/22.The surface mean magnetic field decreases from 128 G to 64 G.Most prominent is the hint of an ongoing polarity reversal from negative to positive radial field taking place in the last season.
To test whether the symmetric component of the mean profile for season 2020/21 indeed results from an uneven phase coverage, we simulate 24 evenly phased Stokes  LSD profiles from the 2020/21 ZDI map (see Fig. 3 middle column) and determine the corresponding mean profile and its symmetric and antisymmetric components (see Fig. C1).The symmetric component disappears with even phase sampling, confirming that the mean LSD profile provides no observational hint for a large-scale axisymmetric toroidal field at the surface of Gl 905.
The reconstructed surface averaged toroidal field ⟨ tor ⟩ is lower than 10 G for each season.To derive a 1 error bar on the simplest possible large-scale axisymmetric toroidal field (described with spherical harmonics coefficients ℓ = 1 and  = 0) at the stellar surface, we proceed in the following way: (1) artificially add an axisymmetric toroidal field of strength ⟨ tor ⟩ to the reconstructed ZDI map, (2) simulate the corresponding Stokes  LSD profiles with the phase coverage and SNR of the actual observations and (3) compute the new  2 with respect to the observed LSD profiles and raise ⟨ tor ⟩ until  2 is increased by 1 with respect to the optimal ZDI fit.We find 1 error bars ranging from 180 G in 2019 down to 55 G in 2020/21.

GJ 1289
The fully convective M dwarf GJ 1289 ( = 0.21 ± 0.02 M ⊙ , Cristofari et al. 2022) is the next star in our sample.SPIRou observed GJ 1289 from 2019 Sept until 2022 June providing a time series of 204 LSD profiles split into three seasons (2019 June -Dec, 2020 May -2021 Jan, 2021 June -2022 Jan) for the per-season analysis.As for Gl 905, the 14 profiles collected in 2022 June at the beginning of a new season, only covering 16% of a rotation cycle, were left out of the per-season analysis.

PCA analysis of GJ 1289
The mean profile is perfectly antisymmetric with respect to the line centre indicating a dominant axisymmetric poloidal component (see Fig. 4a).Both the first and second eigenvectors are found to be antisymmetric with respect to the line centre, which is a strong hint of a non-axisymmetric poloidal component (see Fig. 4b).All further eigenvectors trace noise.
The QP GPR fit applied to  1 (see Fig.In Fig. 4c, we show the mean profile and the phase-folded coefficients split by season.Comparing the mean profiles of the three seasons, the axisymmetric component grows in amplitude and stays always poloidal.As the amplitude of the coefficients increases as well, the magnetic field becomes in general stronger.
The phase-folded coefficient curves indicate a rapidly evolving and complex large-scale field, as we see variations from one rotation cycle to the next (see season 2020/21) and trends that are more complex than sine waves (e.g. for 2019 and 2020/21).Furthermore, Table 3. Magnetic properties of Gl 905 extracted from the ZDI maps per season: the start and end month of the observations used for the ZDI maps, the surface averaged unsigned magnetic field ⟨  ⟩[G], the surface average unsigned dipole magnetic field ⟨ dip ⟩[G], the typical 1 error bar on the ZDI reconstructed surface averaged toroidal field  ⟨tor⟩ , the fractional energy of the poloidal field, the fractional energy of the axisymmetric field (only  = 0 modes), the fractional energy of the dipole component, the tilt angle of the dipole (ℓ = 1) with respect to the negative pole, the phase at which the dipole field faces the observer, the reduced  2 values for the Stokes  profiles ( 2  , ; corresponding to a B = 0 G fit), for the ZDI fit of the Stokes  profiles ( 2  ,,ZDI ) and for the Null profiles (   season 2020/21 stands out, with the second eigenvector contributing significantly to the Stokes  signal before disappearing again for the last season 2021/22, for which  1 shows a simpler trend.

ZDI reconstructions of GJ 1289
We were able to fit the Stokes  profiles for all seasons down to  2  ≈ 1.0 assuming  rot = 73.66d,   sin  = 0.14 km s −1 ,  = 60 • ,   = 0.1.In the first season, the data set only includes about half of the observations of those from the two other seasons and the achieved  2  , = 1.46 is several times lower than the  2  , of the following seasons (see Tab. 4).In the first season, ZDI reveals a weak marginally complex field topology.In 2020/21, the surface averaged field ⟨  ⟩ becomes twice as large due to a growing axisymmetric poloidal dipole and ZDI reconstructs a more complex azimuthal field, featuring a quadrupolar non-axisymmetric azimuthal structure (see Fig. 6).In the last season 2021/22, the dipole tilts more strongly to 39 • and dominates the field topology (see Tab. 4).The toroidal field of the ZDI maps varies between 7 − 25 G, which is again lower than the typical 1 error bar that we derive ranging from 40 to 100 G.

PCA analysis of GJ 1151
The mean profile is close to zero indicating a strongly nonaxisymmetric topology (see Fig. 7a).Only the first eigenvector of the mean-subtracted Stokes  profiles significantly differs from the noise and features an antisymmetric signal with respect to the line centre (see Fig. 7b).
The QP GPR model fits  1 down to a  2  = 0.99 (see Fig. 8 top).We fix the decay time to 300 d similar to D23 and find a  rot = The mean profile of the first two seasons is antisymmetric with respect to the line centre (axisymmetric poloidal field) and is relatively weak (see Fig. 7c).The coefficient  1 shows no obvious trend with phase for the first season 2019/20 and just start to display a weak variation with phase for 2020/21.The low amplitude of the mean profiles and coefficients indicate that the magnetic field must be very weak during the first two seasons.For the last season 2021/22, the amplitude of the mean profile is twice as high as before and also  1 shows a higher amplitude ( 2  = 7.9), indicating that the magnetic field increases significantly for 2021/22.We also notice that the sign of the mean profile (and therefore the projected main polarity of the large-scale magnetic field) changed from negative to positive for the last season, hence why the mean profile over the whole time series is close to zero (see Fig. 7a).
For the first two seasons, the Stokes  profiles are weak (see Fig. C6) and so is the reconstructed field, with a dominant negative polarity in the upper hemisphere that is consistent with the corresponding mean profiles.We see a small increase of ⟨  ⟩ in the second season of 2020/21 explaining the higher amplitude of  1 seen in the PCA analysis (see Fig. 7c).For the last season 2021/22, the ZDI map shows a strongly tilted dipole (tilt angle = 55 • ), that flipped  polarity, and the surface averaged field is more than twice as high as before (⟨  ⟩ = 63 G).To the best of our knowledge this is the first polarity reversal seen in the vector magnetic field map of an M dwarf.
The reconstructed toroidal field ⟨ tor ⟩ is lower than 8 G for all seasons of GJ 1151 and we find that the 1 error bar on the toroidal field ranges between 370 and 450 G.

PCA analysis of GJ 1286
The mean profile of GJ 1286 is antisymmetric with respect to the line centre and appears more noisy than usual but clearly indicates a purely axisymmetric poloidal field (see Fig. 10a).The first eigenvector has an antisymmetric shape, too, and is the only one that emerges from the noise (see Fig. 10b).
The best-fitting model of the QP GPR for  1 finds a  rot = 186.8+9.−13 d).The mean profile of season 2020 is nearly twice as high as for 2021 (see Fig. 10c).Also  1 shows a lower amplitude for 2021.We can therefore conclude, that the surface averaged field decreases for 2021.The field topology becomes simpler as the  1 curve gets less complex with phase for 2021.
As concluded from the PCA analysis, the ZDI maps confirm that the topology becomes simpler and weaker: the fractional energy of the dipole component  dip increases from 0.73 to 0.79, while ⟨  ⟩ decreases by almost half from 113 G to 71 G (see Fig. 12 and Tab. 6).The reconstructed toroidal field ⟨ tor ⟩ of the seasons are 9 G and 6 G, respectively, while the typical 1 error bars on the axisymmetric toroidal field are 325 G and 300 G, respectively.

GL 617B
Gl 617B (EW Dra, HIP 79762, LHS 3176) is a partly convective M dwarf with  = 0.45 ± 0.02 M ⊙ (Cristofari et al. 2022) and was observed between 2019 Sept and 2022 June with SPIRou.Our following analysis is based on 144 LSD Stokes spectra, which we split into three seasons (2020 Feb -Oct, 2021 Jan -July, 2022 Mar -June) for the per-season analysis.As for the other stars, the first 15 spectra collected in 2019 were left out of the per-season analysis.

PCA analysis of Gl 617B
We find, that the mean profile is much larger than the mean-subtracted Stokes  profiles indicating that the axisymmetric component of the magnetic field is dominant (see Fig. 13a).The mean profile is again antisymmetric with respect to the line centre, indicating an axisymmetric poloidal field.The first eigenvector is already very noisy and is the only one that shows a clear signal, confirming that the field is indeed dominantly axisymmetric (see Fig. 13b).
The QP GPR model fits  1 down to a  2  = 0.66 finding a rotation period of 37.8 +8.The mean profile is antisymmetric to the line centre and therefore poloidal dominated for all three seasons, but varies in amplitude (see Fig. 13c, left column).Nonetheless,  1 traces a varying nonaxisymmetric component.Season 2020 shows the highest range in amplitude of  1 , indicating the largest dipole tilt angle of all three seasons, although it will still be small (< 20 • ) due to the predominantly axisymmetric topology of Gl 617B.

ZDI reconstructions of Gl 617B
We could fit Gl 617B down to  2  ≈ 1.0 assuming  rot = 40.4d,   = 0.1,   sin  = 0.50 km s −1 ,  = 60 • .The ZDI maps confirm a very axisymmetric, poloidal configuration.The axisymmetry is always equal to or greater than 97%, so variations of the non-axisymmetric field are difficult to see, but appear largest in 2021 (see Tab. 7).The data set in season 2020 shows the largest tilt angle (7 • ) as predicted by the PCA analysis.The reconstructed toroidal field ⟨ tor ⟩ reaches 9 G for 2020, 6 G for 2021 and 2 G for 2022, while the estimated 1 error bar on the axisymmetric toroidal field is about 7 G in 2020, 13 G in 2021 and 6 G in 2022, which is a lower uncertainty than for the other M dwarfs thanks to the higher   sin  = 0.5 km s −1 of Gl 617B.
We see ⟨  ⟩ changing by approximately ±25 G for the three seasons, otherwise the main properties of the maps are similar (see Tab. 7).
For highly axisymmetric topologies, it is difficult to infer the inclination .It may be that  is actually lower for Gl 617B.We provide the ZDI maps for an inclination  = 30 • and   sin  = 0.29 km s −1 , while otherwise using the same parameters (see Fig. C2).The  2   for the per-season analysis.As for the other stars, the first 17 spectra collected in early 2019 were left out of the per-season analysis.

PCA analysis of Gl 408
Gl 408 has the strongest mean profile compared to the meansubtracted profile in our sample (see Fig. 16a and C9).The first eigenvector of Gl 408 (the only one showing a signal) is already noisy, a strong indication of a very axisymmetric field topology (see Fig. 16b).Fig. 17 (top) presents the QP GPR fit of  1 , which gives a  rot = 175 +12 −14 d similar to D23 determining  rot = 171.0± 8.4 d.Fitting  ℓ with our GP routines, we derive a  rot = 170.7 +7.1 −9.8 d (see Fig. 17 bottom).The decay time was fixed at 200 d for both variables following D23.However, we find a decay time of ≈ 200 ± 70 d but higher  2 for GPR fits without fixing the decay time.The APERO reduced spectra of Gl 408 did not allow Fouqué et al. (2023) to determine a rotation period.
In Fig. 16c, we see that  1 is mostly flat for all three seasons, again indicating a highly axisymmetric topology.All mean profiles are antisymmetric with respect to the line centre and show an axisymmetric poloidal large-scale field.

ZDI reconstructions of Gl 408
All three seasons could be fitted down to  2 ≈ 1.0 assuming  rot = 171.0d,   = 0.1,   sin  = 0.10 km s −1 ,  = 60 • .The topology changes little over the three seasons and is characterised by a strong,  axisymmetric, poloidal dipole of negative polarity (see Fig. 18).It is the most stable topology in our sample and only ⟨  ⟩ varies marginally between 106 − 130 G (see Tab. 8).We find a 1 error bar on the axisymmetric toroidal field of 55 − 81 G for Gl 408 whereas the reconstructed ⟨ tor ⟩ ranges between 4 − 13 G.Similar to Gl 617B, we also determine the ZDI maps for an inclination of  = 30 • and   sin  = 0.06 km s −1 (see Fig. C3).The  2  values of the ZDI fits are again slightly higher for the lower inclination  = 30 • than for  = 60 • .

SUMMARY, DISCUSSION & CONCLUSIONS
In this paper, we study the large-scale magnetic field of six slowly rotating mid to late M dwarfs observed with SPIRou at the CFHT as part of the SLS from 2019 to 2022.The 3.5-yr time series, including ≈ 100 − 200 polarimetric spectra for each of our six M dwarfs, allowed us to confirm their rotation periods and to investigate their magnetic field topology using both our PCA analysis and ZDI.We use the reduced observations from D23 but different analysis tools to redetermine the rotation period.Our estimate of the rotation periods using  1 , i.e., the scaled and translated first coefficient of the PCA analysis (see Sec. 3.2), agrees with the results of D23 and Fouqué et al. (2023).We confirm that both Gl 617B and Gl 408, for which Fouqué et al. (2023) did not recover a rotation period, host very axisymmetric topologies with  axi ≥ 0.97 between 2019 and 2022.The higher the axisymmetry of the large-scale field, the smaller are the variations of  ℓ or  1 with time, and the harder it is to determine a rotation period.For the highly axisymmetric topologies, we find that the  2  of the GPR fits increases, reflecting that in such cases, the  ℓ curves are more sensitive to intrinsic variability and less to rotational modulation, and thereby reducing the ability at measuring rotation periods (see e.g.Fig. A5 or A6 and Tab.B1).
Using the PCA analysis, we derive information about axisymmetry and complexity directly from the LSD Stokes  time series, which are in agreement with the results obtained from the ZDI maps for all six M dwarfs, while PCA does not rely on any assumptions about stellar parameters such as   sin , inclinations, etc.
We find evidence for a polarity reversal of the large-scale field (via sign changes of  ℓ ,  1 or in the mean profiles) taking place on GJ 1151 and possibly also on Gl 905, for which the axisymmetric component collapsed during the last season (to be confirmed with new, ongoing, observations).For most stars, PCA traces the time-evolving field topologies using only the first eigenvector.For GJ 1289 we even detect two evolving field components directly from the Stokes  time series.This highlights that we are able to reliably detect topological complexity in the magnetic fields of slowly rotating M dwarfs directly from the observed LSD Stokes  profiles.The lower the   sin , the higher the 1 error bar on the toroidal field.The typical 1 error bars on the toroidal field ranges from 6 to 450 G depending on SNR and   sin .
We determined the ZDI maps for each season of our targets, obtaining a total of 17 vector magnetic field maps.The ZDI maps of GJ 1151 and Gl 905 confirm the polarity switches that were diagnosed with PCA, and further show that GJ 1151 may have been in a magnetically quiescent state until it became more magnetic in 2022, switching polarity at the same time.
The slowly rotating M dwarfs of our sample show large-scale field strengths in the range ⟨  ⟩ ≈ 20−200 G.They show similar ⟨  ⟩ to faster rotating M dwarfs in the saturated regime.We add our sample to  100 days.Besides, it implies a harsher interplanetary environment for potential close-in planets (e.g., Kavanagh et al. 2021).We stress that our paper focused on the most magnetic M dwarfs of the SLS sample (e.g., D23), whereas the other (less magnetic) stars of this sample will presumably be more in line with (and fill the gap between) the high- stars of the Vidotto et al. ( 2014) sample.This will be the subject of forthcoming studies.Beside, our results suggests that the large-scale fields of the very slowly rotating M dwarfs of our sample are likely generated through dynamo processes operating in a different regime than those of the faster rotators that have been magnetically characterized so far.Fig. 20 summarises the properties of the large-scale magnetic field topology for our six M dwarfs displaying all seasons on top of each other.It can be seen that the two partly convective M dwarfs (Gl 617B and Gl 408) show a smaller range of variations compared to the fully convective stars.The fully convective M dwarfs host large-scale fields that evolve on timescales comparable to their rotation periods.Our small sample suggests that fully-convective, slowly-rotating M dwarfs tend to have large-scale fields that are less axisymmetric than their more massive counterparts.
In conclusion, we have analysed six slowly rotating M dwarfs observed by the SLS over 3.5 years.We find, that the large-scale magnetic field of these M dwarfs is unusually strong despite their slow rotation (40-190 d) and suggest that the efficiency of the dynamo for mid and late M dwarfs depends on  in a different way than that reported in the literature for faster rotators.Furthermore, we find that the large-scale magnetic field topology of the fully convective M dwarfs exhibit a larger range of variations than those of the two partly convective targets of our sample.Given this, it may be useful in the future to apply the time-dependent ZDI (Finociety & Donati 2022), which has only been tested for faster rotating stars up to now.We detected a polarity reversal on one (GJ 1151) and possibly two (Gl 905) of the 4 fully-convective stars of our sample, suggesting that magnetic cycles may indeed be occurring in such stars, as initially suggested by Route (2016) from radio observations.Further longterm observations of the same type are needed to document in a more systematic fashion the long-term evolution of the large-scale magnetic fields of M dwarfs, and whether these field topologies are varying cyclically like for the Sun or in a more random fashion.MNRAS 000, 1-16 (2022) 7 +3.0 −3.2 d and decay time of  = 133 +18 −22 d very similar to the values derived by D23 ( rot = 114.3± 2.8 d and  = 129 +25 −21 d) and Fouqué et al. (2023) ( rot = 109.5+4.9 −5.4 d and  = 149 +26 −25 d) and also consistent with our GP fit of D23's  ℓ values ( rot = 114.4+3.5 −2.4 d and  = 130 +25 −32 d), see Tab. 2 and B1.For consistency, we use the rotation periods found by D23 to determine the rotation phase (see Tab. A1) and to model the ZDI maps for all six stars in our sample (see Tab. 1).

Figure 1 .
Figure 1.The PCA analysis for Gl 905.a.The mean profile (red) for all observations and its decomposition in the antisymmetric (blue dashed) and symmetric (yellow dotted) components (with respect to the line centre) related to the poloidal and toroidal axisymmetric field, respectively.This mean profile is used to determine the mean-subtracted Stokes  profiles to which we apply PCA, yielding the eigenvectors and coefficients shown in panels b and c. b.The first two eigenvectors of the mean-subtracted Stokes  profiles.c.The mean profile (left column),  1 (the scaled and translated first PCA coefficient introduced in Sec.3.2, middle column) and the coefficients of the second eigenvector (right column) for each season (one season per row).The mean profiles of the individual seasons are plotted in the same format as above.The coefficients are colour-coded by rotation cycle.

Figure 2 .
Figure 2. Temporal variations of  1 (top) and longitudinal field  ℓ (bottom) for Gl 905.We show the QP GPR fit and its 1 area as green shaded region in the top panel and the residuals in the bottom panel for both variables.The plot symbols are colour-coded by rotation cycle.The vertical grey lines separate the analysed seasons.The grey shaded region indicates a season for which not enough data were available for a reliable PCA and ZDI analysis.

Figure 3 .
Figure 3.The magnetic field maps of Gl 905 shown in a flattened polar view for the radial (top row), meridional (middle row) and azimuthal component (bottom row).In each plot, the visible north pole is in the centre, the thick line depicts the equator and the dashed line the latitudes in 30 • step.The ticks outside the plot illustrate the observing phases.The different seasons are shown next to each other (one season per column).The colour bar below the third row is used for all maps and indicates the magnetic field strength in G.The bottom panel summarises the main characteristics of the large-scale field of Gl 905 and its evolution with time.For each season, the symbol size indicates the magnetic energy, the symbol shape the fractional energy in the axisymmetric component and the symbol colour the fractional energy stored into the poloidal component of the field (see legend to the right).

Figure 9 .
Figure 9. Same as Fig. 3 for GJ 1151.Note, that the magnetic field switches polarity.
5 −2.6 d in agreement with the results of D23 ( rot = 40.4± 3.0 d) and our GP fit of  ℓ ( rot = 40.6 +2.1 −4.4 d, see Fig. 14).However, the decay time for the GP fit of  1 with  = 35 +8 −4 d is shorter than the results determined from the  ℓ curves ( = 69 +35 −23 d for D23 and  = 82 +45 −30 d from our own fit of  ℓ ).Fouqué et al. (2023) found no clear periodic variation using the APERO pipeline reduced spectra of Gl 617B.

Figure 19 .
Figure 19.The averaged unsigned magnetic field strength ⟨  ⟩ versus Rossby number  as shown by Vidotto et al. (2014) (in grey scales) including our sample of slowly rotating M dwarfs (red circles).Note that for this figure only, we determined  using Wright et al. (2011) for consistency with Vidotto et al. (2014).For further details and the coloured version of the original symbols and annotations see Fig. 4a of Vidotto et al. (2014).

Figure A1 .
Figure A1.The posterior density resulting from the MCMC sampling of the QP best-fitting GPR model for  1 (top) and  ℓ (bottom) for Gl 905.The concentric circle within each panel indicate the 1, 2 and 3 contours of the distribution.The blue lines mark the mode (maximum) of the posterior distribution and the black dashed lines the median and the 16% and 84% percentils of the posterior probability density function (PDF).

Figure C1 .
Figure C1.The mean profile and its decompositions obtained from 24 uniformly phased synthetic Stokes  LSD profiles of the 2020/21 ZDI map of Gl 905.The same format as in Fig. 1a is used.

Figure C2 .
Figure C2.The magnetic fields maps of Gl 617B using an inclination  = 30 • presented in the same format as in Fig. 3.

Figure C3 .
Figure C3.The magnetic fields maps of Gl 408 using an inclination  = 30 • presented in the same format as in Fig. 3.

Figure C4 .
Figure C4.The SPIRou observed LSD Stokes  profiles (black) and the ZDI fit (red) for Gl 905.The rotation phase is indicated to the right of the profile and the error to the left.Each panel corresponds to one season.

Table 1 .
The stellar characteristics of our sample are from 2  ,  ) and the number of observations per season (nb.obs).

Table 4 .
Same as Table3for GJ 1289.The tilt angle now refers to the positive pole.

Table 6 .
Same as Table 4 for GJ 1286.