# SPIRou reveals unusually strong magnetic fields of slowly rotating M dwarfs

L. Lehmann, J. -F. Donati, P. Fouque, C. Moutou, S. Bellotti, X. Delfosse, P. Petit, A. Carmona, Julien Morin, A. Vidotto, et al.

► To cite this version:

L. Lehmann, J. -F. Donati, P. Fouque, C. Moutou, S. Bellotti, et al. SPIRou reveals unusually strong magnetic fields of slowly rotating M dwarfs. Monthly Notices of the Royal Astronomical Society, 2024, 527, pp.4330-4352. 10.1093/mnras/stad3472. hal-04287657

HAL Id: hal-04287657

https://hal.science/hal-04287657

Submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# SPIRou reveals unusually strong magnetic fields of slowly rotating M dwarfs

L. Lehmann, J. -F. Donati, P. Fouque, C. Moutou, S. Bellotti, X. Delfosse, P. Petit, A. Carmona, Julien Morin, A. Vidotto, et al.

► To cite this version:

L. Lehmann, J. -F. Donati, P. Fouque, C. Moutou, S. Bellotti, et al. SPIRou reveals unusually strong magnetic fields of slowly rotating M dwarfs. Monthly Notices of the Royal Astronomical Society, 2024, 527, pp.4330-4352. 10.1093/mnras/stad3472. hal-04287657

HAL Id: hal-04287657

https://hal.science/hal-04287657v1

Submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# SPIRou reveals unusually strong magnetic fields of slowly rotating M dwarfs

L. T. Lehmann<sup>1,★</sup> J.-F. Donati<sup>1</sup> P. Fouqué,<sup>1</sup> C. Moutou,<sup>1</sup> S. Bellotti,<sup>1,2</sup> X. Delfosse,<sup>3</sup> P. Petit,<sup>1</sup> A. Carmona,<sup>3</sup> J. Morin<sup>1,4</sup> and A. A. Vidotto<sup>1,5</sup>

the SLS consortium

Accepted 2023 November 7. Received 2023 November 5; in original form 2023 August 7

## ABSTRACT

In this paper, we study six slowly rotating mid-to-late M dwarfs (rotation period $P_{\text{rot}} \approx 40-190$ d) by analysing spectropolarimetric data collected with SpectroPolarimetre InfraRouge (SPIRou) at the Canada-France-Hawaii Telescope as part of the SPIRou Legacy Survey from 2019 to 2022. From ≈100–200 least-squares-deconvolved (LSD) profiles of circularly polarized 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 polarized 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, and GJ 1286) show a larger amount of temporal variations (mainly in field strength and axisymmetry) than the two partly convective ones (Gl 617B and Gl 408). Surprisingly, the six M dwarfs show large-scale field strengths in the range between 20 and 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.

**Key words:** techniques: polarimetric – stars: imaging – stars: low-mass – stars: magnetic field – stars: rotation.

## 1 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 non-axisymmetric 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 polarized and unpolarized Zeeman signatures on line 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 $Ro$, 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 $Ro$ until saturation occurs at $Ro \sim 0.1$ and below. Whereas large-scale fields of M dwarfs featuring $Ro < 1$ have been extensively studied, very little is known about those of very slow rotators with $Ro \sim 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 near-infrared 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

© 2023 The Author(s).

<footnote n="1"><sup>1</sup>IRAP, Université de Toulouse, CNRS / UMR 5277, UPS-OMP, 14 Avenue E. Belin, F-31400 Toulouse, France</footnote>

<footnote n="2"><sup>2</sup>Science Division, Directorate of Science, European Space Research and Technology Centre (ESA/ESTEC), Keplerlaan 1, NL-2201 AZ Noordwijk, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4"><sup>4</sup>Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS, F-34095 Montpellier, France</footnote>

<footnote n="5"><sup>5</sup>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote>* E-mail: lisa.lehmann@irap.omp.eu</footnote>
# SPIRou reveals unusually strong magnetic fields of slowly rotating M dwarfs

L. T. Lehmann<sup>1,★</sup> J.-F. Donati<sup>1</sup>, P. Fouqué,<sup>1</sup> C. Moutou,<sup>1</sup> S. Bellotti,<sup>1,2</sup> X. Delfosse,<sup>3</sup> P. Petit,<sup>1</sup> A. Carmona,<sup>3</sup> J. Morin<sup>1,4</sup> and A. A. Vidotto<sup>1,5</sup>

the SLS consortium

Accepted 2023 November 7. Received 2023 November 5; in original form 2023 August 7

## ABSTRACT

In this paper, we study six slowly rotating mid-to-late M dwarfs (rotation period $P_{\text{rot}} \approx 40-190$ d) by analysing spectropolarimetric data collected with SpectroPolarimetre InfraRouge (SPIRou) at the Canada-France-Hawaii Telescope as part of the SPIRou Legacy Survey from 2019 to 2022. From ≈100–200 least-squares-deconvolved (LSD) profiles of circularly polarized 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 polarized 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, and GJ 1286) show a larger amount of temporal variations (mainly in field strength and axisymmetry) than the two partly convective ones (Gl 617B and Gl 408). Surprisingly, the six M dwarfs show large-scale field strengths in the range between 20 and 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.

**Key words:** techniques: polarimetric – stars: imaging – stars: low-mass – stars: magnetic field – stars: rotation.

## 1 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 non-axisymmetric 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 polarized and unpolarized Zeeman signatures on line 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 $Ro$, 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 $Ro$ until saturation occurs at $Ro \sim 0.1$ and below. Whereas large-scale fields of M dwarfs featuring $Ro < 1$ have been extensively studied, very little is known about those of very slow rotators with $Ro \sim 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 near-infrared 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

<footnote n="1"><sup>1</sup>IRAP, Université de Toulouse, CNRS / UMR 5277, UPS-OMP, 14 Avenue E. Belin, F-31400 Toulouse, France</footnote>

<footnote n="2"><sup>2</sup>Science Division, Directorate of Science, European Space Research and Technology Centre (ESA/ESTEC), Keplerlaan 1, NL-2201 AZ Noordwijk, the Netherlands</footnote>

<footnote n="3"><sup>3</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4"><sup>4</sup>Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS, F-34095 Montpellier, France</footnote>

<footnote n="5"><sup>5</sup>Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote>* E-mail: lisa.lehmann@irap.omp.eu</footnote>
<table><caption>**Table 1.** The stellar characteristics of our sample are from Cristofari et al. (2022) including spectral type, effective temperature $T_{\text{eff}}$, stellar mass, stellar radius $R_*$, metallicity [M/H], and log g. The rotation period $P_{\text{rot}}$ is copied from D23. For the Rossby number $Ro = P_{\text{rot}}/\tau$, we use $P_{\text{rot}}$ (eighth column) and determine the convective turnover time $\tau$ via Wright et al. (2018, equation 6). In the last column, we give the projected equatorial velocity $v_e \sin i = \frac{2\pi R_*}{P_{\text{rot}}} \sin i$, determined from $R_*$ (fifth column), $P_{\text{rot}}$ (eighth column), and an assumed inclination of $i = 60^{\circ}$ between the stellar rotation axis and the line of sight.</caption><thead><tr><th>Star</th><th>Spectral type</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>Mass<br>(M<sub>☉</sub>)</th><th>Radius<br>(R<sub>☉</sub>)</th><th>[M/H]</th><th>log g</th><th><i>P</i><sub>rot</sub><br>(d)</th><th><i>Ro</i></th><th><i>v</i><sub>e</sub> Sin i<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Gl 905</td><td>M5.0V</td><td>3069 ± 31</td><td>0.15 ± 0.02</td><td>0.165 ± 0.004</td><td>0.05 ± 0.11</td><td>4.78 ± 0.08</td><td>114.3 ± 2.8</td><td>0.88</td><td>0.06</td></tr><tr><td>GJ 1289</td><td>M4.5V</td><td>3238 ± 32</td><td>0.21 ± 0.02</td><td>0.233 ± 0.005</td><td>0.05 ± 0.10</td><td>5.00 ± 0.07</td><td>73.66 ± 0.92</td><td>0.67</td><td>0.14</td></tr><tr><td>GJ 1151</td><td>M4.5V</td><td>3178 ± 31</td><td>0.17 ± 0.02</td><td>0.193 ± 0.004</td><td>-0.04 ± 0.10</td><td>4.71 ± 0.06</td><td>175.6 ± 4.9</td><td>1.45</td><td>0.05</td></tr><tr><td>GJ 1286</td><td>M5.5V</td><td>2961 ± 33</td><td>0.12 ± 0.02</td><td>0.142 ± 0.004</td><td>-0.23 ± 0.10</td><td>4.55 ± 0.12</td><td>178 ± 15</td><td>1.25</td><td>0.03</td></tr><tr><td>Gl 617B</td><td>M3V</td><td>3525 ± 31</td><td>0.45 ± 0.02</td><td>0.460 ± 0.008</td><td>0.20 ± 0.10</td><td>4.84 ± 0.06</td><td>40.4 ± 3.0</td><td>0.77</td><td>0.50</td></tr><tr><td>Gl 408</td><td>M2.5V</td><td>3487 ± 31</td><td>0.38 ± 0.02</td><td>0.390 ± 0.007</td><td>-0.09 ± 0.10</td><td>4.79 ± 0.05</td><td>171.0 ± 8.4</td><td>2.68</td><td>0.10</td></tr></tbody></table>

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 seven 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 (Bellotti et al. 2023; Donati et al. 2023, hereinafter D23; Fouqué et al. 2023).

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 Sections 4–9. We conclude and discuss our results in Section 10.

## 2 SPIROU OBSERVATIONS

We analyse here a total of 986 circularly polarized spectra collected with SPIRou. The spectra span a wavelength range of 0.95–2.5 µm in the near-infrared with a resolving power of $R = 70000$. 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 ESPaDOnS at CFHT optimized for spectropolarimetry and specifically adapted for SPIRou (Donati et al. 2020).

We applied least-squares deconvolution (LSD; Donati et al. 1997) to all reduced unpolarized (Stokes I) and circularly polarized (Stokes V) spectra using an M3 mask constructed from outputs of the VALD-3 data base (Ryabchikova et al. 2015) assuming a temperature $T_{\text{eff}} = 3500$ K, a logarithmic surface gravity $log g = 5$, and a solar metallicity [M/H]. Although our six stars do not have the exact same atmospheric properties (see Table 1), we none the less 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 per cent and resulting in 575 lines for the mask. For further details see D23, section 2. The ephemeris used to calculate the phase and rotation cycle in this paper are summarized in Table 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 (QP) Gaussian Process regression

(GPR) to times series of their longitudinal fields $B_\ell$, 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 (optimized for RV precision; Cook et al. 2022) were 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.

## 3 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 V profiles.

### 3.1 PCA analysis of the LSD Stokes V profiles

Lehmann & Donati (2022) proposed a method to retrieve key information about the large-scale stellar magnetic field directly from time series of Stokes V LSD profiles without the need of an elaborate model of the field topology or several stellar parameters (e.g. the projected equatorial velocity $v_e \sin i$ or the inclination of the stellar rotation axis $i$). The method provides information about the axisymmetry, the poloidal/toroidal fraction of the axisymmetric component, the field complexity, and their evolution with time.

One first determines the mean profile of the whole Stokes V 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 V 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 large-scale 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 V time series removing the signal of axisymmetric field. The resulting mean-subtracted Stokes V 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
<table><caption>**Table 1.** The stellar characteristics of our sample are from Cristofari et al. (2022) including spectral type, effective temperature $T_{\text{eff}}$, stellar mass, stellar radius $R_*$, metallicity [M/H], and log g. The rotation period $P_{\text{rot}}$ is copied from D23. For the Rossby number $Ro = P_{\text{rot}}/\tau$, we use $P_{\text{rot}}$ (eighth column) and determine the convective turnover time $\tau$ via Wright et al. (2018, equation 6). In the last column, we give the projected equatorial velocity $v_e \sin i = \frac{2\pi R_*}{P_{\text{rot}}} \sin i$, determined from $R_*$ (fifth column), $P_{\text{rot}}$ (eighth column), and an assumed inclination of $i = 60^{\circ}$ between the stellar rotation axis and the line of sight.</caption><thead><tr><th>Star</th><th>Spectral type</th><th><i>T</i><sub>eff</sub><br>(K)</th><th>Mass<br>(M<sub>☉</sub>)</th><th>Radius<br>(R<sub>☉</sub>)</th><th>[M/H]</th><th>log g</th><th><i>P</i><sub>rot</sub><br>(d)</th><th><i>Ro</i></th><th><i>v</i><sub>e</sub> Sin i<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Gl 905</td><td>M5.0V</td><td>3069 ± 31</td><td>0.15 ± 0.02</td><td>0.165 ± 0.004</td><td>0.05 ± 0.11</td><td>4.78 ± 0.08</td><td>114.3 ± 2.8</td><td>0.88</td><td>0.06</td></tr><tr><td>GJ 1289</td><td>M4.5V</td><td>3238 ± 32</td><td>0.21 ± 0.02</td><td>0.233 ± 0.005</td><td>0.05 ± 0.10</td><td>5.00 ± 0.07</td><td>73.66 ± 0.92</td><td>0.67</td><td>0.14</td></tr><tr><td>GJ 1151</td><td>M4.5V</td><td>3178 ± 31</td><td>0.17 ± 0.02</td><td>0.193 ± 0.004</td><td>-0.04 ± 0.10</td><td>4.71 ± 0.06</td><td>175.6 ± 4.9</td><td>1.45</td><td>0.05</td></tr><tr><td>GJ 1286</td><td>M5.5V</td><td>2961 ± 33</td><td>0.12 ± 0.02</td><td>0.142 ± 0.004</td><td>-0.23 ± 0.10</td><td>4.55 ± 0.12</td><td>178 ± 15</td><td>1.25</td><td>0.03</td></tr><tr><td>Gl 617B</td><td>M3V</td><td>3525 ± 31</td><td>0.45 ± 0.02</td><td>0.460 ± 0.008</td><td>0.20 ± 0.10</td><td>4.84 ± 0.06</td><td>40.4 ± 3.0</td><td>0.77</td><td>0.50</td></tr><tr><td>Gl 408</td><td>M2.5V</td><td>3487 ± 31</td><td>0.38 ± 0.02</td><td>0.390 ± 0.007</td><td>-0.09 ± 0.10</td><td>4.79 ± 0.05</td><td>171.0 ± 8.4</td><td>2.68</td><td>0.10</td></tr></tbody></table>

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 seven 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 (Bellotti et al. 2023; Donati et al. 2023, hereinafter D23; Fouqué et al. 2023).

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 Sections 4–9. We conclude and discuss our results in Section 10.

## 2 SPIROU OBSERVATIONS

We analyse here a total of 986 circularly polarized spectra collected with SPIRou. The spectra span a wavelength range of 0.95–2.5 µm in the near-infrared with a resolving power of $R = 70000$. 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 ESPaDOnS at CFHT optimized for spectropolarimetry and specifically adapted for SPIRou (Donati et al. 2020).

We applied least-squares deconvolution (LSD; Donati et al. 1997) to all reduced unpolarized (Stokes I) and circularly polarized (Stokes V) spectra using an M3 mask constructed from outputs of the VALD-3 data base (Ryabchikova et al. 2015) assuming a temperature $T_{\text{eff}} = 3500$ K, a logarithmic surface gravity $\log g = 5$, and a solar metallicity [M/H]. Although our six stars do not have the exact same atmospheric properties (see Table 1), we none the less 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 per cent and resulting in 575 lines for the mask. For further details see D23, section 2. The ephemeris used to calculate the phase and rotation cycle in this paper are summarized in Table 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 (QP) Gaussian Process regression

(GPR) to times series of their longitudinal fields $B_\ell$, 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 (optimized for RV precision; Cook et al. 2022) were 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.

## 3 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 V profiles.

### 3.1 PCA analysis of the LSD Stokes V profiles

Lehmann & Donati (2022) proposed a method to retrieve key information about the large-scale stellar magnetic field directly from time series of Stokes V LSD profiles without the need of an elaborate model of the field topology or several stellar parameters (e.g. the projected equatorial velocity $v_e \sin i$ or the inclination of the stellar rotation axis $i$). The method provides information about the axisymmetry, the poloidal/toroidal fraction of the axisymmetric component, the field complexity, and their evolution with time.

One first determines the mean profile of the whole Stokes V 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 V 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 large-scale 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 V time series removing the signal of axisymmetric field. The resulting mean-subtracted Stokes V 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
<figure><img src="image_3.png" /><figcaption>**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 V 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 V profiles. **c.** The mean profile (left column), $c_1$ (the scaled and translated first PCA coefficient introduced in Section 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.</figcaption></figure>

coefficients. In the weighted PCA, the Stokes V 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 seven 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 V 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. 1 c left column). To study the evolution of the non-axisymmetric field, we compare the coefficients of the different seasons (see e.g. Fig. 1 c middle and right column). We caution that the coefficients are derived from the weighted PCA of the mean-subtracted Stokes V 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. 1 c left column). The usage of the weighted mean Stokes V 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 non-axisymmetric 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 $v_e \sin i$. As all our stars have $v_e \sin i \le 0.5$ km s⁻¹, 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.

## 3.2 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-subtracted Stokes V 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 $B_\ell$ (presented by D23) for the individual stars, we scale and translate the first coefficient, which we call $c_1$, using a linear model (scaling factor and offset), that minimizes the distance between the first coefficients and $B_\ell$ taking into account the measurement errors on $B_\ell$.

We can re-determine the stellar rotation period $P_{\text{rot}}$ of our six M dwarfs using a QP GPR fit to $c_1$ allowing us a direct comparison with the QP GPR results of $B_\ell$ presented by D23.

In contrast to D23, we use the PYTHON model presented by Martioli et al. (2022) based on (0:sc)george(0:sc) (Ambikasaran et al. 2015). Our adapted covariance function (or kernel) is given by

$$k(t_{ij}) = \alpha^2 \exp \left[ -\frac{t_{ij}^2}{2l^2} - \frac{1}{2\beta^2} \sin^2 \left( \frac{\pi t_{ij}}{P_{\text{rot}}} \right) \right], \quad (1)$$

where $t_{ij} = t_i - t_j$ is the time difference between the observations i and j, $\alpha$ is the amplitude of the Gaussian Process (GP), l is the decay time describing the typical time-scale on which the modulation pattern evolves, $\beta$ is the smoothing factor indicating the harmonic complexity of the QP modulation (lower values indicating higher harmonic complexity), and $P_{\text{rot}}$ is our new estimate of the stellar rotation period. The GP model parameters are fitted by maximizing the following likelihood function $\mathcal{L}$ using the PYTHON package SCIPY.OPTIMIZE:

$$\log \mathcal{L} = -\frac{1}{2} (N \log 2\pi + \log |\mathbf{K} + \mathbf{\Sigma} + \mathbf{S}| + y^T (\mathbf{K} + \mathbf{\Sigma} + \mathbf{S})^{-1} y), \quad (2)$$

where $\mathbf{K}$ is the QP kernel covariance matrix, $\mathbf{\Sigma}$ is the diagonal variance matrix of $c_1$, $\mathbf{S}$ is the diagonal matrix $\sigma^2 \mathbf{I}$ (with $\sigma$ an added amount of uncorrelated white noise (Angus et al. 2018) and $\mathbf{I}$ the identity matrix), $N$ is the number of observations, and $y$ corresponds to $c_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. Table 2 provides a summary of the results for all $c_1$ GPR fits in this study. For three stars, the decay time l was fixed as in D23 (see Table 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 $B_\ell$ values of D23 (see appendix B). This allows a direct comparison of the GP results for the same $B_\ell$ data set with our GP routine and the GP
<figure><img src="image_2.png" /><figcaption>**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 V 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 V profiles. **c.** The mean profile (left column), $c_1$ (the scaled and translated first PCA coefficient introduced in Section 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.</figcaption></figure>

coefficients. In the weighted PCA, the Stokes V 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 seven 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 V 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. 1 c left column). To study the evolution of the non-axisymmetric field, we compare the coefficients of the different seasons (see e.g. Fig. 1 c middle and right column). We caution that the coefficients are derived from the weighted PCA of the mean-subtracted Stokes V 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. 1 c left column). The usage of the weighted mean Stokes V 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 non-axisymmetric 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 $v_e \sin i$. As all our stars have $v_e \sin i \le 0.5$ km s$^{-1}$, we are likely to miss large-scale toroidal fields. We provide a typical $1\sigma$ error bar for the axisymmetric toroidal field for each star and each observing season of our sample.

## 3.2 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-subtracted Stokes V 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 $B_\ell$ (presented by D23) for the individual stars, we scale and translate the first coefficient, which we call $c_1$, using a linear model (scaling factor and offset), that minimizes the distance between the first coefficients and $B_\ell$ taking into account the measurement errors on $B_\ell$.

We can re-determine the stellar rotation period $P_{\text{rot}}$ of our six M dwarfs using a QP GPR fit to $c_1$ allowing us a direct comparison with the QP GPR results of $B_\ell$ presented by D23.

In contrast to D23, we use the PYTHON model presented by Martioli et al. (2022) based on (0:sc)george(0:sc) (Ambikasaran et al. 2015). Our adapted covariance function (or kernel) is given by

$$k(t_{ij}) = \alpha^2 \exp \left[ -\frac{t_{ij}^2}{2l^2} - \frac{1}{2\beta^2} \sin^2 \left( \frac{\pi t_{ij}}{P_{\text{rot}}} \right) \right], \quad (1)$$

where $t_{ij} = t_i - t_j$ is the time difference between the observations $i$ and $j$, $\alpha$ is the amplitude of the Gaussian Process (GP), $l$ is the decay time describing the typical time-scale on which the modulation pattern evolves, $\beta$ is the smoothing factor indicating the harmonic complexity of the QP modulation (lower values indicating higher harmonic complexity), and $P_{\text{rot}}$ is our new estimate of the stellar rotation period. The GP model parameters are fitted by maximizing the following likelihood function $\mathcal{L}$ using the PYTHON package SCIPY.OPTIMIZE:

$$\log \mathcal{L} = -\frac{1}{2} (N \log 2\pi + \log |\mathbf{K} + \mathbf{\Sigma} + \mathbf{S}| + y^T (\mathbf{K} + \mathbf{\Sigma} + \mathbf{S})^{-1} y), \quad (2)$$

where $\mathbf{K}$ is the QP kernel covariance matrix, $\mathbf{\Sigma}$ is the diagonal variance matrix of $c_1$, $\mathbf{S}$ is the diagonal matrix $\sigma^2 \mathbf{I}$ (with $\sigma$ an added amount of uncorrelated white noise (Angus et al. 2018) and $\mathbf{I}$ the identity matrix), $N$ is the number of observations, and $y$ corresponds to $c_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. Table 2 provides a summary of the results for all $c_1$ GPR fits in this study. For three stars, the decay time $l$ was fixed as in D23 (see Table 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 $B_\ell$ values of D23 (see appendix B). This allows a direct comparison of the GP results for the same $B_\ell$ data set with our GP routine and the GP
<table><caption>**Table 2.** Summary of the best-fitting parameters of the QP GPR fits applied to $c_1$ for the six M dwarfs in our sample, where rms is the root-mean-square of the residuals and $\chi_r^2$ is 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 $B_\ell$ data, and with those of D23, is given in Table B1.</caption><thead><tr><th>Star</th><th>Rotation period<br>P<sub>rot</sub> (d)</th><th>Decay time<br>l (d)</th><th>Smoothing factor<br>β</th><th>Amplitude<br>α (G)</th><th>White noise<br>σ (G)</th><th>rms<br>(G)</th><th>χ<sub>r</sub><sup>2</sup></th></tr></thead><tbody><tr><td>Gl 905</td><td>111.7<sup>+3.0</sup><sub>-3.2</sub></td><td>133<sup>+18</sup><sub>-22</sub></td><td>0.50<sup>+0.09</sup><sub>-0.07</sub></td><td>12.9<sup>+3.1</sup><sub>-2.1</sub></td><td>0.6<sup>+0.5</sup><sub>-0.4</sub></td><td>3.8</td><td>0.79</td></tr><tr><td>GJ 1289</td><td>75.62<sup>+0.85</sup><sub>-0.79</sub></td><td>129<sup>+26</sup><sub>-24</sub></td><td>0.40 ± 0.05</td><td>41.9<sup>+7.9</sup><sub>-6.0</sub></td><td>1.2<sup>+1.1</sup><sub>-0.8</sub></td><td>6.9</td><td>0.72</td></tr><tr><td>GJ 1151</td><td>175.8<sup>+3.2</sup><sub>-3.4</sub></td><td>300</td><td>0.40<sup>+0.10</sup><sub>-0.09</sub></td><td>12.7<sup>+3.2</sup><sub>-2.4</sub></td><td>3.3 ± 0.7</td><td>5.8</td><td>0.99</td></tr><tr><td>GJ 1286</td><td>186.8<sup>+9.5</sup><sub>-5.8</sub></td><td>300</td><td>0.23<sup>+0.05</sup><sub>-0.04</sub></td><td>18.0<sup>+4.3</sup><sub>-3.1</sub></td><td>2.7<sup>+1.5</sup><sub>-1.7</sub></td><td>7.5</td><td>1.02</td></tr><tr><td>Gl 617B</td><td>37.8<sup>+8.5</sup><sub>-2.6</sub></td><td>35<sup>+8</sup><sub>-4</sub></td><td>0.47<sup>+0.09</sup><sub>-0.05</sub></td><td>5.9<sup>+1.2</sup><sub>-0.8</sub></td><td>0.7<sup>+0.6</sup><sub>-0.5</sub></td><td>2.2</td><td>0.66</td></tr><tr><td>Gl 408</td><td>175<sup>+12</sup><sub>-14</sub></td><td>200</td><td>0.18<sup>+0.10</sup><sub>-0.05</sub></td><td>4.2<sup>+0.9</sup><sub>-0.8</sub></td><td>2.2<sup>+0.5</sup><sub>-0.6</sub></td><td>3.8</td><td>1.19</td></tr></tbody></table>

routine used by D23, as well as a comparison of the GP results for $c_1$ and $B_\ell$ obtained by the same GP routine. In general, we find that $c_1$ has lower rms, often shows lower $\sigma$ values, and provides smaller errors for $P_{rot}$ when the topology is not highly axisymmetric.

## 3.3 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 $\chi_r^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 $\alpha_{l,m}$ and $\beta_{l,m}$ coefficients of the poloidal component are modified as indicated in Lehmann & Donati (2022, equation B1). To compute synthetic Stokes profiles, the stellar surface is decomposed into a grid of 1000 cells. For each cell, the local Stokes V and I profiles are determined using Unno-Rachkovsky's analytical solution to the equations of the polarized 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 I LSD profiles beyond those attributable to radial velocity variations, so that we only use the Stokes V profiles for determining the magnetic field map via ZDI. Nevertheless, we make sure that the synthetic Stokes I profiles computed with ZDI agree well with the averaged observed Stokes I profile, especially in terms of width and depth. We assumed a fraction $f_V$ of each grid cell, which actually contributes to the Stokes V profile. This fraction $f_V$ is called filling factor of the large-scale field and is set equally to all cells, see also Morin et al. (2010) and Kochukhov (2021). For each star, we set $f_V = 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 $f_V = 0.1$ by finding lower $\chi_r^2$ values with $f_V = 0.1$ compared to $f_V = 1$ for each season of the different M dwarfs. The filling factor for the Stokes I profiles is set to $f_I = 1.0$ in consistency with the literature (Morin et al. 2010; Kochukhov 2021). The $v_e \sin i$ of our sample is $\le 0.5$ km s$^{-1}$ (see Table 1) and prevents us from reliably determining

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 $B_\ell$ and $c_1$ for most targets that cannot be obtained for pole-on viewed stars. Another reason is that higher values of i are intrinsically more likely than smaller ones. We restrict the spherical harmonics of the ZDI reconstructions to $\ell = 7$, as we see little magnetic energy stored in $\ell \sim 5-7$.

## 4 GL 905

The first star in our sample is the M5.5V dwarf Gl 905 (HH And, Ross 248) with a mass of $0.15 \pm 0.02$ M$_\odot$, (Cristofari et al. 2022). For our analysis, we use 219 Stokes IV LSD profiles observed with SPIRou between 2019 April and 2022 June and split the data in three seasons (2019 April–December, 2020 May–2021 January, and 2021 June–2022 January) for the per-season analysis. The 15 profiles collected in 2022 June at the beginning of a new season, only covering 9 per cent of a rotation cycle, were left out of the per-season analysis.

## 4.1 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 V 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 ($\chi_r^2 = 1.4$) but is likely due to an uneven phase coverage in the season 2020/21 (see Section 4.2). In Fig. 1(b), we show the first two eigenvectors of the mean-subtracted Stokes V 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 $c_1$, which appears very similar to the temporal evolution of $B_\ell$ 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 $c_1$ curve mimics that of $B_\ell$ (Lehmann & Donati 2022).

Our QP GPR model of $c_1$ finds a rotation period of $P_{rot} = 111.7^{+3.0}_{-3.2}$ d and decay time of $l = 133^{+18}_{-22}$ d very similar to the values derived by D23 ($P_{rot} = 114.3 \pm 2.8$ d and $l = 129^{+25}_{-21}$ d) and Fouqué et al. (2023) ($P_{rot} = 109.5^{+4.9}_{-5.4}$ d and $l = 149^{+26}_{-25}$ d) and also consistent with our GP fit of D23's $B_\ell$ values ($P_{rot} = 114.4^{+3.5}_{-2.4}$ d and $l = 130^{+25}_{-32}$ d), see Tables 2 and B1. For consistency, we use the rotation periods found by D23 to determine the rotation phase (see Table A1) and to model the ZDI maps for all six stars in our sample (see Table 1).
<table><caption>**Table 2.** Summary of the best-fitting parameters of the QP GPR fits applied to $c_1$ for the six M dwarfs in our sample, where rms is the root-mean-square of the residuals and $\chi_r^2$ is 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 $B_\ell$ data, and with those of D23, is given in Table B1.</caption><thead><tr><th>Star</th><th>Rotation period<br>P<sub>rot</sub> (d)</th><th>Decay time<br>l (d)</th><th>Smoothing factor<br>β</th><th>Amplitude<br>α (G)</th><th>White noise<br>σ (G)</th><th>rms<br>(G)</th><th>χ<sub>r</sub><sup>2</sup></th></tr></thead><tbody><tr><td>Gl 905</td><td>111.7<sup>+3.0</sup><sub>-3.2</sub></td><td>133<sup>+18</sup><sub>-22</sub></td><td>0.50<sup>+0.09</sup><sub>-0.07</sub></td><td>12.9<sup>+3.1</sup><sub>-2.1</sub></td><td>0.6<sup>+0.5</sup><sub>-0.4</sub></td><td>3.8</td><td>0.79</td></tr><tr><td>GJ 1289</td><td>75.62<sup>+0.85</sup><sub>-0.79</sub></td><td>129<sup>+26</sup><sub>-24</sub></td><td>0.40 ± 0.05</td><td>41.9<sup>+7.9</sup><sub>-6.0</sub></td><td>1.2<sup>+1.1</sup><sub>-0.8</sub></td><td>6.9</td><td>0.72</td></tr><tr><td>GJ 1151</td><td>175.8<sup>+3.2</sup><sub>-3.4</sub></td><td>300</td><td>0.40<sup>+0.10</sup><sub>-0.09</sub></td><td>12.7<sup>+3.2</sup><sub>-2.4</sub></td><td>3.3 ± 0.7</td><td>5.8</td><td>0.99</td></tr><tr><td>GJ 1286</td><td>186.8<sup>+9.5</sup><sub>-5.8</sub></td><td>300</td><td>0.23<sup>+0.05</sup><sub>-0.04</sub></td><td>18.0<sup>+4.3</sup><sub>-3.1</sub></td><td>2.7<sup>+1.5</sup><sub>-1.7</sub></td><td>7.5</td><td>1.02</td></tr><tr><td>Gl 617B</td><td>37.8<sup>+8.5</sup><sub>-2.6</sub></td><td>35<sup>+8</sup><sub>-4</sub></td><td>0.47<sup>+0.09</sup><sub>-0.05</sub></td><td>5.9<sup>+1.2</sup><sub>-0.8</sub></td><td>0.7<sup>+0.6</sup><sub>-0.5</sub></td><td>2.2</td><td>0.66</td></tr><tr><td>Gl 408</td><td>175<sup>+12</sup><sub>-14</sub></td><td>200</td><td>0.18<sup>+0.10</sup><sub>-0.05</sub></td><td>4.2<sup>+0.9</sup><sub>-0.8</sub></td><td>2.2<sup>+0.5</sup><sub>-0.6</sub></td><td>3.8</td><td>1.19</td></tr></tbody></table>

routine used by D23, as well as a comparison of the GP results for $c_1$ and $B_\ell$ obtained by the same GP routine. In general, we find that $c_1$ has lower rms, often shows lower $\sigma$ values, and provides smaller errors for $P_{rot}$ when the topology is not highly axisymmetric.

## 3.3 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 $\chi_r^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 $\alpha_{l,m}$ and $\beta_{l,m}$ coefficients of the poloidal component are modified as indicated in Lehmann & Donati (2022, equation B1). To compute synthetic Stokes profiles, the stellar surface is decomposed into a grid of 1000 cells. For each cell, the local Stokes V and I profiles are determined using Unno-Rachkovsky's analytical solution to the equations of the polarized 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 I LSD profiles beyond those attributable to radial velocity variations, so that we only use the Stokes V profiles for determining the magnetic field map via ZDI. Nevertheless, we make sure that the synthetic Stokes I profiles computed with ZDI agree well with the averaged observed Stokes I profile, especially in terms of width and depth. We assumed a fraction $f_V$ of each grid cell, which actually contributes to the Stokes V profile. This fraction $f_V$ is called filling factor of the large-scale field and is set equally to all cells, see also Morin et al. (2010) and Kochukhov (2021). For each star, we set $f_V = 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 $f_V = 0.1$ by finding lower $\chi_r^2$ values with $f_V = 0.1$ compared to $f_V = 1$ for each season of the different M dwarfs. The filling factor for the Stokes I profiles is set to $f_I = 1.0$ in consistency with the literature (Morin et al. 2010; Kochukhov 2021). The $v_e \sin i$ of our sample is $\le 0.5$ km s$^{-1}$ (see Table 1) and prevents us from reliably determining

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 $B_\ell$ and $c_1$ for most targets that cannot be obtained for pole-on viewed stars. Another reason is that higher values of $i$ are intrinsically more likely than smaller ones. We restrict the spherical harmonics of the ZDI reconstructions to $\ell = 7$, as we see little magnetic energy stored in $\ell \sim 5-7$.

## 4 GL 905

The first star in our sample is the M5.5V dwarf Gl 905 (HH And, Ross 248) with a mass of $0.15 \pm 0.02$ M$_\odot$, (Cristofari et al. 2022). For our analysis, we use 219 Stokes IV LSD profiles observed with SPIRou between 2019 April and 2022 June and split the data in three seasons (2019 April–December, 2020 May–2021 January, and 2021 June–2022 January) for the per-season analysis. The 15 profiles collected in 2022 June at the beginning of a new season, only covering 9 per cent of a rotation cycle, were left out of the per-season analysis.

## 4.1 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 V 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 ($\chi_r^2 = 1.4$) but is likely due to an uneven phase coverage in the season 2020/21 (see Section 4.2). In Fig. 1(b), we show the first two eigenvectors of the mean-subtracted Stokes V 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 $c_1$, which appears very similar to the temporal evolution of $B_\ell$ 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 $c_1$ curve mimics that of $B_\ell$ (Lehmann & Donati 2022).

Our QP GPR model of $c_1$ finds a rotation period of $P_{rot} = 111.7^{+3.0}_{-3.2}$ d and decay time of $l = 133^{+18}_{-22}$ d very similar to the values derived by D23 ($P_{rot} = 114.3 \pm 2.8$ d and $l = 129^{+25}_{-21}$ d) and Fouqué et al. (2023) ($P_{rot} = 109.5^{+4.9}_{-5.4}$ d and $l = 149^{+26}_{-25}$ d) and also consistent with our GP fit of D23's $B_\ell$ values ($P_{rot} = 114.4^{+3.5}_{-2.4}$ d and $l = 130^{+25}_{-32}$ d), see Tables 2 and B1. For consistency, we use the rotation periods found by D23 to determine the rotation phase (see Table A1) and to model the ZDI maps for all six stars in our sample (see Table 1).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Temporal variations of $c_1$ (top) and longitudinal field $B_\ell$ (bottom) for Gl 905. We show the QP GPR fit and its $1\sigma$ 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.</figcaption></figure>

In Fig. 1(c), we plot the mean profiles (left column) and the phase-folded 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 V 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 Sections 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 per cent of the observations concentrating between phase 0.3 and 0.75).

For the first season 2019, $c_1$ features a roughly sinusoidal behaviour indicating a mainly dipolar configuration. For 2020/21, $c_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 $c_1$ reaches its maximum is shifted, with $c_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

<figure><img src="image_5.png" /><figcaption>**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 summarizes 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).</figcaption></figure>

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.

## 4.2 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 summarized in Table 3. We were able to fit all three ZDI maps down to $\chi_r^2 \approx 1.0$ assuming $P_{\text{rot}} = 114.3$ d, $v_e \sin i = 0.06$ km s$^{-1}$, $i = 60^\circ$, and $f_V = 0.1$.

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 per cent to 4 per cent for the last season 2021/22. The surface
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Temporal variations of $c_1$ (top) and longitudinal field $B_\ell$ (bottom) for Gl 905. We show the QP GPR fit and its $1\sigma$ 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.</figcaption></figure>

In Fig. 1(c), we plot the mean profiles (left column) and the phase-folded 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 V 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 Sections 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 per cent of the observations concentrating between phase 0.3 and 0.75).

For the first season 2019, $c_1$ features a roughly sinusoidal behaviour indicating a mainly dipolar configuration. For 2020/21, $c_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 $c_1$ reaches its maximum is shifted, with $c_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

<figure><img src="image_4.png" /><figcaption>**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 summarizes 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).</figcaption></figure>

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.

## 4.2 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 summarized in Table 3. We were able to fit all three ZDI maps down to $\chi_r^2 \approx 1.0$ assuming $P_{\text{rot}} = 114.3$ d, $v_e \sin i = 0.06$ km s$^{-1}$, $i = 60^\circ$, and $f_V = 0.1$.

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 per cent to 4 per cent for the last season 2021/22. The surface
<table><caption>**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 $\langle B_V \rangle$[G], the surface average unsigned dipole magnetic field $\langle B_{dip} \rangle$[G], the typical $1\sigma$ error bar on the ZDI reconstructed surface averaged toroidal field $\sigma_{(B_{tor})}$, the fractional energy of the poloidal field, the fractional energy of the axisymmetric field (only $m = 0$ modes), the fractional energy of the dipole component, the tilt angle of the dipole ($\ell = 1$) with respect to the negative pole, the phase at which the dipole field faces the observer, the reduced $\chi^2$ values for the Stokes V profiles ($\chi^2_V$, corresponding to a $B = 0$ G fit), for the ZDI fit of the Stokes V profiles ($\chi^2_{V,ZDI}$) and for the Null profiles ($\chi^2_{r,N}$), and the number of observations per season (nb. obs).</caption><thead><tr><th>Season</th><th>2019</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 Apr</th><th>2020 May</th><th>2021 June</th></tr><tr><th>End</th><th>2019 Dec</th><th>2021 Jan</th><th>2022 Jan</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>128</td><td>89</td><td>64</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>124</td><td>80</td><td>64</td></tr><tr><td>σ<sub>(B<sub>tor</sub>)</sub> (G)</td><td>179</td><td>55</td><td>84</td></tr><tr><td>f<sub>pol</sub></td><td>1.0</td><td>0.99</td><td>1.0</td></tr><tr><td>f<sub>axi</sub></td><td>0.68</td><td>0.7</td><td>0.04</td></tr><tr><td>f<sub>dip</sub></td><td>0.93</td><td>0.79</td><td>0.94</td></tr><tr><td>Dipole tilt angle</td><td>33°</td><td>22°</td><td>83°</td></tr><tr><td>Pointing phase</td><td>0.29</td><td>0.12</td><td>0.96</td></tr><tr><td>χ<sup>2</sup><sub>r,V</sub></td><td>1.94</td><td>3.81</td><td>1.47</td></tr><tr><td>χ<sup>2</sup><sub>r,V,ZDI</sub></td><td>1.04</td><td>1.12</td><td>1.02</td></tr><tr><td>χ<sup>2</sup><sub>r,N</sub></td><td>1.22</td><td>1.10</td><td>1.04</td></tr><tr><td>Nb. obs</td><td>43</td><td>84</td><td>77</td></tr></tbody></table>

mean magnetic field decreases from 128 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 V 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 $\langle B_{tor} \rangle$ is lower than 10 G for each season. To derive a $1\sigma$ error bar on the simplest possible large-scale axisymmetric toroidal field (described with spherical harmonics coefficients $\ell = 1$ and $m = 0$) at the stellar surface, we proceed in the following way: (1) artificially add an axisymmetric toroidal field of strength $\langle B_{tor} \rangle$ to the reconstructed ZDI map, (2) simulate the corresponding Stokes V LSD profiles with the phase coverage and SNR of the actual observations, and (3) compute the new $\chi^2$ with respect to the observed LSD profiles and raise $\langle B_{tor} \rangle$ until $\chi^2$ is increased by 1 with respect to the optimal ZDI fit. We find $1\sigma$ error bars ranging from 180 G in 2019 down to 55 G in 2020/21.

## 5 GJ 1289

The fully convective M dwarf GJ 1289 ($M = 0.21 \pm 0.02$ M$_{\odot}$, Cristofari et al. 2022) is the next star in our sample. SPIRou observed GJ 1289 from 2019 September until 2022 June providing a time series of 204 LSD profiles split into three seasons (2019 June–December, 2020 May–2021 January, and 2021 June–2022 January) for the per-season analysis. As for Gl 905, the 14 profiles collected in 2022

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Same as Fig. 1 for GJ 1289.</figcaption></figure>

June at the beginning of a new season, only covering 16 per cent of a rotation cycle, were left out of the per-season analysis.

## 5.1 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 $c_1$ (see Fig. 5 top) finds $P_{rot} = 75.62^{+0.85}_{-0.79}$ d with a decay time of $l = 129^{+26}_{-24}$ d fitting all five parameters with a $\chi^2 = 0.72$ (see Table 2). The rotation period and decay time agree with the values found from the GPR fits of $B_\ell$ found by D23 ($P_{rot} = 73.66 \pm 0.92$ d and $l = 152^{+32}_{-27}$ d) and Fouqué et al. (2023) ($P_{rot} = 74.0^{+1.5}_{-1.3}$ d and $l = 142^{+33}_{-26}$ d) and are also consistent with our GP fit of the $B_\ell$ values of D23 ($P_{rot} = 73.67^{+1.01}_{-0.91}$ d and $l = 152^{+30}_{-28}$ d, see Fig. 5 bottom).

In Fig. 4(c), 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, season 2020/21 stands out, with the second eigenvector contributing significantly to the Stokes V signal before disappearing again for the last season 2021/22, for which $c_1$ shows a simpler trend.
<table><caption>**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 $\langle B_V \rangle$[G], the surface average unsigned dipole magnetic field $\langle B_{dip} \rangle$[G], the typical $1\sigma$ error bar on the ZDI reconstructed surface averaged toroidal field $\sigma_{(B_{tor})}$, the fractional energy of the poloidal field, the fractional energy of the axisymmetric field (only $m = 0$ modes), the fractional energy of the dipole component, the tilt angle of the dipole ($\ell = 1$) with respect to the negative pole, the phase at which the dipole field faces the observer, the reduced $\chi^2$ values for the Stokes V profiles ($\chi^2_V$, corresponding to a $B = 0$ G fit), for the ZDI fit of the Stokes V profiles ($\chi^2_{V,ZDI}$) and for the Null profiles ($\chi^2_{r,N}$), and the number of observations per season (nb. obs).</caption><thead><tr><th>Season</th><th>2019</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 Apr</th><th>2020 May</th><th>2021 June</th></tr><tr><th>End</th><th>2019 Dec</th><th>2021 Jan</th><th>2022 Jan</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>128</td><td>89</td><td>64</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>124</td><td>80</td><td>64</td></tr><tr><td>σ<sub>(B<sub>tor</sub>)</sub> (G)</td><td>179</td><td>55</td><td>84</td></tr><tr><td>f<sub>pol</sub></td><td>1.0</td><td>0.99</td><td>1.0</td></tr><tr><td>f<sub>axi</sub></td><td>0.68</td><td>0.7</td><td>0.04</td></tr><tr><td>f<sub>dip</sub></td><td>0.93</td><td>0.79</td><td>0.94</td></tr><tr><td>Dipole tilt angle</td><td>33°</td><td>22°</td><td>83°</td></tr><tr><td>Pointing phase</td><td>0.29</td><td>0.12</td><td>0.96</td></tr><tr><td>χ<sup>2</sup><sub>r,V</sub></td><td>1.94</td><td>3.81</td><td>1.47</td></tr><tr><td>χ<sup>2</sup><sub>r,V,ZDI</sub></td><td>1.04</td><td>1.12</td><td>1.02</td></tr><tr><td>χ<sup>2</sup><sub>r,N</sub></td><td>1.22</td><td>1.10</td><td>1.04</td></tr><tr><td>Nb. obs</td><td>43</td><td>84</td><td>77</td></tr></tbody></table>

mean magnetic field decreases from 128 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 V 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 $\langle B_{tor} \rangle$ is lower than 10 G for each season. To derive a $1\sigma$ error bar on the simplest possible large-scale axisymmetric toroidal field (described with spherical harmonics coefficients $\ell = 1$ and $m = 0$) at the stellar surface, we proceed in the following way: (1) artificially add an axisymmetric toroidal field of strength $\langle B_{tor} \rangle$ to the reconstructed ZDI map, (2) simulate the corresponding Stokes V LSD profiles with the phase coverage and SNR of the actual observations, and (3) compute the new $\chi^2$ with respect to the observed LSD profiles and raise $\langle B_{tor} \rangle$ until $\chi^2$ is increased by 1 with respect to the optimal ZDI fit. We find $1\sigma$ error bars ranging from 180 G in 2019 down to 55 G in 2020/21.

## 5 GJ 1289

The fully convective M dwarf GJ 1289 ($M = 0.21 \pm 0.02$ M$_{\odot}$, Cristofari et al. 2022) is the next star in our sample. SPIRou observed GJ 1289 from 2019 September until 2022 June providing a time series of 204 LSD profiles split into three seasons (2019 June–December, 2020 May–2021 January, and 2021 June–2022 January) for the per-season analysis. As for Gl 905, the 14 profiles collected in 2022

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Same as Fig. 1 for GJ 1289.</figcaption></figure>

June at the beginning of a new season, only covering 16 per cent of a rotation cycle, were left out of the per-season analysis.

## 5.1 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 $c_1$ (see Fig. 5 top) finds $P_{rot} = 75.62_{-0.79}^{+0.85}$ d with a decay time of $l = 129_{-24}^{+26}$ d fitting all five parameters with a $\chi^2 = 0.72$ (see Table 2). The rotation period and decay time agree with the values found from the GPR fits of $B_\ell$ found by D23 ($P_{rot} = 73.66 \pm 0.92$ d and $l = 152_{-27}^{+32}$ d) and Fouqué et al. (2023) ($P_{rot} = 74.0_{-1.3}^{+1.5}$ d and $l = 142_{-26}^{+33}$ d) and are also consistent with our GP fit of the $B_\ell$ values of D23 ($P_{rot} = 73.67_{-0.91}^{+1.01}$ d and $l = 152_{-28}^{+30}$ d, see Fig. 5 bottom).

In Fig. 4(c), 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, season 2020/21 stands out, with the second eigenvector contributing significantly to the Stokes V signal before disappearing again for the last season 2021/22, for which $c_1$ shows a simpler trend.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Same as Fig. 2 for GJ 1289.</figcaption></figure>

<table><caption>**Table 4.** Same as Table 3 for GJ 1289. The tilt angle now refers to the positive pole.</caption><thead><tr><th>Season</th><th>2019</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 June</th><th>2020 May</th><th>2021 June</th></tr><tr><th>End</th><th>2019 Dec</th><th>2021 Jan</th><th>2022 Jan</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>83</td><td>199</td><td>214</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>79</td><td>194</td><td>200</td></tr><tr><td>σ<sub>(B<sub>rot</sub>)</sub> (G)</td><td>42</td><td>57</td><td>100</td></tr><tr><td>f<sub>pol</sub></td><td>0.99</td><td>0.98</td><td>0.99</td></tr><tr><td>f<sub>axi</sub></td><td>0.93</td><td>0.86</td><td>0.57</td></tr><tr><td>f<sub>dip</sub></td><td>0.67</td><td>0.79</td><td>0.82</td></tr><tr><td>Dipole tilt angle</td><td>11°</td><td>8°</td><td>39°</td></tr><tr><td>Pointing phase</td><td>0.19</td><td>0.18</td><td>0.74</td></tr><tr><td>χ<sup>2</sup><sub>r,V</sub></td><td>1.46</td><td>3.60</td><td>8.34</td></tr><tr><td>χ<sup>2</sup><sub>r,V,ZDI</sub></td><td>0.98</td><td>1.04</td><td>0.98</td></tr><tr><td>χ<sup>2</sup><sub>r,N</sub></td><td>0.95</td><td>0.94</td><td>0.97</td></tr><tr><td>Nb. obs</td><td>35</td><td>80</td><td>75</td></tr></tbody></table>

## 5.2 ZDI reconstructions of GJ 1289

We were able to fit the Stokes V profiles for all seasons down to $\chi_r^2 \approx 1.0$ assuming $P_{\text{rot}} = 73.66$ d, $v_e \sin i = 0.14$ km s$^{-1}$, $i = 60^\circ$, and $f_V = 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 $\chi_{r,V}^2 = 1.46$ is several times lower than the $\chi_{r,V}^2$ of the following seasons (see Table 4). In the first season, ZDI reveals a weak marginally complex field topology. In 2020/21, the surface averaged field $\langle B_V \rangle$ becomes twice as large due to a growing

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Same as Fig. 3 for GJ 1289.</figcaption></figure>

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 Table 4). The toroidal field of the ZDI maps varies between 7 and 25 G, which is again lower than the typical $1\sigma$ error bar that we derive ranging from 40 to 100 G.

## 6 GJ 1151

Our next star, GJ 1151, is also a fully convective M dwarf ($M = 0.17 \pm 0.02$ M$_\odot$, Cristofari et al. 2022) and was observed between 2019 December and 2022 June with SPIRou providing us 158 LSD profiles (seasons: 2019 December–2020 July, 2020 December–2021 July, 2021 December–2022 June).

### 6.1 PCA analysis of GJ 1151

The mean profile is close to zero indicating a strongly non-axisymmetric topology (see Fig. 7a). Only the first eigenvector of the mean-subtracted Stokes V 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 $c_1$ down to a $\chi_r^2 = 0.99$ (see Fig. 8 top). We fix the decay time to 300 d similar to D23 and find a $P_{\text{rot}} = 175.8^{+3.2}_{-3.4}$ d similar to the results of D23 ($P_{\text{rot}} = 175.6 \pm 4.9$ d) and our own GP fit of $B_\ell$ ($P_{\text{rot}} = 176.1^{+3.6}_{-4.1}$ d), see also Fig. 8 bottom. Our rotation period is a bit higher than the one found by Fouqué et al. (2023) ($P_{\text{rot}} = 158 \pm 12$ d) but compatible at $1\sigma$.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Fig. 2 for GJ 1289.</figcaption></figure>

<table><caption>**Table 4.** Same as Table 3 for GJ 1289. The tilt angle now refers to the positive pole.</caption><thead><tr><th>Season</th><th>2019</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 June</th><th>2020 May</th><th>2021 June</th></tr><tr><th>End</th><th>2019 Dec</th><th>2021 Jan</th><th>2022 Jan</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>83</td><td>199</td><td>214</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>79</td><td>194</td><td>200</td></tr><tr><td>σ<sub>(B<sub>tor</sub>)</sub> (G)</td><td>42</td><td>57</td><td>100</td></tr><tr><td>f<sub>pol</sub></td><td>0.99</td><td>0.98</td><td>0.99</td></tr><tr><td>f<sub>axi</sub></td><td>0.93</td><td>0.86</td><td>0.57</td></tr><tr><td>f<sub>dip</sub></td><td>0.67</td><td>0.79</td><td>0.82</td></tr><tr><td>Dipole tilt angle</td><td>11°</td><td>8°</td><td>39°</td></tr><tr><td>Pointing phase</td><td>0.19</td><td>0.18</td><td>0.74</td></tr><tr><td>χ<sup>2</sup><sub>r,V</sub></td><td>1.46</td><td>3.60</td><td>8.34</td></tr><tr><td>χ<sup>2</sup><sub>r,V,ZDI</sub></td><td>0.98</td><td>1.04</td><td>0.98</td></tr><tr><td>χ<sup>2</sup><sub>r,N</sub></td><td>0.95</td><td>0.94</td><td>0.97</td></tr><tr><td>Nb. obs</td><td>35</td><td>80</td><td>75</td></tr></tbody></table>

## 5.2 ZDI reconstructions of GJ 1289

We were able to fit the Stokes V profiles for all seasons down to $\chi_r^2 \approx 1.0$ assuming $P_{\text{rot}} = 73.66$ d, $v_e \sin i = 0.14$ km s$^{-1}$, $i = 60^\circ$, and $f_V = 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 $\chi_{r,V}^2 = 1.46$ is several times lower than the $\chi_{r,V}^2$ of the following seasons (see Table 4). In the first season, ZDI reveals a weak marginally complex field topology. In 2020/21, the surface averaged field $\langle B_V \rangle$ becomes twice as large due to a growing

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Same as Fig. 3 for GJ 1289.</figcaption></figure>

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^\circ$ and dominates the field topology (see Table 4). The toroidal field of the ZDI maps varies between 7 and 25 G, which is again lower than the typical $1\sigma$ error bar that we derive ranging from 40 to 100 G.

## 6 GJ 1151

Our next star, GJ 1151, is also a fully convective M dwarf ($M = 0.17 \pm 0.02$ M$_\odot$, Cristofari et al. 2022) and was observed between 2019 December and 2022 June with SPIRou providing us 158 LSD profiles (seasons: 2019 December–2020 July, 2020 December–2021 July, 2021 December–2022 June).

### 6.1 PCA analysis of GJ 1151

The mean profile is close to zero indicating a strongly non-axisymmetric topology (see Fig. 7a). Only the first eigenvector of the mean-subtracted Stokes V 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 $c_1$ down to a $\chi_r^2 = 0.99$ (see Fig. 8 top). We fix the decay time to 300 d similar to D23 and find a $P_{\text{rot}} = 175.8^{+3.2}_{-3.4}$ d similar to the results of D23 ($P_{\text{rot}} = 175.6 \pm 4.9$ d) and our own GP fit of $B_\ell$ ($P_{\text{rot}} = 176.1^{+3.6}_{-4.1}$ d), see also Fig. 8 bottom. Our rotation period is a bit higher than the one found by Fouqué et al. (2023) ($P_{\text{rot}} = 158 \pm 12$ d) but compatible at $1\sigma$.
<figure><img src="image_10.png" /><figcaption>Figure 7. Same as Fig. 1 for GJ 1151.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Same as Fig. 2 for GJ 1151.</figcaption></figure>

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 $c_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 $c_1$ shows a higher amplitude ($\chi_r^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).

## 6.2 ZDI reconstructions of GJ 1151

We could fit the Stokes V profiles down to $\chi_r^2 \approx 1.0$ assuming $P_{\text{rot}} = 175.6$ d, $f_V = 0.1$, $v_e \sin i = 0.05$ km s$^{-1}$, and inclination of $i = 60^\circ$.

For the first two seasons, the Stokes V profiles are weak (see Fig. C6) and so is the reconstructed field (see Fig. 9), with a dominant negative polarity in the upper hemisphere that is consistent with the corresponding mean profiles. We see a small increase of $\langle B_V \rangle$ in the second season of 2020/21 explaining the higher amplitude of $c_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$^\circ$), that flipped polarity, and the surface averaged field is more than twice as high as before ($\langle B_V \rangle = 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 $\langle B_{\text{tor}} \rangle$ is lower than 8 G for all seasons of GJ 1151 and we find that the $1\sigma$ error bar on the toroidal field ranges between 370 and 450 G (see Table 5).

## 7 GJ 1286

GJ 1286 (LHS 546) is the lowest mass M dwarf in our sample ($M = 0.12 \pm 0.02$ M$_\odot$, Cristofari et al. 2022). We analyse here 104 observations, which we split into two seasons (2020 June–December and 2021 August–December) for the per-season analysis. As the first and last season (2019 September–December and 2022 June) do not contain enough observations, we once more left 21 LSD profiles out of the per-season analysis.

### 7.1 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 $c_1$ finds a $P_{\text{rot}} = 186.8_{-5.8}^{+9.5}$ d for a fixed decay time of 300 d reaching a $\chi_r^2 = 1.02$ similar to the GPR fit of $B_\ell$ (see Fig. 11). The rotation period agrees with the values found by D23 and Fouqué et al. (2023) ($P_{\text{rot}} = 178 \pm 15$ d and $P_{\text{rot}} = 203_{-21}^{+14}$ d, respectively) and our own GP result for $B_\ell$ ($P_{\text{rot}} = 181_{-13}^{+18}$ d).

The mean profile of season 2020 is nearly twice as high as for 2021 (see Fig. 10c). Also $c_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 $c_1$ curve gets less complex with phase for 2021.
<figure><img src="image_9.png" /><figcaption>Figure 7. Same as Fig. 1 for GJ 1151.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Same as Fig. 2 for GJ 1151.</figcaption></figure>

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 $c_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 $c_1$ shows a higher amplitude ($\chi_r^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).

## 6.2 ZDI reconstructions of GJ 1151

We could fit the Stokes V profiles down to $\chi_r^2 \approx 1.0$ assuming $P_{\text{rot}} = 175.6$ d, $f_V = 0.1$, $v_e \sin i = 0.05$ km s$^{-1}$, and inclination of $i = 60^\circ$.

For the first two seasons, the Stokes V profiles are weak (see Fig. C6) and so is the reconstructed field (see Fig. 9), with a dominant negative polarity in the upper hemisphere that is consistent with the corresponding mean profiles. We see a small increase of $\langle B_V \rangle$ in the second season of 2020/21 explaining the higher amplitude of $c_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$^\circ$), that flipped polarity, and the surface averaged field is more than twice as high as before ($\langle B_V \rangle = 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 $\langle B_{\text{tor}} \rangle$ is lower than 8 G for all seasons of GJ 1151 and we find that the $1\sigma$ error bar on the toroidal field ranges between 370 and 450 G (see Table 5).

## 7 GJ 1286

GJ 1286 (LHS 546) is the lowest mass M dwarf in our sample ($M = 0.12 \pm 0.02$ M$_\odot$, Cristofari et al. 2022). We analyse here 104 observations, which we split into two seasons (2020 June–December and 2021 August–December) for the per-season analysis. As the first and last season (2019 September–December and 2022 June) do not contain enough observations, we once more left 21 LSD profiles out of the per-season analysis.

### 7.1 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 $c_1$ finds a $P_{\text{rot}} = 186.8_{-5.8}^{+9.5}$ d for a fixed decay time of 300 d reaching a $\chi_r^2 = 1.02$ similar to the GPR fit of $B_\ell$ (see Fig. 11). The rotation period agrees with the values found by D23 and Fouqué et al. (2023) ($P_{\text{rot}} = 178 \pm 15$ d and $P_{\text{rot}} = 203_{-21}^{+14}$ d, respectively) and our own GP result for $B_\ell$ ($P_{\text{rot}} = 181_{-13}^{+18}$ d).

The mean profile of season 2020 is nearly twice as high as for 2021 (see Fig. 10c). Also $c_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 $c_1$ curve gets less complex with phase for 2021.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Same as Fig. 3 for GJ 1151. Note that the magnetic field switches polarity.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Same as Fig. 1 for GJ 1286.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Same as Fig. 2 for GJ 1286.</figcaption></figure>

<table><caption>**Table 5.** Same as Table 3 for GJ 1151.</caption><thead><tr><th>Season</th><th>2019/20</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 Dec</th><th>2020 Dec</th><th>2021 Dec</th></tr><tr><th>End</th><th>2020 July</th><th>2021 July</th><th>2022 June</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>26</td><td>35</td><td>63</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>23</td><td>32</td><td>62</td></tr><tr><td>σ⟨B<sub>tor</sub>⟩ (G)</td><td>448</td><td>378</td><td>372</td></tr><tr><td>f<sub>pol</sub></td><td>0.99</td><td>0.99</td><td>0.98</td></tr><tr><td>f<sub>axi</sub></td><td>0.84</td><td>0.64</td><td>0.38</td></tr><tr><td>f<sub>dip</sub></td><td>0.64</td><td>0.72</td><td>0.75</td></tr><tr><td>Dipole tilt angle</td><td>4°</td><td>32°</td><td>-55°</td></tr><tr><td>Pointing phase</td><td>0.44</td><td>0.05</td><td>0.33</td></tr><tr><td>χ<sub>r</sub><sup>2</sup>, V</td><td>1.16</td><td>1.25</td><td>2.57</td></tr><tr><td>χ<sub>r</sub><sup>2</sup>, V, ZDI</td><td>1.00</td><td>0.87</td><td>0.91</td></tr><tr><td>χ<sub>r</sub><sup>2</sup>, N</td><td>1.00</td><td>1.01</td><td>0.88</td></tr><tr><td>Nb. obs</td><td>38</td><td>53</td><td>67</td></tr></tbody></table>

## 7.2 ZDI reconstructions of GJ 1286

We fitted the LSD Stokes V profiles of the two seasons for GJ 1286 down to $\chi_r^2 \approx 1.0$, assuming $P_{rot} = 178$ d, $f_V = 0.1$, $v_e \sin i = 0.03$ km s$^{-1}$, and $i = 60^\circ$.

As concluded from the PCA analysis, the ZDI maps confirm that the topology becomes simpler and weaker: the fractional energy of the dipole component $f_{dip}$ increases from 0.73 to 0.79, while $\langle B_V \rangle$ decreases by almost half from 113 to 71 G (see Fig. 12 and Table 6). The reconstructed toroidal field $\langle B_{tor} \rangle$ of the seasons are 9 and 6 G,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Same as Fig. 3 for GJ 1151. Note that the magnetic field switches polarity.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Fig. 1 for GJ 1286.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Same as Fig. 2 for GJ 1286.</figcaption></figure>

<table><caption>**Table 5.** Same as Table 3 for GJ 1151.</caption><thead><tr><th>Season</th><th>2019/20</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 Dec</th><th>2020 Dec</th><th>2021 Dec</th></tr><tr><th>End</th><th>2020 July</th><th>2021 July</th><th>2022 June</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>26</td><td>35</td><td>63</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>23</td><td>32</td><td>62</td></tr><tr><td>σ⟨B<sub>tor</sub>⟩ (G)</td><td>448</td><td>378</td><td>372</td></tr><tr><td>f<sub>pol</sub></td><td>0.99</td><td>0.99</td><td>0.98</td></tr><tr><td>f<sub>axi</sub></td><td>0.84</td><td>0.64</td><td>0.38</td></tr><tr><td>f<sub>dip</sub></td><td>0.64</td><td>0.72</td><td>0.75</td></tr><tr><td>Dipole tilt angle</td><td>4°</td><td>32°</td><td>-55°</td></tr><tr><td>Pointing phase</td><td>0.44</td><td>0.05</td><td>0.33</td></tr><tr><td>χ<sub>r</sub><sup>2</sup>, V</td><td>1.16</td><td>1.25</td><td>2.57</td></tr><tr><td>χ<sub>r</sub><sup>2</sup>, V, ZDI</td><td>1.00</td><td>0.87</td><td>0.91</td></tr><tr><td>χ<sub>r</sub><sup>2</sup>, N</td><td>1.00</td><td>1.01</td><td>0.88</td></tr><tr><td>Nb. obs</td><td>38</td><td>53</td><td>67</td></tr></tbody></table>

## 7.2 ZDI reconstructions of GJ 1286

We fitted the LSD Stokes V profiles of the two seasons for GJ 1286 down to $\chi_r^2 \approx 1.0$, assuming $P_{rot} = 178$ d, $f_V = 0.1$, $v_e \sin i = 0.03$ km s$^{-1}$, and $i = 60^\circ$.

As concluded from the PCA analysis, the ZDI maps confirm that the topology becomes simpler and weaker: the fractional energy of the dipole component $f_{dip}$ increases from 0.73 to 0.79, while $\langle B_V \rangle$ decreases by almost half from 113 to 71 G (see Fig. 12 and Table 6). The reconstructed toroidal field $\langle B_{tor} \rangle$ of the seasons are 9 and 6 G,
<table><caption>**Table 6.** Same as Table 4 for GJ 1286.</caption><thead><tr><th>Season</th><th>2020</th><th>2021</th></tr><tr><th>Start</th><th>2020 June</th><th>2021 Aug</th></tr><tr><th>End</th><th>2020 Dec</th><th>2021 Dec</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>113</td><td>71</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>103</td><td>67</td></tr><tr><td>σ<sub>(B<sub>tor</sub>)</sub> (G)</td><td>325</td><td>300</td></tr><tr><td>f<sub>pol</sub></td><td>0.99</td><td>0.99</td></tr><tr><td>f<sub>axi</sub></td><td>0.67</td><td>0.79</td></tr><tr><td>f<sub>dip</sub></td><td>0.73</td><td>0.79</td></tr><tr><td>Dipole tilt angle</td><td>28°</td><td>25°</td></tr><tr><td>Pointing phase</td><td>0.03</td><td>0.95</td></tr><tr><td>χ<sup>2</sup><sub>R,V</sub></td><td>2.21</td><td>1.65</td></tr><tr><td>χ<sup>2</sup><sub>R,V,ZDI</sub></td><td>0.87</td><td>0.94</td></tr><tr><td>χ<sup>2</sup><sub>R,N</sub></td><td>1.05</td><td>0.82</td></tr><tr><td>Nb. obs</td><td>38</td><td>45</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Same as Fig. 3 for GJ 1286.</figcaption></figure>

respectively, while the typical $1\sigma$ error bars on the axisymmetric toroidal field are 325 and 300 G, respectively.

## 8 GL 617B

Gl 617B (EW Dra, HIP 79762, LHS 3176) is a partly convective M dwarf with $M = 0.45 \pm 0.02$ M☉ (Cristofari et al. 2022) and was observed between 2019 September and 2022 June with SPIRou. Our following analysis is based on 144 LSD Stokes spectra, which we split into three seasons (2020 February–October, 2021 January–July, and 2022 March–June) for the per-season analysis. As for the

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Same as Fig. 1 for Gl 617B.</figcaption></figure>

other stars, the first 15 spectra collected in 2019 were left out of the per-season analysis.

## 8.1 PCA analysis of Gl 617B

We find that the mean profile is much larger than the mean-subtracted Stokes V 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 $c_1$ down to a $\chi^2 = 0.66$ finding a rotation period of $37.8_{-2.6}^{+8.5}$ d in agreement with the results of D23 ($P_{\text{rot}} = 40.4 \pm 3.0$ d) and our GP fit of $B_\ell$ ($P_{\text{rot}} = 40.6_{-4.4}^{+2.1}$ d, see Fig. 14). However, the decay time for the GP fit of $c_1$ with $l = 35_{-4}^{+8}$ d is shorter than the results determined from the $B_\ell$ curves ($l = 69_{-23}^{+35}$ d for D23 and $l = 82_{-30}^{+45}$ d from our own fit of $B_\ell$). Fouqué et al. (2023) found no clear periodic variation using the APERO pipeline reduced spectra of Gl 617B.

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). None the less, $c_1$ traces a varying non-axisymmetric component. Season 2020 shows the highest range in amplitude of $c_1$, indicating the largest dipole tilt angle of all three seasons, although it will still be small ($<20^\circ$) due to the predominantly axisymmetric topology of Gl 617B.
<table><caption>**Table 6.** Same as Table 4 for GJ 1286.</caption><thead><tr><th>Season</th><th>2020</th><th>2021</th></tr><tr><th>Start</th><th>2020 June</th><th>2021 Aug</th></tr><tr><th>End</th><th>2020 Dec</th><th>2021 Dec</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>113</td><td>71</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>103</td><td>67</td></tr><tr><td>σ<sub>(B<sub>tor</sub>)</sub> (G)</td><td>325</td><td>300</td></tr><tr><td>f<sub>pol</sub></td><td>0.99</td><td>0.99</td></tr><tr><td>f<sub>axi</sub></td><td>0.67</td><td>0.79</td></tr><tr><td>f<sub>dip</sub></td><td>0.73</td><td>0.79</td></tr><tr><td>Dipole tilt angle</td><td>28°</td><td>25°</td></tr><tr><td>Pointing phase</td><td>0.03</td><td>0.95</td></tr><tr><td>χ<sup>2</sup><sub>R,V</sub></td><td>2.21</td><td>1.65</td></tr><tr><td>χ<sup>2</sup><sub>R,V,ZDI</sub></td><td>0.87</td><td>0.94</td></tr><tr><td>χ<sup>2</sup><sub>R,N</sub></td><td>1.05</td><td>0.82</td></tr><tr><td>Nb. obs</td><td>38</td><td>45</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Same as Fig. 3 for GJ 1286.</figcaption></figure>

respectively, while the typical $1\sigma$ error bars on the axisymmetric toroidal field are 325 and 300 G, respectively.

## 8 GL 617B

Gl 617B (EW Dra, HIP 79762, LHS 3176) is a partly convective M dwarf with $M = 0.45 \pm 0.02$ M☉ (Cristofari et al. 2022) and was observed between 2019 September and 2022 June with SPIRou. Our following analysis is based on 144 LSD Stokes spectra, which we split into three seasons (2020 February–October, 2021 January–July, and 2022 March–June) for the per-season analysis. As for the

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Same as Fig. 1 for Gl 617B.</figcaption></figure>

other stars, the first 15 spectra collected in 2019 were left out of the per-season analysis.

## 8.1 PCA analysis of Gl 617B

We find that the mean profile is much larger than the mean-subtracted Stokes V 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 $c_1$ down to a $\chi^2 = 0.66$ finding a rotation period of $37.8_{-2.6}^{+8.5}$ d in agreement with the results of D23 ($P_{\text{rot}} = 40.4 \pm 3.0$ d) and our GP fit of $B_\ell$ ($P_{\text{rot}} = 40.6_{-4.4}^{+2.1}$ d, see Fig. 14). However, the decay time for the GP fit of $c_1$ with $l = 35_{-4}^{+8}$ d is shorter than the results determined from the $B_\ell$ curves ($l = 69_{-23}^{+35}$ d for D23 and $l = 82_{-30}^{+45}$ d from our own fit of $B_\ell$). Fouqué et al. (2023) found no clear periodic variation using the APERO pipeline reduced spectra of Gl 617B.

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). None the less, $c_1$ traces a varying non-axisymmetric component. Season 2020 shows the highest range in amplitude of $c_1$, indicating the largest dipole tilt angle of all three seasons, although it will still be small ($<20^\circ$) due to the predominantly axisymmetric topology of Gl 617B.
<figure><img src="image_17.png" /><figcaption>Figure 14. Same as Fig. 2 for Gl 617B.</figcaption></figure>

<table><caption>Table 7. Same as Table 4 for Gl 617B.</caption><thead><tr><th>Season</th><th>2020</th><th>2021</th><th>2022</th></tr><tr><th>Start</th><th>2020 Feb</th><th>2021 Jan</th><th>2022 Mar</th></tr><tr><th>End</th><th>2020 Oct</th><th>2021 July</th><th>2022 June</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>53</td><td>75</td><td>36</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>52</td><td>73</td><td>35</td></tr><tr><td>σ<sub>(B<sub>tor</sub>)</sub> (G)</td><td>7</td><td>13</td><td>6</td></tr><tr><td>f<sub>pol</sub></td><td>0.98</td><td>0.99</td><td>1.0</td></tr><tr><td>f<sub>axi</sub></td><td>0.98</td><td>0.97</td><td>0.98</td></tr><tr><td>f<sub>dip</sub></td><td>0.67</td><td>0.74</td><td>0.71</td></tr><tr><td>Dipole tilt angle</td><td>7°</td><td>4°</td><td>3°</td></tr><tr><td>Pointing phase</td><td>0.74</td><td>0.47</td><td>0.87</td></tr><tr><td>χ<sup>2</sup><sub>R,V</sub></td><td>2.37</td><td>3.08</td><td>1.59</td></tr><tr><td>χ<sup>2</sup><sub>R,V,ZDI</sub></td><td>0.99</td><td>0.95</td><td>1.00</td></tr><tr><td>χ<sup>2</sup><sub>R,N</sub></td><td>1.00</td><td>1.06</td><td>0.89</td></tr><tr><td>Nb. obs</td><td>70</td><td>26</td><td>33</td></tr></tbody></table>

## 8.2 ZDI reconstructions of Gl 617B

We could fit Gl 617B down to $\chi_r^2 \approx 1.0$ assuming $P_{rot} = 40.4$ d, $f_V = 0.1$, $v_e \sin i = 0.50$ km s$^{-1}$, and $i = 60^\circ$.

The ZDI maps confirm a very axisymmetric, poloidal configuration (see Fig. 15). The axisymmetry is always equal to or greater than 97 per cent, so variations of the non-axisymmetric field are difficult to see, but appear largest in 2021 (see Table 7). The data set in season 2020 shows the largest tilt angle ($7^\circ$) as predicted by the PCA analysis. The reconstructed toroidal field $\langle B_{tor} \rangle$ reaches 9 G for 2020, 6 G for 2021, and 2 G for 2022, while the estimated $1\sigma$

<figure><img src="image_18.png" /><figcaption>Figure 15. Same as Fig. 3 for Gl 617B.</figcaption></figure>

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 $v_e \sin i = 0.5$ km s$^{-1}$ of Gl 617B.

We see $\langle B_V \rangle$ changing by approximately $\pm 25$ G for the three seasons, otherwise the main properties of the maps are similar (see Table 7).

For highly axisymmetric topologies, it is difficult to infer the inclination $i$. It may be that $i$ is actually lower for Gl 617B. We provide the ZDI maps for an inclination $i = 30^\circ$ and $v_e \sin i = 0.29$ km s$^{-1}$, while otherwise using the same parameters (see Fig. C2). The $\chi_r^2$ values reached for the ZDI fits are slightly higher for $i = 30^\circ$ than for $i = 60^\circ$.

## 9 GL 408

Gl 408 (Ross 104, HIP 53767, LHS 6193) is another partly convective star, with $M = 0.38 \pm 0.02$ M$_\odot$, (Cristofari et al. 2022). SPIRou observed Gl 408 between 2019 April and 2022 June. We use 157 Stokes V profiles for the following analysis split into three seasons (2019 October–2020 June, 2020 October–2021 July, and 2021 November–2022 June) 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.

## 9.1 PCA analysis of Gl 408

Gl 408 has the strongest mean profile compared to the mean-subtracted profile in our sample (see Figs 16 a and C9). The first eigenvector of Gl 408 (the only one showing a signal) is already
<figure><img src="image_16.png" /><figcaption>Figure 14. Same as Fig. 2 for Gl 617B.</figcaption></figure>

<table><caption>Table 7. Same as Table 4 for Gl 617B.</caption><thead><tr><th>Season</th><th>2020</th><th>2021</th><th>2022</th></tr><tr><th>Start</th><th>2020 Feb</th><th>2021 Jan</th><th>2022 Mar</th></tr><tr><th>End</th><th>2020 Oct</th><th>2021 July</th><th>2022 June</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩ (G)</td><td>53</td><td>75</td><td>36</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>52</td><td>73</td><td>35</td></tr><tr><td>σ<sub>(B<sub>tor</sub>)</sub> (G)</td><td>7</td><td>13</td><td>6</td></tr><tr><td>f<sub>pol</sub></td><td>0.98</td><td>0.99</td><td>1.0</td></tr><tr><td>f<sub>axi</sub></td><td>0.98</td><td>0.97</td><td>0.98</td></tr><tr><td>f<sub>dip</sub></td><td>0.67</td><td>0.74</td><td>0.71</td></tr><tr><td>Dipole tilt angle</td><td>7°</td><td>4°</td><td>3°</td></tr><tr><td>Pointing phase</td><td>0.74</td><td>0.47</td><td>0.87</td></tr><tr><td>χ<sup>2</sup><sub>R,V</sub></td><td>2.37</td><td>3.08</td><td>1.59</td></tr><tr><td>χ<sup>2</sup><sub>R,V,ZDI</sub></td><td>0.99</td><td>0.95</td><td>1.00</td></tr><tr><td>χ<sup>2</sup><sub>R,N</sub></td><td>1.00</td><td>1.06</td><td>0.89</td></tr><tr><td>Nb. obs</td><td>70</td><td>26</td><td>33</td></tr></tbody></table>

## 8.2 ZDI reconstructions of Gl 617B

We could fit Gl 617B down to $\chi_r^2 \approx 1.0$ assuming $P_{rot} = 40.4$ d, $f_V = 0.1$, $v_e \sin i = 0.50$ km s$^{-1}$, and $i = 60^\circ$.

The ZDI maps confirm a very axisymmetric, poloidal configuration (see Fig. 15). The axisymmetry is always equal to or greater than 97 per cent, so variations of the non-axisymmetric field are difficult to see, but appear largest in 2021 (see Table 7). The data set in season 2020 shows the largest tilt angle ($7^\circ$) as predicted by the PCA analysis. The reconstructed toroidal field $\langle B_{tor} \rangle$ reaches 9 G for 2020, 6 G for 2021, and 2 G for 2022, while the estimated $1\sigma$

<figure><img src="image_17.png" /><figcaption>Figure 15. Same as Fig. 3 for Gl 617B.</figcaption></figure>

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 $v_e \sin i = 0.5$ km s$^{-1}$ of Gl 617B.

We see $\langle B_V \rangle$ changing by approximately $\pm 25$ G for the three seasons, otherwise the main properties of the maps are similar (see Table 7).

For highly axisymmetric topologies, it is difficult to infer the inclination $i$. It may be that $i$ is actually lower for Gl 617B. We provide the ZDI maps for an inclination $i = 30^\circ$ and $v_e \sin i = 0.29$ km s$^{-1}$, while otherwise using the same parameters (see Fig. C2). The $\chi_r^2$ values reached for the ZDI fits are slightly higher for $i = 30^\circ$ than for $i = 60^\circ$.

## 9 GL 408

Gl 408 (Ross 104, HIP 53767, LHS 6193) is another partly convective star, with $M = 0.38 \pm 0.02$ M<sub>☉</sub>, (Cristofari et al. 2022). SPIRou observed Gl 408 between 2019 April and 2022 June. We use 157 Stokes V profiles for the following analysis split into three seasons (2019 October–2020 June, 2020 October–2021 July, and 2021 November–2022 June) 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.

## 9.1 PCA analysis of Gl 408

Gl 408 has the strongest mean profile compared to the mean-subtracted profile in our sample (see Figs 16 a and C9). The first eigenvector of Gl 408 (the only one showing a signal) is already
<figure><img src="image_19.png" /><figcaption>Figure 16. Same as Fig. 1 for Gl 408.</figcaption></figure>

noisy, a strong indication of a very axisymmetric field topology (see Fig. 16b).

Fig. 17 (top) presents the QP GPR fit of $c_1$, which gives a $P_{\text{rot}} = 175_{-14}^{+12}$ d similar to D23 determining $P_{\text{rot}} = 171.0 \pm 8.4$ d. Fitting $B_V$ with our GP routines, we derive a $P_{\text{rot}} = 170.7_{-9.8}^{+7.1}$ 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 $\approx 200 \pm 70$ d but higher $\chi^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. 16(c), we see that $c_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.

## 9.2 ZDI reconstructions of Gl 408

All three seasons could be fitted down to $\chi^2 \approx 1.0$ assuming $P_{\text{rot}} = 171.0$ d, $f_V = 0.1$, $v_e \sin i = 0.10 \text{ km s}^{-1}$, and $i = 60^\circ$. The topology changes little over the three seasons and is characterized by a strong, axisymmetric, poloidal dipole of negative polarity (see Fig. 18). It is the most stable topology in our sample and only $\langle B_V \rangle$ varies marginally between 106 and 130 G (see Table 8). We find a $1\sigma$ error bar on the axisymmetric toroidal field of 55–81 G for Gl 408 whereas the reconstructed $\langle B_{\text{tor}} \rangle$ ranges between 4 and 13 G.

Similar to Gl 617B, we also determine the ZDI maps for an inclination of $i = 30^\circ$ and $v_e \sin i = 0.06 \text{ km s}^{-1}$ (see Fig. C3). The $\chi^2$ values of the ZDI fits are again slightly higher for the lower inclination $i = 30^\circ$ than for $i = 60^\circ$.

<figure><img src="image_20.png" /><figcaption>Figure 17. Same as Fig. 2 for Gl 408.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 18. Same as Fig. 3 for Gl 408.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 16. Same as Fig. 1 for Gl 408.</figcaption></figure>

noisy, a strong indication of a very axisymmetric field topology (see Fig. 16b).

Fig. 17 (top) presents the QP GPR fit of $c_1$, which gives a $P_{\text{rot}} = 175_{-14}^{+12}$ d similar to D23 determining $P_{\text{rot}} = 171.0 \pm 8.4$ d. Fitting $B_v$ with our GP routines, we derive a $P_{\text{rot}} = 170.7_{-9.8}^{+7.1}$ 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 $\approx 200 \pm 70$ d but higher $\chi^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. 16(c), we see that $c_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.

## 9.2 ZDI reconstructions of Gl 408

All three seasons could be fitted down to $\chi^2 \approx 1.0$ assuming $P_{\text{rot}} = 171.0$ d, $f_v = 0.1$, $v_e \sin i = 0.10 \text{ km s}^{-1}$, and $i = 60^\circ$. The topology changes little over the three seasons and is characterized by a strong, axisymmetric, poloidal dipole of negative polarity (see Fig. 18). It is the most stable topology in our sample and only $\langle B_v \rangle$ varies marginally between 106 and 130 G (see Table 8). We find a $1\sigma$ error bar on the axisymmetric toroidal field of 55–81 G for Gl 408 whereas the reconstructed $\langle B_{\text{tor}} \rangle$ ranges between 4 and 13 G.

Similar to Gl 617B, we also determine the ZDI maps for an inclination of $i = 30^\circ$ and $v_e \sin i = 0.06 \text{ km s}^{-1}$ (see Fig. C3). The $\chi^2$ values of the ZDI fits are again slightly higher for the lower inclination $i = 30^\circ$ than for $i = 60^\circ$.

<figure><img src="image_19.png" /><figcaption>Figure 17. Same as Fig. 2 for Gl 408.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 18. Same as Fig. 3 for Gl 408.</figcaption></figure>
<table><caption>**Table 8.** Same as Table 3 for Gl 408.</caption><thead><tr><th>Season</th><th>2019/20</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 Oct</th><th>2020 Oct</th><th>2021 Nov</th></tr><tr><th>End</th><th>2020 June</th><th>2021 July</th><th>2022 June</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩</td><td>130</td><td>106</td><td>120</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>129</td><td>104</td><td>117</td></tr><tr><td>σ<sub>⟨B<sub>tor</sub>⟩</sub> (G)</td><td>81</td><td>55</td><td>57</td></tr><tr><td>f<sub>pol</sub></td><td>1.0</td><td>1.0</td><td>0.99</td></tr><tr><td>f<sub>axi</sub></td><td>0.98</td><td>0.98</td><td>0.98</td></tr><tr><td>f<sub>dip</sub></td><td>0.82</td><td>0.77</td><td>0.78</td></tr><tr><td>Dipole tilt angle</td><td>6°</td><td>4°</td><td>4°</td></tr><tr><td>Pointing phase</td><td>0.53</td><td>0.44</td><td>0.47</td></tr><tr><td>χ<sup>2</sup><sub>R,V</sub></td><td>5.63</td><td>4.13</td><td>5.66</td></tr><tr><td>χ<sup>2</sup><sub>R,V,ZDI</sub></td><td>1.01</td><td>1.02</td><td>1.01</td></tr><tr><td>χ<sup>2</sup><sub>R,N</sub></td><td>0.89</td><td>0.99</td><td>1.02</td></tr><tr><td>Nb. obs</td><td>31</td><td>62</td><td>47</td></tr></tbody></table>

## 10 SUMMARY, DISCUSSION, AND 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 $c_1$, i.e. the scaled and translated first coefficient of the PCA analysis (see Section 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 $f_{axi} \ge 0.97$ between 2019 and 2022. The higher the axisymmetry of the large-scale field, the smaller the variations of $B_\ell$ or $c_1$ with time, and the harder it is to determine a rotation period. For the highly axisymmetric topologies, we find that the $\chi^2$ of the GPR fits increases, reflecting that in such cases, the $B_\ell$ curves are more sensitive to intrinsic variability and less to rotational modulation, and thereby reducing the ability at measuring rotation periods (see e.g. Figs A5 or A6 and Table B1).

Using the PCA analysis, we derive information about axisymmetry and complexity directly from the LSD Stokes V 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 $v_e \sin i$, inclinations, etc.

We find evidence for a polarity reversal of the large-scale field (via sign changes of $B_\ell$, $c_1$, or in the mean profiles) taking place on GJ 1151 and possibly also on Gl 905, for which the axisym- metric 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 V 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 V profiles. The lower the $v_e \sin i$, the higher the $1\sigma$ error bar on the toroidal field. The typical $1\sigma$ error bars on the toroidal field ranges from 6 to 450 G depending on SNR and $v_e \sin i$.

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

<figure><img src="image_22.png" /><figcaption>**Figure 19.** The averaged unsigned magnetic field strength $\langle B_V \rangle$ versus Rossby number $Ro$ 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 $Ro$ 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. 4 a of Vidotto et al. (2014).</figcaption></figure>

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 $\langle B_V \rangle \approx 20-200$ G. They show similar $\langle B_V \rangle$ to faster rotating M dwarfs in the saturated regime. We add our sample to fig. 4 a of Vidotto et al. (2014), which originally shows $\langle B_V \rangle$ versus Rossby number $Ro$ for 73 stars, including stars in the mass range ~0.1–2 M$_\odot$ (see Fig. 19). The grey open circles depict the mid- and late-type rapidly rotating M dwarfs of Vidotto et al. (2014), while the red circles show our slow-rotating M dwarfs. The solar-like G–K dwarfs (grey diamonds and pentagons) follow a decreasing trend with increasing $Ro$, while our M dwarfs show stronger $\langle B_V \rangle$ values than expected given their $Ro$. This is in agreement with the results of Medina et al. (2022), showing that M dwarfs can remain extremely active (flaring) even when their rotation period increases beyond 100 d. 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-$Ro$ stars of the Vidotto et al. (2014) sample. This will be the subject of forthcoming studies. Besides, our results suggest 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 summarizes 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 time-scales 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.
<table><caption>**Table 8.** Same as Table 3 for Gl 408.</caption><thead><tr><th>Season</th><th>2019/20</th><th>2020/21</th><th>2021/22</th></tr><tr><th>Start</th><th>2019 Oct</th><th>2020 Oct</th><th>2021 Nov</th></tr><tr><th>End</th><th>2020 June</th><th>2021 July</th><th>2022 June</th></tr></thead><tbody><tr><td>⟨B<sub>V</sub>⟩</td><td>130</td><td>106</td><td>120</td></tr><tr><td>⟨B<sub>dip</sub>⟩ (G)</td><td>129</td><td>104</td><td>117</td></tr><tr><td>σ<sub>⟨B<sub>tor</sub>⟩</sub> (G)</td><td>81</td><td>55</td><td>57</td></tr><tr><td>f<sub>pol</sub></td><td>1.0</td><td>1.0</td><td>0.99</td></tr><tr><td>f<sub>axi</sub></td><td>0.98</td><td>0.98</td><td>0.98</td></tr><tr><td>f<sub>dip</sub></td><td>0.82</td><td>0.77</td><td>0.78</td></tr><tr><td>Dipole tilt angle</td><td>6°</td><td>4°</td><td>4°</td></tr><tr><td>Pointing phase</td><td>0.53</td><td>0.44</td><td>0.47</td></tr><tr><td>χ<sub>τ,V</sub><sup>2</sup></td><td>5.63</td><td>4.13</td><td>5.66</td></tr><tr><td>χ<sub>τ,V,ZDI</sub><sup>2</sup></td><td>1.01</td><td>1.02</td><td>1.01</td></tr><tr><td>χ<sub>τ,N</sub><sup>2</sup></td><td>0.89</td><td>0.99</td><td>1.02</td></tr><tr><td>Nb. obs</td><td>31</td><td>62</td><td>47</td></tr></tbody></table>

## 10 SUMMARY, DISCUSSION, AND 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 $c_1$, i.e. the scaled and translated first coefficient of the PCA analysis (see Section 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 $f_{axi} \ge 0.97$ between 2019 and 2022. The higher the axisymmetry of the large-scale field, the smaller the variations of $B_\ell$ or $c_1$ with time, and the harder it is to determine a rotation period. For the highly axisymmetric topologies, we find that the $\chi^2$ of the GPR fits increases, reflecting that in such cases, the $B_\ell$ curves are more sensitive to intrinsic variability and less to rotational modulation, and thereby reducing the ability at measuring rotation periods (see e.g. Figs A5 or A6 and Table B1).

Using the PCA analysis, we derive information about axisymmetry and complexity directly from the LSD Stokes V 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 $v_e \sin i$, inclinations, etc.

We find evidence for a polarity reversal of the large-scale field (via sign changes of $B_\ell$, $c_1$, or in the mean profiles) taking place on GJ 1151 and possibly also on Gl 905, for which the axisym- metric 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 V 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 V profiles. The lower the $v_e \sin i$, the higher the $1\sigma$ error bar on the toroidal field. The typical $1\sigma$ error bars on the toroidal field ranges from 6 to 450 G depending on SNR and $v_e \sin i$.

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

<figure><img src="image_21.png" /><figcaption>**Figure 19.** The averaged unsigned magnetic field strength $\langle B_V \rangle$ versus Rossby number $Ro$ 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 $Ro$ 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. 4 a of Vidotto et al. (2014).</figcaption></figure>

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 $\langle B_V \rangle \approx 20-200$ G. They show similar $\langle B_V \rangle$ to faster rotating M dwarfs in the saturated regime. We add our sample to fig. 4 a of Vidotto et al. (2014), which originally shows $\langle B_V \rangle$ versus Rossby number $Ro$ for 73 stars, including stars in the mass range ~0.1–2 M$_\odot$ (see Fig. 19). The grey open circles depict the mid- and late-type rapidly rotating M dwarfs of Vidotto et al. (2014), while the red circles show our slow-rotating M dwarfs. The solar- like G–K dwarfs (grey diamonds and pentagons) follow a decreasing trend with increasing $Ro$, while our M dwarfs show stronger $\langle B_V \rangle$ values than expected given their $Ro$. This is in agreement with the results of Medina et al. (2022), showing that M dwarfs can remain extremely active (flaring) even when their rotation period increases beyond 100 d. Besides, it implies a harsher interplan- etary 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-$Ro$ stars of the Vidotto et al. (2014) sample. This will be the subject of forthcoming studies. Besides, our results suggest 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 summarizes 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 time-scales 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.
<figure><img src="image_23.png" /><figcaption>**Figure 20.** A summary of the magnetic properties of our M dwarf sample as a function of stellar mass and rotation period. The symbol size indicates the magnetic energy, the colour the poloidal fractional energy $f_{pol}$, and the shape the axisymmetric fractional energy $f_{axi}$ (see the legend on the right). For each star, we display two to three symbols representing the individual seasons (blue, grey, and orange border for the first, second, and third season) of each M dwarf. The thin grey line indicates $Ro = 1$ determined using the empirical relation of Wright et al. (2018). The thick grey line marks the fully convection limit at $M \approx 0.35 M_{\odot}$.</figcaption></figure>

In conclusion, we have analysed six slowly rotating M dwarfs observed by the SLS over 3.5 yr. 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 *Ro* 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 (GI 905) of the four 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 long-term 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.

## ACKNOWLEDGEMENTS

We acknowledge funding from the European Research Council (ERC) under the H2020 research & innovation programme (grant agreement #740651 NewWorlds). AC acknowledges funding from the French ANR under contract number ANR18CE310019 (SPlaSH). This work is supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the 'Origin of Life' project of the Grenoble-Alpes University. AAV acknowledges funding from the European Research Council (ERC) under the European

Union's Horizon 2020 research and innovation programme (grant agreement no. 817540, ASTROFLOW). Our study is based on data obtained at the CFHT, which is operated by the CNRC (Canada), INSU/CNRS (France), and the University of Hawaii. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are very fortunate to have the opportunity to conduct observations from this mountain and gratefully acknowledge the CFHT QSO observers. This work has used ASTROPY (Astropy Collaboration 2013, 2018, 2022), NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), SCIPY (Virtanen et al. 2020), and benefited from the SIMBAD CDS data base and the ADS system.

## DATA AVAILABILITY

All data underlying this paper are part of the SLS, and will be publicly available from the Canadian Astronomy Data Center by 2024 February.

## REFERENCES

Ambikasaran S., Foreman-Mackey D., Greengard L., Hogg D. W., O'Neil M., 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 2

Angus R., Morton T., Aigrain S., Foreman-Mackey D., Rajpaul V., 2018, MNRAS, 474, 2094

Astropy Collaboration, 2013, A&A, 558, A33

Astropy Collaboration, 2018, AJ, 156, 123

Astropy Collaboration, 2022, ApJ, 935, 167

Bellotti S. et al., 2023, A&A, 676, A56

Cook N. J. et al., 2022, Astrophysics Source Code Library, record ascl:2211.019
<figure><img src="image_22.png" /><figcaption>**Figure 20.** A summary of the magnetic properties of our M dwarf sample as a function of stellar mass and rotation period. The symbol size indicates the magnetic energy, the colour the poloidal fractional energy $f_{pol}$, and the shape the axisymmetric fractional energy $f_{axi}$ (see the legend on the right). For each star, we display two to three symbols representing the individual seasons (blue, grey, and orange border for the first, second, and third season) of each M dwarf. The thin grey line indicates $Ro = 1$ determined using the empirical relation of Wright et al. (2018). The thick grey line marks the fully convection limit at $M \approx 0.35 M_{\odot}$.</figcaption></figure>

In conclusion, we have analysed six slowly rotating M dwarfs observed by the SLS over 3.5 yr. 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 *Ro* 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 (GI 905) of the four 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 long-term 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.

## ACKNOWLEDGEMENTS

We acknowledge funding from the European Research Council (ERC) under the H2020 research & innovation programme (grant agreement #740651 NewWorlds). AC acknowledges funding from the French ANR under contract number ANR18CE310019 (SPlaSH). This work is supported by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the 'Origin of Life' project of the Grenoble-Alpes University. AAV acknowledges funding from the European Research Council (ERC) under the European

Union's Horizon 2020 research and innovation programme (grant agreement no. 817540, ASTROFLOW). Our study is based on data obtained at the CFHT, which is operated by the CNRC (Canada), INSU/CNRS (France), and the University of Hawaii. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are very fortunate to have the opportunity to conduct observations from this mountain and gratefully acknowledge the CFHT QSO observers. This work has used ASTROPY (Astropy Collaboration 2013, 2018, 2022), NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), SCIPY (Virtanen et al. 2020), and benefited from the SIMBAD CDS data base and the ADS system.

## DATA AVAILABILITY

All data underlying this paper are part of the SLS, and will be publicly available from the Canadian Astronomy Data Center by 2024 February.

## REFERENCES

* Ambikasaran S., Foreman-Mackey D., Greengard L., Hogg D. W., O'Neil M., 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 2

* Angus R., Morton T., Aigrain S., Foreman-Mackey D., Rajpaul V., 2018, MNRAS, 474, 2094

* Astropy Collaboration, 2013, A&A, 558, A33

* Astropy Collaboration, 2018, AJ, 156, 123

* Astropy Collaboration, 2022, ApJ, 935, 167

* Bellotti S. et al., 2023, A&A, 676, A56

* Cook N. J. et al., 2022, Astrophysics Source Code Library, record ascl:2211.019
Cristofari P. I. et al., 2022, MNRAS, 516, 3802

Delchambre L., 2015, MNRAS, 446, 3545

Donati J.-F., Brown S. F., 1997, A&A, 326, 1135

Donati J.-F., Landstreet J., 2009, ARA&A, 47, 333

Donati J.-F., Semel M., Carter B. D., Rees D. E., Collier Cameron A., 1997, MNRAS, 291, 658

Donati J.-F., Forveille T., Collier Cameron A., Barnes J. R., Delfosse X., Jardine M. M., Valenti J. A., 2006, Science, 311, 633

Donati J.-F. et al., 2008, MNRAS, 390, 545

Donati J. F. et al., 2020, MNRAS, 498, 5684

Donati J. F., Lehmann L. T., Cristofari P. I., Fouqué P., Moutou M., Charpentier P., Ould-Elhkim M., SLS consortium, 2023, MNRAS, 525, 2015(D23)

Finociety B., Donati J. F., 2022, MNRAS, 516, 5887

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASF, 125, 306

Fouqué P. et al., 2023, A&A, 672, A52

Harris C. R. et al., 2020, Nature, 585, 357

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Kavanagh R. D., Vidotto A. A., Klein B., Jardine M. M., Donati J.-F., Ó Fionnagáin D., 2021, MNRAS, 504, 1511

Klein B. et al., 2021, MNRAS, 500, 1844

Kochukhov O., 2021, A&AR, 29, 1

Landi Degl'Innocenti E., Landolfi M., 2004, in Polarization in Spectral Lines, Vol. 307, University of Fierenze, Kluwer Academic Publishers, Dordrecht

Lehmann L. T., Donati J. F., 2022, MNRAS, 514, 2333

Martioli E. et al., 2022, A&A, 660, A86

Medina A. A., Winters J. G., Irwin J. M., Charbonneau D., 2022, ApJ, 935, 104

Morin J. et al., 2008, MNRAS, 390, 567

Morin J., Donati J.-F., Petit P., Delfosse X., Forveille T., Jardine M. M., 2010, MNRAS, 407, 2269

Moutou C. et al., 2017, MNRAS, 472, 4563

Noyes R. W., Hartmann L. W., Baliunas S. L., Duncan D. K., Vaughan A. H., 1984, ApJ, 279, 763

Reiners A. et al., 2022, A&A, 662, A41

Route M., 2016, ApJ, 830, L27

Ryabchikova T., Piskunov N., Kurucz R. L., Stempels H. C., Heiter U., Pakhomov Y., Barklem P. S., 2015, Phys. Scr., 90, 054005

Saar S. H., 1996, in Strassmeier K. G., Linsky J. L., eds, Proc. IAU Symp. 176, Stellar Surface Structure. Kluwer, Dordrecht, p. 237

See V. et al., 2015, MNRAS, 453, 4301

Skilling J., Bryan R. K., 1984, MNRAS, 211, 111

Vidotto A. A. et al., 2014, MNRAS, 441, 2361

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Wright N. J., Drake J. J., Mamajek E. E., Henry G. W., 2011, ApJ, 743, 48

Wright N. J., Newton E. R., Williams P. K. G., Drake J. J., Yadav R. K., 2018, MNRAS, 479, 2351

<figure><img src="image_24.png" /><figcaption>**Figure A1.** The posterior density resulting from the MCMC sampling of the QP best-fitting GPR model for $c_1$ (top) and $B_ℓ$ (bottom) for Gl 905. The concentric circle within each panel indicate the $1\sigma$, $2\sigma$, and $3\sigma$ contours of the distribution. The blue solid lines mark the mode (maximum) of the posterior distribution and the black dashed lines the median and the 16 per cent and 84 per cent percentiles of the posterior probability density function (PDF).</figcaption></figure>

## APPENDIX A: ADDITIONAL GPR FITS
INFORMATION AND FIGURES

Table A1 presents the ephemeris used for our M dwarf sample.

For the GPR fits we used uniform $U$, normal $N$, and Jeffreys $J$ prior distributions. Table A2 indicates the prior used for $P_{\text{rot}}$ and the decay time $l$. For the smoothing factor, we use in general $U(0.1, 1.0)$ beside for the GPR fits of GJ 1151 ($U(0.2, 0.6)$) and Gl 617B ($U(0.4, 0.8)$). For the amplitude $\alpha$, we applied $U(0, \infty)$ with an exception for the GPR fits of GJ 1289 ($N(25, 10)$) and GJ 1286 ($N(16.4, 5)$). For the white noise, we used in general $U(0, \infty)$ beside for Gl 905, where we applied a Jeffreys prior distribution $J(1.5, \infty)$ for the GP fits.

In following, we present the corner plots of all QP GPR fits applied in this paper (see Fig. A1–A6).
Cristofari P. I. et al., 2022, MNRAS, 516, 3802

Delchambre L., 2015, MNRAS, 446, 3545

Donati J.-F., Brown S. F., 1997, A&A, 326, 1135

Donati J.-F., Landstreet J., 2009, ARA&A, 47, 333

Donati J.-F., Semel M., Carter B. D., Rees D. E., Collier Cameron A., 1997, MNRAS, 291, 658

Donati J.-F., Forveille T., Collier Cameron A., Barnes J. R., Delfosse X., Jardine M. M., Valenti J. A., 2006, Science, 311, 633

Donati J.-F. et al., 2008, MNRAS, 390, 545

Donati J. F. et al., 2020, MNRAS, 498, 5684

Donati J. F., Lehmann L. T., Cristofari P. I., Fouqué P., Moutou M., Charpentier P., Ould-Elhkim M., SLS consortium, 2023, MNRAS, 525, 2015(D23)

Finociety B., Donati J. F., 2022, MNRAS, 516, 5887

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASF, 125, 306

Fouqué P. et al., 2023, A&A, 672, A52

Harris C. R. et al., 2020, Nature, 585, 357

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Kavanagh R. D., Vidotto A. A., Klein B., Jardine M. M., Donati J.-F., Ó Fionnagáin D., 2021, MNRAS, 504, 1511

Klein B. et al., 2021, MNRAS, 500, 1844

Kochukhov O., 2021, A&AR, 29, 1

Landi Degl'Innocenti E., Landolfi M., 2004, in Polarization in Spectral Lines, Vol. 307, University of Fierenze, Kluwer Academic Publishers, Dordrecht

Lehmann L. T., Donati J. F., 2022, MNRAS, 514, 2333

Martioli E. et al., 2022, A&A, 660, A86

Medina A. A., Winters J. G., Irwin J. M., Charbonneau D., 2022, ApJ, 935, 104

Morin J. et al., 2008, MNRAS, 390, 567

Morin J., Donati J.-F., Petit P., Delfosse X., Forveille T., Jardine M. M., 2010, MNRAS, 407, 2269

Moutou C. et al., 2017, MNRAS, 472, 4563

Noyes R. W., Hartmann L. W., Baliunas S. L., Duncan D. K., Vaughan A. H., 1984, ApJ, 279, 763

Reiners A. et al., 2022, A&A, 662, A41

Route M., 2016, ApJ, 830, L27

Ryabchikova T., Piskunov N., Kurucz R. L., Stempels H. C., Heiter U., Pakhomov Y., Barklem P. S., 2015, Phys. Scr., 90, 054005

Saar S. H., 1996, in Strassmeier K. G., Linsky J. L., eds, Proc. IAU Symp. 176, Stellar Surface Structure. Kluwer, Dordrecht, p. 237

See V. et al., 2015, MNRAS, 453, 4301

Skilling J., Bryan R. K., 1984, MNRAS, 211, 111

Vidotto A. A. et al., 2014, MNRAS, 441, 2361

Virtanen P. et al., 2020, Nat. Methods, 17, 261

Wright N. J., Drake J. J., Mamajek E. E., Henry G. W., 2011, ApJ, 743, 48

Wright N. J., Newton E. R., Williams P. K. G., Drake J. J., Yadav R. K., 2018, MNRAS, 479, 2351

<figure><img src="image_23.png" /><figcaption>**Figure A1.** The posterior density resulting from the MCMC sampling of the QP best-fitting GPR model for $c_1$ (top) and $B_ℓ$ (bottom) for Gl 905. The concentric circle within each panel indicate the $1\sigma$, $2\sigma$, and $3\sigma$ contours of the distribution. The blue solid lines mark the mode (maximum) of the posterior distribution and the black dashed lines the median and the 16 per cent and 84 per cent percentiles of the posterior probability density function (PDF).</figcaption></figure>

## APPENDIX A: ADDITIONAL GPR FITS
INFORMATION AND FIGURES

Table A1 presents the ephemeris used for our M dwarf sample.

For the GPR fits we used uniform $U$, normal $N$, and Jeffreys $J$ prior distributions. Table A2 indicates the prior used for $P_{\text{rot}}$ and the decay time $l$. For the smoothing factor, we use in general $U(0.1, 1.0)$ beside for the GPR fits of GJ 1151 ($U(0.2, 0.6)$) and Gl 617B ($U(0.4, 0.8)$). For the amplitude $\alpha$, we applied $U(0, \infty)$ with an exception for the GPR fits of GJ 1289 ($N(25, 10)$) and GJ 1286 ($N(16.4, 5)$). For the white noise, we used in general $U(0, \infty)$ beside for Gl 905, where we applied a Jeffreys prior distribution $J(1.5, \infty)$ for the GP fits.

In following, we present the corner plots of all QP GPR fits applied in this paper (see Fig. A1–A6).
<figure><img src="image_25.png" /><figcaption>Figure A2. Same as Fig. A1 for GJ 1289.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure A3. Same as Fig. A1 for GJ 1151.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure A2. Same as Fig. A1 for GJ 1289.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure A3. Same as Fig. A1 for GJ 1151.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure A4. Same as Fig. A1 for GJ 1286.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure A5. Same as Fig. A1 for Gl 617B.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure A4.** Same as Fig. A1 for GJ 1286.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure A5.** Same as Fig. A1 for GJ 617B.</figcaption></figure>
<table><caption>**Table A1.** The ephemeris data used to determine the phases (phase = $T_0 - T_{obs}/P_{rot}$) for the six M dwarfs of our sample. For $T_0$ we used the barycentric Julian date of the first SPIRou observation of the target (second column). The last column indicates the rotation period $P_{rot}$ used for the phase determination.</caption><thead><tr><td>Star</td><td>$T_0$ (JD)</td><td>$P_{rot}$ (d)</td></tr></thead><tbody><tr><td>Gl 905</td><td>2458600.1348284</td><td>114.3</td></tr><tr><td>GJ 1289</td><td>2458648.9872293</td><td>73.66</td></tr><tr><td>GJ 1151</td><td>2458828.1234878</td><td>175.6</td></tr><tr><td>GJ 1286</td><td>2458745.8452793</td><td>178.0</td></tr><tr><td>Gl 617B</td><td>2458744.743105</td><td>40.4</td></tr><tr><td>Gl 408</td><td>2458590.0269338</td><td>171.0</td></tr></tbody></table>

<table><caption>**Table A2.** The prior distribution used for the QP GPR fits for the rotation period $P_{rot}$ and decay time $l$.</caption><thead><tr><td>Star</td><td>Prior $P_{rot}$ (d)</td><td>Prior $l$ (d)</td></tr></thead><tbody><tr><td>Gl 905</td><td>$\mathcal{N}$(115, 20)</td><td>$\mathcal{N}$(130, 25)</td></tr><tr><td>GJ 1289</td><td>$\mathcal{N}$(75, 5)</td><td>$\mathcal{U}$(50, 1000)</td></tr><tr><td>GJ 1151</td><td>$\mathcal{N}$(175, 20)</td><td>300</td></tr><tr><td>GJ 1286</td><td>$\mathcal{N}$(170, 20)</td><td>300</td></tr><tr><td>Gl 617B</td><td>$\mathcal{N}$(40.4, 5)</td><td>$\mathcal{U}$(30, 1000)</td></tr><tr><td>Gl 408</td><td>$\mathcal{N}$(170, 15)</td><td>200</td></tr></tbody></table>

## APPENDIX B: THE GPR FITS OF THE LONGITUDINAL FIELD VALUES

Table B1 provides a comparison of the GP results for $c_1$ and $B_\ell$ obtained with the GPR framework presented in Section 3.2 and the GP results of D23 for the $B_\ell$ values (marked by D23 in the second column). For the GP fits the mean value is fixed to the $B_\ell$ mean values determined by D23 for both $c_1$ and $B_\ell$. Within the error, all three GP results agree for $P_{rot}$, $l$, $\beta$, and $\alpha$ for each M dwarf, confirming $P_{rot}$, $l$, $\beta$, and $\alpha$ using two different variables and calculation routines with one exception for Gl 617B, where the decay time determined from $c_1$ differs from the GP results of $B_\ell$. The additional white noise $\sigma$ is also consistent for both $B_\ell$, while $c_1$ often shows lower $\sigma$ values. Our GP routines fit the $B_\ell$ data with a slightly higher or equal $\chi^2$ compared to the results of D23, while obtaining lower $\sigma$ values. Comparing the results for $c_1$ and $B_\ell$, we find that $c_1$ has a lower rms. The rms of $c_1$ and $B_\ell$ is for all M dwarfs lower than the corresponding averaged error of $c_1$ and $B_\ell$. For M dwarfs with significant non-axisymmetric field (e.g. Gl 905, GJ 1289, GJ 1286), $c_1$ is fitted with lower or equal $\chi^2$ and can provide smaller errors for $P_{rot}$ (e.g. GJ 1289, GJ 1286). For M dwarfs with highly axisymmetric fields (e.g. Gl 617B, Gl 408), $B_\ell$ gives smaller errors on $P_{rot}$.

<figure><img src="image_29.png" /><figcaption>**Figure A6.** Same as Fig. A1 for Gl 408.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Figure A6.** Same as Fig. A1 for Gl 408.</figcaption></figure>

<table><caption>**Table A1.** The ephemeris data used to determine the phases (phase = $T_0 - T_{obs}/P_{rot}$) for the six M dwarfs of our sample. For $T_0$ we used the barycentric Julian date of the first SPIRou observation of the target (second column). The last column indicates the rotation period $P_{rot}$ used for the phase determination.</caption><thead><tr><td>Star</td><td>$T_0$ (JD)</td><td>$P_{rot}$ (d)</td></tr></thead><tbody><tr><td>Gl 905</td><td>2458600.1348284</td><td>114.3</td></tr><tr><td>GJ 1289</td><td>2458648.9872293</td><td>73.66</td></tr><tr><td>GJ 1151</td><td>2458828.1234878</td><td>175.6</td></tr><tr><td>GJ 1286</td><td>2458745.8452793</td><td>178.0</td></tr><tr><td>Gl 617B</td><td>2458744.743105</td><td>40.4</td></tr><tr><td>Gl 408</td><td>2458590.0269338</td><td>171.0</td></tr></tbody></table>

<table><caption>**Table A2.** The prior distribution used for the QP GPR fits for the rotation period $P_{rot}$ and decay time $l$.</caption><thead><tr><td>Star</td><td>Prior $P_{rot}$ (d)</td><td>Prior $l$ (d)</td></tr></thead><tbody><tr><td>Gl 905</td><td>$\mathcal{N}$(115, 20)</td><td>$\mathcal{N}$(130, 25)</td></tr><tr><td>GJ 1289</td><td>$\mathcal{N}$(75, 5)</td><td>$\mathcal{U}$(50, 1000)</td></tr><tr><td>GJ 1151</td><td>$\mathcal{N}$(175, 20)</td><td>300</td></tr><tr><td>GJ 1286</td><td>$\mathcal{N}$(170, 20)</td><td>300</td></tr><tr><td>Gl 617B</td><td>$\mathcal{N}$(40.4, 5)</td><td>$\mathcal{U}$(30, 1000)</td></tr><tr><td>Gl 408</td><td>$\mathcal{N}$(170, 15)</td><td>200</td></tr></tbody></table>

## APPENDIX B: THE GPR FITS OF THE LONGITUDINAL FIELD VALUES

Table B1 provides a comparison of the GP results for $c_1$ and $B_\ell$ obtained with the GPR framework presented in Section 3.2 and the GP results of D23 for the $B_\ell$ values (marked by D23 in the second column). For the GP fits the mean value is fixed to the $B_\ell$ mean values determined by D23 for both $c_1$ and $B_\ell$. Within the error, all three GP results agree for $P_{rot}$, $l$, $\beta$, and $\alpha$ for each M dwarf, confirming $P_{rot}$, $l$, $\beta$, and $\alpha$ using two different variables and calculation routines with one exception for Gl 617B, where the decay time determined from $c_1$ differs from the GP results of $B_\ell$. The additional white noise $\sigma$ is also consistent for both $B_\ell$, while $c_1$ often shows lower $\sigma$ values. Our GP routines fit the $B_\ell$ data with a slightly higher or equal $\chi^2$ compared to the results of D23, while obtaining lower $\sigma$ values. Comparing the results for $c_1$ and $B_\ell$, we find that $c_1$ has a lower rms. The rms of $c_1$ and $B_\ell$ is for all M dwarfs lower than the corresponding averaged error of $c_1$ and $B_\ell$. For M dwarfs with significant non-axisymmetric field (e.g. Gl 905, GJ 1289, GJ 1286), $c_1$ is fitted with lower or equal $\chi^2$ and can provide smaller errors for $P_{rot}$ (e.g. GJ 1289, GJ 1286). For M dwarfs with highly axisymmetric fields (e.g. Gl 617B, Gl 408), $B_\ell$ gives smaller errors on $P_{rot}$.
<table><caption>**Table B1.** Summary of the best-fitting parameters of the QP GPR fits applied to $c_1$ and $B_ℓ$ for the six M dwarfs in our sample and the comparison with the corresponding GPR results for $B_ℓ$ by D23 marked as D23 in the second column.</caption><thead><tr><th>Star</th><th>Data</th><th>Rotation period<br>P<sub>rot</sub> (d)</th><th>Decay time<br>l (d)</th><th>Smoothing factor<br>β</th><th>Amplitude<br>α (G)</th><th>White noise<br>σ (G)</th><th>Rms<br>(G)</th><th>χ<sup>2</sup><sub>r</sub></th></tr></thead><tbody><tr><td>Gl 905</td><td>c<sub>1</sub></td><td>111.7<sup>+3.0</sup><sub>-3.2</sub></td><td>133<sup>+18</sup><sub>-22</sub></td><td>0.50<sup>+0.09</sup><sub>-0.07</sub></td><td>12.9<sup>+3.1</sup><sub>-2.1</sub></td><td>0.6<sup>+0.5</sup><sub>-0.4</sub></td><td>3.8</td><td>0.79</td></tr><tr><td>Gl 905</td><td>B<sub>ℓ</sub></td><td>114.4<sup>+3.5</sup><sub>-2.4</sub></td><td>130<sup>+25</sup><sub>-32</sub></td><td>0.44<sup>+0.12</sup><sub>-0.09</sub></td><td>14.1<sup>+3.5</sup><sub>-2.4</sub></td><td>1.2<sup>+1.1</sup><sub>-0.8</sub></td><td>6.9</td><td>0.87</td></tr><tr><td>Gl 905</td><td>D23</td><td>114.3 ± 2.8</td><td>129<sup>+25</sup><sub>-21</sub></td><td>0.43 ± 0.09</td><td>13.3<sup>+2.5</sup><sub>-2.1</sub></td><td>1.7<sup>+1.2</sup><sub>-0.7</sub></td><td>6.2</td><td>0.84</td></tr><tr><td>GJ 1289</td><td>c<sub>1</sub></td><td>75.62<sup>+0.85</sup><sub>-0.79</sub></td><td>129<sup>+26</sup><sub>-24</sub></td><td>0.40 ± 0.05</td><td>41.9<sup>+7.9</sup><sub>-6.0</sub></td><td>1.2<sup>+1.1</sup><sub>-0.8</sub></td><td>6.9</td><td>0.72</td></tr><tr><td>GJ 1289</td><td>B<sub>ℓ</sub></td><td>73.67<sup>+1.01</sup><sub>-0.91</sub></td><td>152<sup>+30</sup><sub>-28</sub></td><td>0.48 ± 0.08</td><td>53.3<sup>+11.7</sup><sub>-9.4</sub></td><td>3.2<sup>+2.4</sup><sub>-2.1</sub></td><td>14.3</td><td>0.84</td></tr><tr><td>GJ 1289</td><td>D23</td><td>73.66 ± 0.92</td><td>152<sup>+32</sup><sub>-27</sub></td><td>0.48 ± 0.09</td><td>53.2<sup>+12.4</sup><sub>-10.1</sub></td><td>4.2<sup>+2.7</sup><sub>-1.6</sub></td><td>13.9</td><td>0.82</td></tr><tr><td>GJ 1151</td><td>c<sub>1</sub></td><td>175.8<sup>+3.2</sup><sub>-3.4</sub></td><td>300</td><td>0.40<sup>+0.10</sup><sub>-0.09</sub></td><td>12.7<sup>+3.2</sup><sub>-2.4</sub></td><td>3.3 ± 0.7</td><td>5.8</td><td>0.99</td></tr><tr><td>GJ 1151</td><td>B<sub>ℓ</sub></td><td>176.1<sup>+3.6</sup><sub>-4.1</sub></td><td>300</td><td>0.42<sup>+0.10</sup><sub>-0.09</sub></td><td>14.6<sup>+3.4</sup><sub>-2.7</sub></td><td>1.1<sup>+1.1</sup><sub>-0.8</sub></td><td>6.7</td><td>0.73</td></tr><tr><td>GJ 1151</td><td>D23</td><td>175.6 ± 4.9</td><td>300</td><td>0.43 ± 0.11</td><td>14.9<sup>+4.2</sup><sub>-3.3</sub></td><td>1.6<sup>+1.3</sup><sub>-0.7</sub></td><td>6.7</td><td>0.72</td></tr><tr><td>GJ 1286</td><td>c<sub>1</sub></td><td>186.8<sup>+9.5</sup><sub>-5.8</sub></td><td>300</td><td>0.23<sup>+0.05</sup><sub>-0.04</sub></td><td>18.0<sup>+4.3</sup><sub>-3.1</sub></td><td>2.7<sup>+1.5</sup><sub>-1.7</sub></td><td>7.5</td><td>1.02</td></tr><tr><td>GJ 1286</td><td>B<sub>ℓ</sub></td><td>181<sup>+18</sup><sub>-13</sub></td><td>300</td><td>0.28<sup>+0.13</sup><sub>-0.07</sub></td><td>17.2<sup>+6.2</sup><sub>-3.6</sub></td><td>3.9<sup>+2.3</sup><sub>-2.4</sub></td><td>12.0</td><td>1.03</td></tr><tr><td>GJ 1286</td><td>D23</td><td>178 ± 15</td><td>300</td><td>0.29 ± 0.09</td><td>16.7<sup>+4.6</sup><sub>-3.6</sub></td><td>4.6<sup>+2.6</sup><sub>-1.6</sub></td><td>10.1</td><td>1.02</td></tr><tr><td>Gl 617B</td><td>c<sub>1</sub></td><td>37.8<sup>+8.5</sup><sub>-2.6</sub></td><td>35<sup>+8</sup><sub>-4</sub></td><td>0.47<sup>+0.09</sup><sub>-0.05</sub></td><td>5.9<sup>+1.2</sup><sub>-0.8</sub></td><td>0.7<sup>+0.6</sup><sub>-0.5</sub></td><td>2.2</td><td>0.66</td></tr><tr><td>Gl 617B</td><td>B<sub>ℓ</sub></td><td>40.6<sup>+2.1</sup><sub>-4.4</sub></td><td>82<sup>+45</sup><sub>-30</sub></td><td>0.58<sup>+0.13</sup><sub>-0.11</sub></td><td>5.5<sup>+1.3</sup><sub>-1.1</sub></td><td>1.3<sup>+1.0</sup><sub>-0.9</sub></td><td>5.0</td><td>0.86</td></tr><tr><td>Gl 617B</td><td>D23</td><td>40.4 ± 3.0</td><td>69<sup>+35</sup><sub>-23</sub></td><td>0.60 ± 0.22</td><td>5.4<sup>+1.3</sup><sub>-1.0</sub></td><td>1.7<sup>+1.1</sup><sub>-0.7</sub></td><td>4.9</td><td>0.86</td></tr><tr><td>Gl 408</td><td>c<sub>1</sub></td><td>175<sup>+12</sup><sub>-14</sub></td><td>200</td><td>0.18<sup>+0.10</sup><sub>-0.05</sub></td><td>4.2<sup>+0.9</sup><sub>-0.8</sub></td><td>2.2<sup>+0.5</sup><sub>-0.6</sub></td><td>3.9</td><td>1.19</td></tr><tr><td>Gl 408</td><td>B<sub>ℓ</sub></td><td>170.7<sup>+7.1</sup><sub>-9.8</sub></td><td>200</td><td>0.22<sup>+0.13</sup><sub>-0.08</sub></td><td>6.3<sup>+1.6</sup><sub>-1.2</sub></td><td>1.0<sup>+1.0</sup><sub>-0.7</sub></td><td>6.9</td><td>0.72</td></tr><tr><td>Gl 408</td><td>D23</td><td>171.0 ± 8.4</td><td>200</td><td>0.21 ± 0.10</td><td>6.3<sup>+1.5</sup><sub>-1.2</sub></td><td>1.5<sup>+1.2</sup><sub>-0.7</sub></td><td>6.3</td><td>0.66</td></tr></tbody></table>

## APPENDIX C: ADDITIONAL PCA AND ZDI FIGURES

We can confirm that the symmetric component (with respect to the line centre) seen in the mean profile of the 2020/21 season for Gl 905 is due to the irregular phase coverage of the observations and does not reflect an axisymmetric toroidal field. We simulated 24 equally spaced Stokes V LSD profiles from the 2020/21 ZDI map (Fig. 3 middle column). The symmetric component of the resulting mean profile disappears for the uniform phase coverage (see Fig. C1). None the less, we might still miss toroidal field due to the low $v_c \sin i$, see the $\sigma_{(B_{tor})}$ estimation in Section 4.2.

Figs C2 and C3 show the ZDI maps of Gl 617B and Gl 408 for an alternative inclination $i = 30^{\circ}$.

The following figures (Fig. C4–C9) display the observed Stokes V LSD profiles (black) and their ZDI fits (red) split by season for the six M dwarfs.

<figure><img src="image_30.png" /><figcaption>**Figure C1.** The mean profile and its decompositions obtained from 24 uniformly phased synthetic Stokes V LSD profiles of the 2020/21 ZDI map of Gl 905. The same format as in Fig. 1(a) is used.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure C2.** The magnetic field maps of Gl 617B using an inclination $i = 30^{\circ}$ presented in the same format as in Fig. 3.</figcaption></figure>
<table><caption>**Table B1.** Summary of the best-fitting parameters of the QP GPR fits applied to $c_1$ and $B_ℓ$ for the six M dwarfs in our sample and the comparison with the corresponding GPR results for $B_ℓ$ by D23 marked as D23 in the second column.</caption><thead><tr><th>Star</th><th>Data</th><th>Rotation period<br>P<sub>rot</sub> (d)</th><th>Decay time<br>l (d)</th><th>Smoothing factor<br>β</th><th>Amplitude<br>α (G)</th><th>White noise<br>σ (G)</th><th>Rms<br>(G)</th><th>χ<sup>2</sup><sub>r</sub></th></tr></thead><tbody><tr><td>Gl 905</td><td>c<sub>1</sub></td><td>111.7<sup>+3.0</sup><sub>-3.2</sub></td><td>133<sup>+18</sup><sub>-22</sub></td><td>0.50<sup>+0.09</sup><sub>-0.07</sub></td><td>12.9<sup>+3.1</sup><sub>-2.1</sub></td><td>0.6<sup>+0.5</sup><sub>-0.4</sub></td><td>3.8</td><td>0.79</td></tr><tr><td>Gl 905</td><td>B<sub>ℓ</sub></td><td>114.4<sup>+3.5</sup><sub>-2.4</sub></td><td>130<sup>+25</sup><sub>-32</sub></td><td>0.44<sup>+0.12</sup><sub>-0.09</sub></td><td>14.1<sup>+3.5</sup><sub>-2.4</sub></td><td>1.2<sup>+1.1</sup><sub>-0.8</sub></td><td>6.9</td><td>0.87</td></tr><tr><td>Gl 905</td><td>D23</td><td>114.3 ± 2.8</td><td>129<sup>+25</sup><sub>-21</sub></td><td>0.43 ± 0.09</td><td>13.3<sup>+2.5</sup><sub>-2.1</sub></td><td>1.7<sup>+1.2</sup><sub>-0.7</sub></td><td>6.2</td><td>0.84</td></tr><tr><td>GJ 1289</td><td>c<sub>1</sub></td><td>75.62<sup>+0.85</sup><sub>-0.79</sub></td><td>129<sup>+26</sup><sub>-24</sub></td><td>0.40 ± 0.05</td><td>41.9<sup>+7.9</sup><sub>-6.0</sub></td><td>1.2<sup>+1.1</sup><sub>-0.8</sub></td><td>6.9</td><td>0.72</td></tr><tr><td>GJ 1289</td><td>B<sub>ℓ</sub></td><td>73.67<sup>+1.01</sup><sub>-0.91</sub></td><td>152<sup>+30</sup><sub>-28</sub></td><td>0.48 ± 0.08</td><td>53.3<sup>+11.7</sup><sub>-9.4</sub></td><td>3.2<sup>+2.4</sup><sub>-2.1</sub></td><td>14.3</td><td>0.84</td></tr><tr><td>GJ 1289</td><td>D23</td><td>73.66 ± 0.92</td><td>152<sup>+32</sup><sub>-27</sub></td><td>0.48 ± 0.09</td><td>53.2<sup>+12.4</sup><sub>-10.1</sub></td><td>4.2<sup>+2.7</sup><sub>-1.6</sub></td><td>13.9</td><td>0.82</td></tr><tr><td>GJ 1151</td><td>c<sub>1</sub></td><td>175.8<sup>+3.2</sup><sub>-3.4</sub></td><td>300</td><td>0.40<sup>+0.10</sup><sub>-0.09</sub></td><td>12.7<sup>+3.2</sup><sub>-2.4</sub></td><td>3.3 ± 0.7</td><td>5.8</td><td>0.99</td></tr><tr><td>GJ 1151</td><td>B<sub>ℓ</sub></td><td>176.1<sup>+3.6</sup><sub>-4.1</sub></td><td>300</td><td>0.42<sup>+0.10</sup><sub>-0.09</sub></td><td>14.6<sup>+3.4</sup><sub>-2.7</sub></td><td>1.1<sup>+1.1</sup><sub>-0.8</sub></td><td>6.7</td><td>0.73</td></tr><tr><td>GJ 1151</td><td>D23</td><td>175.6 ± 4.9</td><td>300</td><td>0.43 ± 0.11</td><td>14.9<sup>+4.2</sup><sub>-3.3</sub></td><td>1.6<sup>+1.3</sup><sub>-0.7</sub></td><td>6.7</td><td>0.72</td></tr><tr><td>GJ 1286</td><td>c<sub>1</sub></td><td>186.8<sup>+9.5</sup><sub>-5.8</sub></td><td>300</td><td>0.23<sup>+0.05</sup><sub>-0.04</sub></td><td>18.0<sup>+4.3</sup><sub>-3.1</sub></td><td>2.7<sup>+1.5</sup><sub>-1.7</sub></td><td>7.5</td><td>1.02</td></tr><tr><td>GJ 1286</td><td>B<sub>ℓ</sub></td><td>181<sup>+18</sup><sub>-13</sub></td><td>300</td><td>0.28<sup>+0.13</sup><sub>-0.07</sub></td><td>17.2<sup>+6.2</sup><sub>-3.6</sub></td><td>3.9<sup>+2.3</sup><sub>-2.4</sub></td><td>12.0</td><td>1.03</td></tr><tr><td>GJ 1286</td><td>D23</td><td>178 ± 15</td><td>300</td><td>0.29 ± 0.09</td><td>16.7<sup>+4.6</sup><sub>-3.6</sub></td><td>4.6<sup>+2.6</sup><sub>-1.6</sub></td><td>10.1</td><td>1.02</td></tr><tr><td>Gl 617B</td><td>c<sub>1</sub></td><td>37.8<sup>+8.5</sup><sub>-2.6</sub></td><td>35<sup>+8</sup><sub>-4</sub></td><td>0.47<sup>+0.09</sup><sub>-0.05</sub></td><td>5.9<sup>+1.2</sup><sub>-0.8</sub></td><td>0.7<sup>+0.6</sup><sub>-0.5</sub></td><td>2.2</td><td>0.66</td></tr><tr><td>Gl 617B</td><td>B<sub>ℓ</sub></td><td>40.6<sup>+2.1</sup><sub>-4.4</sub></td><td>82<sup>+45</sup><sub>-30</sub></td><td>0.58<sup>+0.13</sup><sub>-0.11</sub></td><td>5.5<sup>+1.3</sup><sub>-1.1</sub></td><td>1.3<sup>+1.0</sup><sub>-0.9</sub></td><td>5.0</td><td>0.86</td></tr><tr><td>Gl 617B</td><td>D23</td><td>40.4 ± 3.0</td><td>69<sup>+35</sup><sub>-23</sub></td><td>0.60 ± 0.22</td><td>5.4<sup>+1.3</sup><sub>-1.0</sub></td><td>1.7<sup>+1.1</sup><sub>-0.7</sub></td><td>4.9</td><td>0.86</td></tr><tr><td>Gl 408</td><td>c<sub>1</sub></td><td>175<sup>+12</sup><sub>-14</sub></td><td>200</td><td>0.18<sup>+0.10</sup><sub>-0.05</sub></td><td>4.2<sup>+0.9</sup><sub>-0.8</sub></td><td>2.2<sup>+0.5</sup><sub>-0.6</sub></td><td>3.9</td><td>1.19</td></tr><tr><td>Gl 408</td><td>B<sub>ℓ</sub></td><td>170.7<sup>+7.1</sup><sub>-9.8</sub></td><td>200</td><td>0.22<sup>+0.13</sup><sub>-0.08</sub></td><td>6.3<sup>+1.6</sup><sub>-1.2</sub></td><td>1.0<sup>+1.0</sup><sub>-0.7</sub></td><td>6.9</td><td>0.72</td></tr><tr><td>Gl 408</td><td>D23</td><td>171.0 ± 8.4</td><td>200</td><td>0.21 ± 0.10</td><td>6.3<sup>+1.5</sup><sub>-1.2</sub></td><td>1.5<sup>+1.2</sup><sub>-0.7</sub></td><td>6.3</td><td>0.66</td></tr></tbody></table>

## APPENDIX C: ADDITIONAL PCA AND ZDI FIGURES

We can confirm that the symmetric component (with respect to the line centre) seen in the mean profile of the 2020/21 season for Gl 905 is due to the irregular phase coverage of the observations and does not reflect an axisymmetric toroidal field. We simulated 24 equally spaced Stokes V LSD profiles from the 2020/21 ZDI map (Fig. 3 middle column). The symmetric component of the resulting mean profile disappears for the uniform phase coverage (see Fig. C1). None the less, we might still miss toroidal field due to the low $v_c \sin i$, see the $\sigma_{(B_{tor})}$ estimation in Section 4.2.

Figs C2 and C3 show the ZDI maps of Gl 617B and Gl 408 for an alternative inclination $i = 30^{\circ}$.

The following figures (Fig. C4–C9) display the observed Stokes V LSD profiles (black) and their ZDI fits (red) split by season for the six M dwarfs.

<figure><img src="image_29.png" /><figcaption>**Figure C1.** The mean profile and its decompositions obtained from 24 uniformly phased synthetic Stokes V LSD profiles of the 2020/21 ZDI map of Gl 905. The same format as in Fig. 1(a) is used.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure C2.** The magnetic field maps of Gl 617B using an inclination $i = 30^{\circ}$ presented in the same format as in Fig. 3.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Figure C3.** The magnetic field maps of Gl 408 using an inclination $i = 30^\circ$ presented in the same format as in Fig. 3.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Figure C4.** The SPIRou observed LSD Stokes V 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.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Figure C3.** The magnetic field maps of Gl 408 using an inclination $i = 30^\circ$ presented in the same format as in Fig. 3.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure C4.** The SPIRou observed LSD Stokes V 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.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure C5. Same as Fig. C4 for GJ 1289.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Figure C6. Same as Fig. C4 for GJ 1151.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Figure C5. Same as Fig. C4 for GJ 1289.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure C6. Same as Fig. C4 for GJ 1151.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>Figure C7. Same as Fig. C4 for GJ 1286.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Figure C8. Same as Fig. C4 for Gl 617B.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure C7. Same as Fig. C4 for GJ 1286.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure C8. Same as Fig. C4 for Gl 617B.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>Figure C9. Same as Fig. C4 for Gl 408.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.
<figure><img src="image_37.png" /><figcaption>Figure C9. Same as Fig. C4 for Gl 408.</figcaption></figure>

This paper has been typeset from a TEX/LATEX file prepared by the author.