# Long-term monitoring of large-scale magnetic fields across optical and near-infrared domains with ESPaDOnS, Narval, and SPIRou

S. Bellotti, J. Morin, L T Lehmann, P. Petit, G a J Hussain, J.-F. Donati, C

P Folsom, A. Carmona, E. Martioli, B. Klein, et al.

► To cite this version:

S. Bellotti, J. Morin, L T Lehmann, P. Petit, G a J Hussain, et al., Long-term monitoring of large-scale magnetic fields across optical and near-infrared domains with ESPaDOnS, Narval, and SPIRou. Astronomy and Astrophysics - A&A, 2024, 686, pp.A66. 10.1051/0004-6361/202348043. hal-04592914

HAL Id: hal-04592914

https://hal.science/hal-04592914

Submitted on 29 May 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
# Long-term monitoring of large-scale magnetic fields across optical and near-infrared domains with ESPaDOnS, Narval, and SPIRou

S. Bellotti, J. Morin, L T Lehmann, P. Petit, G a J Hussain, J.-F. Donati, C

P Folsom, A. Carmona, E. Martioli, B. Klein, et al.

► To cite this version:

S. Bellotti, J. Morin, L T Lehmann, P. Petit, G a J Hussain, et al.. Long-term monitoring of large-scale magnetic fields across optical and near-infrared domains with ESPaDOnS, Narval, and SPIRou. Astronomy & Astrophysics - A&A, 2024, 686, pp.A66. 10.1051/0004-6361/202348043. hal-04592914

HAL Id: hal-04592914

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

Submitted on 29 May 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
# Long-term monitoring of large-scale magnetic fields across optical and near-infrared domains with ESPaDOnS, Narval, and SPIRou

## The cases of EV Lac, DS Leo, and CN Leo*

S. Bellotti<sup>1,2 Ⓐ</sup>, J. Morin<sup>3 Ⓐ</sup>, L. T. Lehmann<sup>1 Ⓐ</sup>, P. Petit<sup>1 Ⓐ</sup>, G. A. J. Hussain<sup>2 Ⓐ</sup>, J.-F. Donati<sup>1 Ⓐ</sup>, C. P. Folsom<sup>4 Ⓐ</sup>, A. Carmona<sup>5 Ⓐ</sup>, E. Martioli<sup>6,7 Ⓐ</sup>, B. Klein<sup>8 Ⓐ</sup>, P. Fouqué<sup>1 Ⓐ</sup>, C. Moutou<sup>1 Ⓐ</sup>, S. Alencar<sup>9</sup>, E. Artigau<sup>10 Ⓐ</sup>, I. Boisse<sup>11 Ⓐ</sup>, F. Bouchy<sup>12 Ⓐ</sup>, J. Bouvier<sup>5 Ⓐ</sup>, N. J. Cook<sup>10 Ⓐ</sup>, X. Delfosse<sup>5 Ⓐ</sup>, R. Doyon<sup>10 Ⓐ</sup>, and G. Hébrard<sup>7,13 Ⓐ</sup>

Received 22 September 2023 / Accepted 5 March 2024

## ABSTRACT

**Context.** Dynamo models describing the generation of stellar magnetic fields for partly and fully convective stars are guided by observational constraints. Zeeman-Doppler imaging has revealed a variety of magnetic field geometries and, for fully convective stars in particular, a dichotomy: either strong, mostly axisymmetric, and dipole-dominated or weak, non-axisymmetric, and multipole-dominated. This dichotomy is explained either by dynamo bistability (i.e., two coexisting and stable dynamo branches) or by long-term magnetic cycles with polarity reversals, but there is no definite conclusion on the matter.

**Aims.** Our aim is to monitor the evolution of the large-scale field for a sample of nearby M dwarfs with masses between 0.1 and 0.6 $M_\odot$, which is of prime interest to inform distinct dynamo theories and explain the variety of magnetic field geometries studied in previous works. This also has the potential to put long-term cyclic variations of the Sun's magnetic field into a broader context.

**Methods.** We analysed optical spectropolarimetric data sets collected with ESPaDOnS and Narval between 2005 and 2016, and near-infrared SPIRou data obtained between 2019 and 2022 for three well-studied, active M dwarfs: EV Lac, DS Leo, and CN Leo. We looked for secular changes in time series of longitudinal magnetic field, width of unpolarised mean-line profiles, and large-scale field topology as retrieved with principal component analysis and Zeeman-Doppler imaging.

**Results.** We retrieved pulsating (EV Lac), stable (DS Leo), and sine-like (CN Leo) long-term trends in longitudinal field. The width of near-infrared mean-line profiles exhibits rotational modulation only for DS Leo, whereas in the optical it is evident for both EV Lac and DS Leo. The line width variations are not necessarily correlated to those of the longitudinal field, suggesting complex relations between small- and large-scale field. We also recorded topological changes in the form of a reduced axisymmetry for EV Lac and transition from a toroidal-dominated to poloidal-dominated regime for DS Leo. For CN Leo, the topology remained predominantly poloidal, dipolar, and axisymmetric, with only an oscillation in field strength.

**Conclusions.** Our results show a peculiar evolution of the magnetic field for each M dwarf individually, with DS Leo and EV Lac manifesting more evident variations than CN Leo. These findings confirm that M dwarfs with distinct masses and rotation periods can undergo magnetic long-term variations and suggest an underlying variety of cyclic behaviours of their magnetic fields.

**Key words.** techniques: polarimetric – stars: activity – stars: magnetic field – stars: individual: EV Lac – stars: individual: DS Leo – stars: individual: CN Leo

## 1. Introduction

The magnetic fields of low-mass stars ($M_* < 1.2 M_\odot$) are powered by dynamos (Schrijver & Zwaan 2000), and their study is

paramount to understanding stellar evolution and activity phenomena. The evolution of a star's rotation is linked to magnetic activity as the magnetic field couples with the stellar wind and results in angular momentum loss over the star's lifetime (Skumanich 1972; Vidotto et al. 2014; Finley & Matt 2017; See et al. 2019). Therefore, one can infer the stellar age and rotation period using magnetic activity as a proxy (Noyes et al.

<footnote n="1"><sup>1</sup> Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, IRAP/UMR 5277, 14 Avenue Edouard Belin, 31400 Toulouse, France
e-mail: stefano.bellotti@irap.omp.eu</footnote>

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

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

<footnote n="4"><sup>4</sup> Tartu Observatory, University of Tartu, Observatooriumi 1, Töravere 61602 Tartumaa, Estonia</footnote>

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

<footnote n="6"><sup>6</sup> Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, 37504-364 Itajubá, MG, Brazil</footnote>

<footnote n="7"><sup>7</sup> Institut d'Astrophysique de Paris, CNRS, UMR 7095, Sorbonne Université, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="8"><sup>8</sup> Sub-department of Astrophysics, Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

<footnote n="9"><sup>9</sup> Universidade Federal de Minas Gerais, Belo Horizonte, MG 31270-901, Brazil</footnote>

<footnote n="10"><sup>10</sup> Université de Montréal, Département de Physique, IREX, Montréal, QC H3C 3J7, Canada</footnote>

<footnote n="11"><sup>11</sup> Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="12"><sup>12</sup> Observatoire de Genève, Université de Genève, Chemin Pegasi, 51, 1290 Sauverny, Switzerland</footnote>

<footnote n="13"><sup>13</sup> Observatoire de Haute Provence, St Michel l'Observatoire, France</footnote>

<footnote>* The reduced data are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/686/A66</footnote>
# Long-term monitoring of large-scale magnetic fields across optical and near-infrared domains with ESPaDOnS, Narval, and SPIRou

## The cases of EV Lac, DS Leo, and CN Leo*

S. Bellotti<sup>1,2 Ⓓ</sup>, J. Morin<sup>3 Ⓓ</sup>, L. T. Lehmann<sup>1 Ⓓ</sup>, P. Petit<sup>1 Ⓓ</sup>, G. A. J. Hussain<sup>2 Ⓓ</sup>, J.-F. Donati<sup>1 Ⓓ</sup>, C. P. Folsom<sup>4 Ⓓ</sup>, A. Carmona<sup>5 Ⓓ</sup>, E. Martioli<sup>6,7 Ⓓ</sup>, B. Klein<sup>8 Ⓓ</sup>, P. Fouqué<sup>1 Ⓓ</sup>, C. Moutou<sup>1 Ⓓ</sup>, S. Alencar<sup>9</sup>, E. Artigau<sup>10 Ⓓ</sup>, I. Boisse<sup>11 Ⓓ</sup>, F. Bouchy<sup>12 Ⓓ</sup>, J. Bouvier<sup>5 Ⓓ</sup>, N. J. Cook<sup>10 Ⓓ</sup>, X. Delfosse<sup>5 Ⓓ</sup>, R. Doyon<sup>10 Ⓓ</sup>, and G. Hébrard<sup>7,13 Ⓓ</sup>

Received 22 September 2023 / Accepted 5 March 2024

## ABSTRACT

**Context.** Dynamo models describing the generation of stellar magnetic fields for partly and fully convective stars are guided by observational constraints. Zeeman-Doppler imaging has revealed a variety of magnetic field geometries and, for fully convective stars in particular, a dichotomy: either strong, mostly axisymmetric, and dipole-dominated or weak, non-axisymmetric, and multipole-dominated. This dichotomy is explained either by dynamo bistability (i.e., two coexisting and stable dynamo branches) or by long-term magnetic cycles with polarity reversals, but there is no definite conclusion on the matter.

**Aims.** Our aim is to monitor the evolution of the large-scale field for a sample of nearby M dwarfs with masses between 0.1 and 0.6 $M_\odot$, which is of prime interest to inform distinct dynamo theories and explain the variety of magnetic field geometries studied in previous works. This also has the potential to put long-term cyclic variations of the Sun's magnetic field into a broader context.

**Methods.** We analysed optical spectropolarimetric data sets collected with ESPaDOnS and Narval between 2005 and 2016, and near-infrared SPIRou data obtained between 2019 and 2022 for three well-studied, active M dwarfs: EV Lac, DS Leo, and CN Leo. We looked for secular changes in time series of longitudinal magnetic field, width of unpolarised mean-line profiles, and large-scale field topology as retrieved with principal component analysis and Zeeman-Doppler imaging.

**Results.** We retrieved pulsating (EV Lac), stable (DS Leo), and sine-like (CN Leo) long-term trends in longitudinal field. The width of near-infrared mean-line profiles exhibits rotational modulation only for DS Leo, whereas in the optical it is evident for both EV Lac and DS Leo. The line width variations are not necessarily correlated to those of the longitudinal field, suggesting complex relations between small- and large-scale field. We also recorded topological changes in the form of a reduced axisymmetry for EV Lac and transition from a toroidal-dominated to poloidal-dominated regime for DS Leo. For CN Leo, the topology remained predominantly poloidal, dipolar, and axisymmetric, with only an oscillation in field strength.

**Conclusions.** Our results show a peculiar evolution of the magnetic field for each M dwarf individually, with DS Leo and EV Lac manifesting more evident variations than CN Leo. These findings confirm that M dwarfs with distinct masses and rotation periods can undergo magnetic long-term variations and suggest an underlying variety of cyclic behaviours of their magnetic fields.

**Key words.** techniques: polarimetric – stars: activity – stars: magnetic field – stars: individual: EV Lac – stars: individual: DS Leo – stars: individual: CN Leo

## 1. Introduction

The magnetic fields of low-mass stars ($M_* < 1.2 M_\odot$) are powered by dynamos (Schrijver & Zwaan 2000), and their study is

paramount to understanding stellar evolution and activity phenomena. The evolution of a star's rotation is linked to magnetic activity as the magnetic field couples with the stellar wind and results in angular momentum loss over the star's lifetime (Skumanich 1972; Vidotto et al. 2014; Finley & Matt 2017; See et al. 2019). Therefore, one can infer the stellar age and rotation period using magnetic activity as a proxy (Noyes et al.

<footnote n="1">¹ Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, IRAP/UMR 5277, 14 Avenue Edouard Belin, 31400 Toulouse, France
e-mail: stefano.bellotti@irap.omp.eu</footnote>

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

<footnote n="3">³ Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS, 34095 Montpellier, France</footnote>

<footnote n="4">⁴ Tartu Observatory, University of Tartu, Observatooriumi 1, Töravere 61602 Tartumaa, Estonia</footnote>

<footnote n="5">⁵ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="6">⁶ Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, 37504-364 Itajubá, MG, Brazil</footnote>

<footnote n="7">⁷ Institut d'Astrophysique de Paris, CNRS, UMR 7095, Sorbonne Université, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="8">⁸ Sub-department of Astrophysics, Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

<footnote n="9">⁹ Universidade Federal de Minas Gerais, Belo Horizonte, MG 31270-901, Brazil</footnote>

<footnote n="10">¹⁰ Université de Montréal, Département de Physique, IREX, Montréal, QC H3C 3J7, Canada</footnote>

<footnote n="11">¹¹ Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="12">¹² Observatoire de Genève, Université de Genève, Chemin Pegasi, 51, 1290 Sauverny, Switzerland</footnote>

<footnote n="13">¹³ Observatoire de Haute Provence, St Michel l'Observatoire, France</footnote>

<footnote>* The reduced data are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/686/A66</footnote>
1984; Lorenzo-Oliveira et al. 2018; Dungee et al. 2022). Moreover, stellar magnetic activity is responsible for inhomogeneities in brightness (spots, faculae, plages) and local velocity fields (suppression of convection) that may prevent the unambiguous detection and characterisation of exoplanets, especially those similar to the Earth (Queloz et al. 2001; Huerta et al. 2008; Meunier 2021). Finally, the stellar magnetism dictates the environment in which exoplanets orbit (e.g., Folsom et al. 2020; Bellotti et al. 2023a), influencing their potential habitability (e.g., Segura et al. 2010; Vidotto et al. 2013; Luger & Barnes 2015; Tilley et al. 2019).

On the low-mass end of the main sequence, M dwarfs are the most common spectral type in the solar neighbourhood (Reid et al. 2004), with masses ranging between 0.08 and 0.57 M☉ (Pecaut & Mamajek 2013). Above 0.35 M☉ (approximately M3.5 type), the internal structure is solar-like, that is, with an inner radiative core and an outer convective envelope separated by the tachocline (Chabrier & Baraffe 1997). The tachocline is an interface region of strong shear, where the magnetic field is thought to be amplified. Below 0.35 M☉, M dwarfs possess fully convective interiors, and the absence of a tachocline challenges dynamo theories relying on a deep-seated interface (Durney et al. 1993; Chabrier & Küker 2006; Browning 2008; Yadav et al. 2015). Overall, M dwarfs represent exquisite laboratories to study dynamo-powered magnetic field generation, under similar and different physical interior conditions to those of the Sun.

Dynamo theories are informed by observations of stellar magnetic fields, which can be measured using two complementary approaches (Reiners 2012; Kochukhov 2021). The modulus of the magnetic field vector is estimated via radiative transfer modelling of Zeeman splitting and magnetic intensification of individual unpolarised spectral lines. Values of field strength between 0.2 and 7 KG have been reported (Shulyak et al. 2017, 2019; Reiners et al. 2022; Cristofari et al. 2023), and they follow the activity-rotation relation exhibiting a (quasi-)saturated and non-saturated regime (Reiners & Basri 2009; Shulyak et al. 2019; Reiners et al. 2022). The geometry of the large-scale field can be inferred by means of tomographic inversion techniques, for example Zeeman-Doppler imaging (ZDI; Semel 1989; Donati & Brown 1997) applied to spectropolarimetric time series of linearly and circularly polarised spectra (Morin et al. 2012).

The application of ZDI to spectropolarimetric data of partly and fully convective low-mass stars have revealed a wide variety of large-scale field geometries (e.g., Donati et al. 2008; Morin et al. 2008, 2010; Fares et al. 2013; Hébrard et al. 2016). Partly convective stars with masses higher than 0.5 M☉ tend to have weak, predominantly toroidal, and non-axisymmetric fields, whereas stars with masses between 0.35 and 0.5 M☉ harbour stronger, poloidal, and axisymmetric fields (Donati et al. 2008; Phan-Bao et al. 2009). For fully convective stars close to the 0.35 M☉ boundary, the large-scale field is strong, mainly poloidal, and axisymmetric, while below M < 0.2 M☉ a dichotomy of topologies co-exist: weak, complex, and non-axisymmetric or strong, simple, and axisymmetric (Morin et al. 2008, 2010). This dichotomy can be explained by either two distinct and independent branches of dynamo known as bistability (Morin et al. 2011; Gastine et al. 2013; Kochukhov & Lavail 2017) or by assuming that fully convective M dwarfs undergo magnetic cycles, and previous ZDI reconstructions captured only a snapshot of a long-term topological variation (Kitchatinov et al. 2014). However, no definitive interpretation has been reached so far.

The Sun is an important benchmark for stellar cycles, because the large-scale dipolar component of the magnetic

field undergoes a polarity reversal in a cyclical fashion every 11 yr, and it is accompanied by an oscillation in the fraction of poloidal-to-toroidal magnetic energy (Sanderson et al. 2003; Charbonneau 2010). More precisely, the poloidal component peaks at cycle minimum, and the toroidal components increase towards cycle maximum. However, our understanding of the cyclic nature of the solar magnetism is still not complete (Charbonneau 2020), and its contextualisation advocates for additional magnetic field observations of distinct stellar types. In this direction, informing dynamo theories requires a long-term spectropolarimetric monitoring of selected M dwarfs, for which we can trace the secular evolution of the large-scale field geometry (Klein et al. 2021, 2022; Bellotti et al. 2023b). There is observational evidence of activity cycles for M dwarfs, from photometric and chromospheric activity monitoring (e.g., Suárez Mascareño et al. 2016, 2018; Fuhrmeister et al. 2023; Mignon et al. 2023), and from radial velocity searches of exoplanets (da Gomes Silva et al. 2012; Wargelin et al. 2017; Lopez-Santiago et al. 2020), but we are only starting to capture long-term behaviour that may resemble solar-like magnetic cycles (Bellotti et al. 2023b; Lehmann et al. 2024). For instance, the recent work of Bellotti et al. (2023b) reports the long-term evolution of the magnetic field of AD Leo, which has similarities with the evolution of the Sun's field.

With this work we investigate the long-term monitoring of three active M dwarfs, EV Lac, DS Leo, and CN Leo, by analysing spectropolarimetric data collected across the optical and near-infrared domain. To investigate the long-term evolution, we analysed the full time series of the longitudinal magnetic field ($B_t$) and the full width at half maximum (FWHM) of unpolarised (Stokes I) mean-line profiles. We reconstructed the large-scale field topology via ZDI, and we inspected temporal changes in the morphology of circular polarisation profiles via principal component analysis (PCA; Lehmann & Donati 2022). The paper is structured as follows. In Sect. 2 we describe the near-infrared and optical spectropolarimetric observations. We outline the analysis of the longitudinal magnetic field in Sect. 3, the analysis FWHM of Stokes I in Sect. 4, the PCA technique in Sect. 5, and the ZDI reconstructions in Sect. 6. Finally, we present our conclusions in Sect. 7.

## 2. Observations

We used spectropolarimetric data collected in near-infrared and optical domains, and performed a large-scale magnetic field monitoring for three M-type stars: EV Lac (GJ 873), DS Leo (GJ 410), and CN Leo (GJ 406). A summary of their properties is given in Table 1.

According to stellar models, the transition between partly and fully convective interiors spans from 0.20 to 0.35 M☉ (Dorman et al. 1989; Chabrier & Baraffe 1997; Rabus et al. 2019). Other factors like age, metallicity and magnetic field strength play a role (Mullan & MacDonald 2001; Maeder & Meynet 2000; van Saders & Pinsonneault 2012; Tanner et al. 2013). With such a range, DS Leo and CN Leo are partly and fully convective, respectively. With a mass of 0.32 M☉, EV Lac is either partly convective or near the transition between the two regimes.

### 2.1. Near-infrared

All near-infrared observations were performed in circular polarisation mode with the SpectroPolarimètre InfraRouge (SPIRou)
1984; Lorenzo-Oliveira et al. 2018; Dungee et al. 2022). Moreover, stellar magnetic activity is responsible for inhomogeneities in brightness (spots, faculae, plages) and local velocity fields (suppression of convection) that may prevent the unambiguous detection and characterisation of exoplanets, especially those similar to the Earth (Queloz et al. 2001; Huerta et al. 2008; Meunier 2021). Finally, the stellar magnetism dictates the environment in which exoplanets orbit (e.g., Folsom et al. 2020; Bellotti et al. 2023a), influencing their potential habitability (e.g., Segura et al. 2010; Vidotto et al. 2013; Luger & Barnes 2015; Tilley et al. 2019).

On the low-mass end of the main sequence, M dwarfs are the most common spectral type in the solar neighbourhood (Reid et al. 2004), with masses ranging between 0.08 and 0.57 M☉ (Pecaut & Mamajek 2013). Above 0.35 M☉ (approximately M3.5 type), the internal structure is solar-like, that is, with an inner radiative core and an outer convective envelope separated by the tachocline (Chabrier & Baraffe 1997). The tachocline is an interface region of strong shear, where the magnetic field is thought to be amplified. Below 0.35 M☉, M dwarfs possess fully convective interiors, and the absence of a tachocline challenges dynamo theories relying on a deep-seated interface (Durney et al. 1993; Chabrier & Küker 2006; Browning 2008; Yadav et al. 2015). Overall, M dwarfs represent exquisite laboratories to study dynamo-powered magnetic field generation, under similar and different physical interior conditions to those of the Sun.

Dynamo theories are informed by observations of stellar magnetic fields, which can be measured using two complementary approaches (Reiners 2012; Kochukhov 2021). The modulus of the magnetic field vector is estimated via radiative transfer modelling of Zeeman splitting and magnetic intensification of individual unpolarised spectral lines. Values of field strength between 0.2 and 7 KG have been reported (Shulyak et al. 2017, 2019; Reiners et al. 2022; Cristofari et al. 2023), and they follow the activity-rotation relation exhibiting a (quasi-)saturated and non-saturated regime (Reiners & Basri 2009; Shulyak et al. 2019; Reiners et al. 2022). The geometry of the large-scale field can be inferred by means of tomographic inversion techniques, for example Zeeman-Doppler imaging (ZDI; Semel 1989; Donati & Brown 1997) applied to spectropolarimetric time series of linearly and circularly polarised spectra (Morin et al. 2012).

The application of ZDI to spectropolarimetric data of partly and fully convective low-mass stars have revealed a wide variety of large-scale field geometries (e.g., Donati et al. 2008; Morin et al. 2008, 2010; Fares et al. 2013; Hébrard et al. 2016). Partly convective stars with masses higher than 0.5 M☉ tend to have weak, predominantly toroidal, and non-axisymmetric fields, whereas stars with masses between 0.35 and 0.5 M☉ harbour stronger, poloidal, and axisymmetric fields (Donati et al. 2008; Phan-Bao et al. 2009). For fully convective stars close to the 0.35 M☉ boundary, the large-scale field is strong, mainly poloidal, and axisymmetric, while below M < 0.2 M☉ a dichotomy of topologies co-exist: weak, complex, and non-axisymmetric or strong, simple, and axisymmetric (Morin et al. 2008, 2010). This dichotomy can be explained by either two distinct and independent branches of dynamo known as bistability (Morin et al. 2011; Gastine et al. 2013; Kochukhov & Lavail 2017) or by assuming that fully convective M dwarfs undergo magnetic cycles, and previous ZDI reconstructions captured only a snapshot of a long-term topological variation (Kitchatinov et al. 2014). However, no definitive interpretation has been reached so far.

The Sun is an important benchmark for stellar cycles, because the large-scale dipolar component of the magnetic

field undergoes a polarity reversal in a cyclical fashion every 11 yr, and it is accompanied by an oscillation in the fraction of poloidal-to-toroidal magnetic energy (Sanderson et al. 2003; Charbonneau 2010). More precisely, the poloidal component peaks at cycle minimum, and the toroidal components increase towards cycle maximum. However, our understanding of the cyclic nature of the solar magnetism is still not complete (Charbonneau 2020), and its contextualisation advocates for additional magnetic field observations of distinct stellar types. In this direction, informing dynamo theories requires a long-term spectropolarimetric monitoring of selected M dwarfs, for which we can trace the secular evolution of the large-scale field geometry (Klein et al. 2021, 2022; Bellotti et al. 2023b). There is observational evidence of activity cycles for M dwarfs, from photometric and chromospheric activity monitoring (e.g., Suárez Mascareño et al. 2016, 2018; Fuhrmeister et al. 2023; Mignon et al. 2023), and from radial velocity searches of exoplanets (da Gomes Silva et al. 2012; Wargelin et al. 2017; Lopez-Santiago et al. 2020), but we are only starting to capture long-term behaviour that may resemble solar-like magnetic cycles (Bellotti et al. 2023b; Lehmann et al. 2024). For instance, the recent work of Bellotti et al. (2023b) reports the long-term evolution of the magnetic field of AD Leo, which has similarities with the evolution of the Sun's field.

With this work we investigate the long-term monitoring of three active M dwarfs, EV Lac, DS Leo, and CN Leo, by analysing spectropolarimetric data collected across the optical and near-infrared domain. To investigate the long-term evolution, we analysed the full time series of the longitudinal magnetic field ($B_t$) and the full width at half maximum (FWHM) of unpolarised (Stokes I) mean-line profiles. We reconstructed the large-scale field topology via ZDI, and we inspected temporal changes in the morphology of circular polarisation profiles via principal component analysis (PCA; Lehmann & Donati 2022). The paper is structured as follows. In Sect. 2 we describe the near-infrared and optical spectropolarimetric observations. We outline the analysis of the longitudinal magnetic field in Sect. 3, the analysis FWHM of Stokes I in Sect. 4, the PCA technique in Sect. 5, and the ZDI reconstructions in Sect. 6. Finally, we present our conclusions in Sect. 7.

## 2. Observations

We used spectropolarimetric data collected in near-infrared and optical domains, and performed a large-scale magnetic field monitoring for three M-type stars: EV Lac (GJ 873), DS Leo (GJ 410), and CN Leo (GJ 406). A summary of their properties is given in Table 1.

According to stellar models, the transition between partly and fully convective interiors spans from 0.20 to 0.35 M☉ (Dorman et al. 1989; Chabrier & Baraffe 1997; Rabus et al. 2019). Other factors like age, metallicity and magnetic field strength play a role (Mullan & MacDonald 2001; Maeder & Meynet 2000; van Saders & Pinsonneault 2012; Tanner et al. 2013). With such a range, DS Leo and CN Leo are partly and fully convective, respectively. With a mass of 0.32 M☉, EV Lac is either partly convective or near the transition between the two regimes.

### 2.1. Near-infrared

All near-infrared observations were performed in circular polarisation mode with the SpectroPolarimètre InfraRouge (SPIRou)
<table><caption>**Table 1.** Properties of the M dwarfs examined.</caption><thead><tr><th>Star</th><th>ID</th><th>Sp type</th><th><i>M</i><sub>*</sub><br>[<i>M</i><sub>☉</sub>]</th><th><i>R</i><sub>*</sub><br>[<i>R</i><sub>☉</sub>]</th><th><i>P</i><sub>rot</sub><br>[d]</th><th><i>v</i><sub>eq</sub> sin <i>i</i><br>[km s<sup>-1</sup>]</th><th><i>i</i><br>[deg]</th><th>log(<i>L</i><sub>X</sub>/<i>L</i><sub>bol</sub>)<br>[dex]</th><th>log<i>R'</i><sub>HK</sub><br>[dex]</th></tr></thead><tbody><tr><td>DS Leo</td><td>GJ 410</td><td>M1.0</td><td>0.58</td><td>0.53</td><td>13.91 ± 0.01<sup>(†)</sup></td><td>2.0</td><td>60</td><td>-3.80</td><td>-4.16</td></tr><tr><td>EV Lac</td><td>GJ 873</td><td>M3.5</td><td>0.32</td><td>0.30</td><td>4.36 ± 0.01<sup>(†)</sup></td><td>4.0</td><td>60</td><td>-1.99</td><td>-3.75</td></tr><tr><td>CN Leo</td><td>GJ 406</td><td>M5.5</td><td>0.10</td><td>0.12</td><td>2.70 ± 0.01<sup>(†)</sup></td><td>2.0</td><td>45 ± 20<sup>(†)</sup></td><td>...</td><td>-4.01</td></tr></tbody></table>

**Notes.** The columns list the following quantities: 1) and 2) stellar name and alternative identifier, 3) spectral type, 4) mass, 5) radius, 6) rotation period at equator, 7) equatorial projected rotational velocity, 8) inclination, 9) X-ray-to-bolometric luminosity ratio (Wright et al. 2011), and 10) CaII H&K index (Noyes et al. 1984; Boro Saikia et al. 2018). Columns 4) and 5) are from Cristofari et al. (2023), and columns 7) and 8) are from Morin et al. (2008), Hébrard et al. (2016), and Cristofari et al. (2023). Columns 9) and 10) are from Wright et al. (2011), Stelzer et al. (2013) and Boro Saikia et al. (2018), respectively. The dagger (†) indicates that the parameter was estimated in this work.

as part of the large programme SPIRou Legacy Survey<fnref n="1" /> (SLS; id P42, PI: Jean-François Donati). SPIRou is a stabilised high-resolution near-infrared spectropolarimeter (Donati et al. 2020) mounted on the 3.6 m Canada-France-Hawaii Telescope (CFHT) atop Maunakea, Hawaii. It provides a quasi-continuous coverage of the near-infrared spectrum from 0.96 to 2.5 µm (YJHK bands) at a spectral resolving power of R ~ 70 000, with a 2-nm gap between 2.4371 and 2.4391 µm (Donati et al. 2020). Optimal extraction of SPIRou spectra was carried out with A PipelinE to Reduce Observations (APER0 v0.6.132), a fully automatic reduction package installed at CFHT (Cook et al. 2022). The journal of the observations is given in Appendix E.

Starting from the polarimetric products of APERO (see Cook et al. 2022, for more details), we computed Stokes I (unpolarised) and V (circularly polarised) mean profiles using least-squares deconvolution (LSD; Donati et al. 1997; Kochukhov et al. 2010).<fnref n="2" /> With this technique the observed spectrum is deconvolved with a line list, namely a series of Dirac delta functions located at each absorption line in the stellar spectrum and with the associated line features such as depth, and sensitivity to Zeeman effect (commonly known as Landé factor and indicated as $g_{eff}$). The deconvolution results in an individual, high-signal-to-noise ratio (S/N) kernel summarising the properties of thousands of spectral lines, and allowing us to extract polarimetric information from the spectrum.

For our stars, we used two line lists corresponding to a local thermodynamic equilibrium model (Gustafsson et al. 2008) characterised by $log g = 5.0$ [cm s⁻<fnref n="2" />], $v_{micro} = 1$ km s⁻<fnref n="1" />, and $T_{eff} = 3000$ K (for CN Leo) and $T_{eff} = 3500$ K (for EV Lac and DS Leo). These two lists contain 1000 and 1400 atomic photospheric lines between 950 and 2600 nm with depth larger than 3% the continuum level. The depth threshold is chosen to remove shallow lines (with low effective S/N), while keeping a large number of lines with which to compute LSD profiles. The masks were synthesised using the Vienna Atomic Line Database<fnref n="3" /> (VALD, Ryabchikova et al. 2015), and contain information on $g_{eff}$ (ranging from 0 to 3), which is not accessible with an empirical or molecular mask. Given the spectral type of EV Lac and DS Leo, we respectively tested a mask associated with 3000 K (590 lines) 4000 K (1550 lines), to check the robustness of our results. For both Stokes I and V, we adopted the weighting factor $dg_{eff}λ/(dn_g g_{eff,n}λ_n)$, with d is the line depth

and λ is the wavelength, whereas the quantities at denominator are the normalisation parameters (see Kochukhov et al. 2010, for more details).

A forest of telluric absorption lines due to the Earth's atmosphere pollute the wavelength domain in which SPIRou operates. As described in Cook et al. (2022) and Artigau et al. (2022), the removal of telluric contamination from science data is performed by APERO with a two-step algorithm. First, a pre-cleaning is carried out by means of a TAPAS (Transmissions of the Atmosphere for Stromical data; Bertaux et al. 2014) absorption model, which is also applied to a set of hot stars to obtain a library of telluric-correction residuals. From this, a telluric-residuals model is constructed and subtracted from the science frames. To account for potential residuals in the telluric correction, we ignored the following intervals of the spectrum when computing the LSD profiles: [950, 979], [1081, 1169], [1328, 1492], [1784, 2029], [2380, 2500] nm. These correspond to spectral windows where the telluric absorption is highest. The final number of spectral lines used in LSD is 588 and 830 for the $T_{eff} = 3000$ K and $T_{eff} = 3500$ K line lists, respectively.

We show an example of LSD Stokes profiles for the examined stars in Fig. 1. As mentioned by Lavail et al. (2018), the continuum of the Stokes I LSD profiles for M dwarfs lies below unity owing to molecular spectral lines that are not accounted in the line lists. We thus re-normalised the profiles to unity by fitting a linear model to the region outside the line, to include residuals of continuum normalisation at the level of the spectra. The Stokes V and N profiles were correspondingly re-scaled with the same fit. The operation of continuum normalisation does not alter the magnetic analyses presented here appreciably. For instance, the longitudinal field values computed with non-normalised LSD profiles are at most 5% larger compared to the values obtained from normalised LSD profiles. In a similar manner, the retrieved field topology and evolution are robust against LSD profiles normalisation, as long as a multi-epoch time series for a particular target is analysed in a consistent way.

EV Lac was observed for 163 nights between September 2019 and November 2021, spanning a period of ~800 days. We split the entire time series in smaller epochs based on the instrumental scheduling and visibility gaps, following Bellotti et al. (2023b). We obtained three epochs: 2019b (35 obs between 2019.86 and 2019.95), 2020b2021a (67 observations between 2020.57 and 2021.02) and 2021b (61 observations between 2021.47 and 2021.72). We recorded a S/N at 1650 nm per spectral element between 31 and 165, with a mean of 130, and an average airmass of 1.2.

DS Leo was observed for 130 nights between November 2020 and June 2022, spanning a period of ~580 days. The time series was split in two epochs: 2020b2021a (61 observations

<footnote n="1">¹ http://spirou.irap.omp.eu/Observations/The-SPIRou-Legacy-Survey</footnote>

<footnote n="2">² We used the python LSD code available at https://github.com/folsomcp/LSDpy</footnote>

<footnote n="3">³ http://vald.astro.uu.se/ using the Montpellier mirror to request locally MARCS model atmospheres.</footnote>
<table><caption>**Table 1.** Properties of the M dwarfs examined.</caption><thead><tr><th>Star</th><th>ID</th><th>Sp type</th><th>M<sub>*</sub><br>[M<sub>☉</sub>]</th><th>R<sub>*</sub><br>[R<sub>☉</sub>]</th><th>P<sub>rot</sub><br>[d]</th><th>v<sub>eq</sub> sin i<br>[km s<sup>-1</sup>]</th><th>i<br>[deg]</th><th>log(L<sub>X</sub>/L<sub>bol</sub>)<br>[dex]</th><th>log R'<sub>HK</sub><br>[dex]</th></tr></thead><tbody><tr><td>DS Leo</td><td>GJ 410</td><td>M1.0</td><td>0.58</td><td>0.53</td><td>13.91 ± 0.01<sup>(†)</sup></td><td>2.0</td><td>60</td><td>-3.80</td><td>-4.16</td></tr><tr><td>EV Lac</td><td>GJ 873</td><td>M3.5</td><td>0.32</td><td>0.30</td><td>4.36 ± 0.01<sup>(†)</sup></td><td>4.0</td><td>60</td><td>-1.99</td><td>-3.75</td></tr><tr><td>CN Leo</td><td>GJ 406</td><td>M5.5</td><td>0.10</td><td>0.12</td><td>2.70 ± 0.01<sup>(†)</sup></td><td>2.0</td><td>45 ± 20<sup>(†)</sup></td><td>...</td><td>-4.01</td></tr></tbody></table>

**Notes.** The columns list the following quantities: 1) and 2) stellar name and alternative identifier, 3) spectral type, 4) mass, 5) radius, 6) rotation period at equator, 7) equatorial projected rotational velocity, 8) inclination, 9) X-ray-to-bolometric luminosity ratio (Wright et al. 2011), and 10) CaII H&K index (Noyes et al. 1984; Boro Saikia et al. 2018). Columns 4) and 5) are from Cristofari et al. (2023), and columns 7) and 8) are from Morin et al. (2008), Hébrard et al. (2016), and Cristofari et al. (2023). Columns 9) and 10) are from Wright et al. (2011), Stelzer et al. (2013) and Boro Saikia et al. (2018), respectively. The dagger (†) indicates that the parameter was estimated in this work.

as part of the large programme SPIRou Legacy Survey<fnref n="1" /> (SLS; id P42, PI: Jean-François Donati). SPIRou is a stabilised high-resolution near-infrared spectropolarimeter (Donati et al. 2020) mounted on the 3.6 m Canada-France-Hawaii Telescope (CFHT) atop Maunakea, Hawaii. It provides a quasi-continuous coverage of the near-infrared spectrum from 0.96 to 2.5 µm (YJHK bands) at a spectral resolving power of R ~ 70 000, with a 2-nm gap between 2.4371 and 2.4391 µm (Donati et al. 2020). Optimal extraction of SPIRou spectra was carried out with A PipelinE to Reduce Observations (APER0 v0.6.132), a fully automatic reduction package installed at CFHT (Cook et al. 2022). The journal of the observations is given in Appendix E.

Starting from the polarimetric products of APERO (see Cook et al. 2022, for more details), we computed Stokes I (unpolarised) and V (circularly polarised) mean profiles using least-squares deconvolution (LSD; Donati et al. 1997; Kochukhov et al. 2010).<fnref n="2" /> With this technique the observed spectrum is deconvolved with a line list, namely a series of Dirac delta functions located at each absorption line in the stellar spectrum and with the associated line features such as depth, and sensitivity to Zeeman effect (commonly known as Landé factor and indicated as $g_{eff}$). The deconvolution results in an individual, high-signal-to-noise ratio (S/N) kernel summarising the properties of thousands of spectral lines, and allowing us to extract polarimetric information from the spectrum.

For our stars, we used two line lists corresponding to a local thermodynamic equilibrium model (Gustafsson et al. 2008) characterised by $log g = 5.0$ [cm s$^{-2}$], $v_{micro} = 1$ km s$^{-1}$, and $T_{eff} = 3000$ K (for CN Leo) and $T_{eff} = 3500$ K (for EV Lac and DS Leo). These two lists contain 1000 and 1400 atomic photospheric lines between 950 and 2600 nm with depth larger than 3% the continuum level. The depth threshold is chosen to remove shallow lines (with low effective S/N), while keeping a large number of lines with which to compute LSD profiles. The masks were synthesised using the Vienna Atomic Line Database<fnref n="3" /> (VALD, Ryabchikova et al. 2015), and contain information on $g_{eff}$ (ranging from 0 to 3), which is not accessible with an empirical or molecular mask. Given the spectral type of EV Lac and DS Leo, we respectively tested a mask associated with 3000 K (590 lines) 4000 K (1550 lines), to check the robustness of our results. For both Stokes I and V, we adopted the weighting factor $dg_{eff}λ/(dn_g g_{eff,n}λ_n)$, with d is the line depth

and λ is the wavelength, whereas the quantities at denominator are the normalisation parameters (see Kochukhov et al. 2010, for more details).

A forest of telluric absorption lines due to the Earth's atmosphere pollute the wavelength domain in which SPIRou operates. As described in Cook et al. (2022) and Artigau et al. (2022), the removal of telluric contamination from science data is performed by APERO with a two-step algorithm. First, a pre-cleaning is carried out by means of a TAPAS (Transmissions of the Atmosphere for Stromical data; Bertaux et al. 2014) absorption model, which is also applied to a set of hot stars to obtain a library of telluric-correction residuals. From this, a telluric-residuals model is constructed and subtracted from the science frames. To account for potential residuals in the telluric correction, we ignored the following intervals of the spectrum when computing the LSD profiles: [950, 979], [1081, 1169], [1328, 1492], [1784, 2029], [2380, 2500] nm. These correspond to spectral windows where the telluric absorption is highest. The final number of spectral lines used in LSD is 588 and 830 for the $T_{eff} = 3000$ K and $T_{eff} = 3500$ K line lists, respectively.

We show an example of LSD Stokes profiles for the examined stars in Fig. 1. As mentioned by Lavail et al. (2018), the continuum of the Stokes I LSD profiles for M dwarfs lies below unity owing to molecular spectral lines that are not accounted in the line lists. We thus re-normalised the profiles to unity by fitting a linear model to the region outside the line, to include residuals of continuum normalisation at the level of the spectra. The Stokes V and N profiles were correspondingly re-scaled with the same fit. The operation of continuum normalisation does not alter the magnetic analyses presented here appreciably. For instance, the longitudinal field values computed with non-normalised LSD profiles are at most 5% larger compared to the values obtained from normalised LSD profiles. In a similar manner, the retrieved field topology and evolution are robust against LSD profiles normalisation, as long as a multi-epoch time series for a particular target is analysed in a consistent way.

EV Lac was observed for 163 nights between September 2019 and November 2021, spanning a period of ~800 days. We split the entire time series in smaller epochs based on the instrumental scheduling and visibility gaps, following Bellotti et al. (2023b). We obtained three epochs: 2019b (35 obs between 2019.86 and 2019.95), 2020b2021a (67 observations between 2020.57 and 2021.02) and 2021b (61 observations between 2021.47 and 2021.72). We recorded a S/N at 1650 nm per spectral element between 31 and 165, with a mean of 130, and an average airmass of 1.2.

DS Leo was observed for 130 nights between November 2020 and June 2022, spanning a period of ~580 days. The time series was split in two epochs: 2020b2021a (61 observations

<footnote n="1">¹ http://spirou.irap.omp.eu/Observations/The-SPIRou-Legacy-Survey</footnote>

<footnote n="2">² We used the python LSD code available at https://github.com/folsomcp/LSDpy</footnote>

<footnote n="3">³ http://vald.astro.uu.se/ using the Montpellier mirror to request locally MARCS model atmospheres.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Examples of Stokes profiles obtained with SPIRou for EV Lac, DS Leo and CN Leo. In each panel, Stokes I (top), V (middle), and null polarisation profile (bottom; see Bagnulo et al. 2009) are illustrated in units of unpolarised continuum ($I_c$). Clear circular polarisation signatures are present, with a S/N of 5000, 5600, and 3080 for EV Lac, DS Leo, and CN Leo, respectively. The different magnetic activity level of the stars is shown by the different widths of the average line profiles in both unpolarised and polarised light.</figcaption></figure>

between 2020.27 and 2021.55) and 2021b2022a (69 observations between 2021.71 and 2022.46). The S/N ranged between 51 and 151, with a mean of 125, and the mean airmass is 1.3.

CN Leo was observed for 169 nights between April 2019 and June 2022, spanning a period of ~1150 days. The time series was split in four epochs: 2019a (19 observations between 2019.29 and 2019.47), 2019b2020a (37 observations between 2019.83 and 2020.44) and 2020b2021a (63 observations between 2020.83 and 2021.50) 2021b2022a (45 observations between 2021.87 and 2022.44). The S/N ranged between 40 and 146, with a mean of 111, and the average airmass is 1.3. Six observations in February 2019 were excluded due to an opti-

cal component not working nominally at very low temperatures, in a similar way to Bellotti et al. (2023b).

## 2.2. Optical

We used archival optical spectropolarimetric observations collected with ESPaDOnS and Narval. ESPaDOnS is the spectropolarimeter on the 3.6 m CFHT located atop Mauna Kea in Hawaii (Donati et al. 2006a), and Narval is the twin instrument on the 2 m Télescope Bernard Lyot (TBL) at the Pic du Midi Observatory in France (Donati 2003). Data reduction was performed with LIBRE-ESPRIT (Donati et al. 1997), and the continuum normalised spectra were retrieved from PolarBase (Petit et al. 2014).

The optical data sets were already examined in previous studies (Donati et al. 2008; Morin et al. 2008, 2010; Hébrard et al. 2016). For EV Lac, there are 79 observations taken between 2005 and 2016. The 2010 observations were affected by a malfunction of the rhombs which affected the circular polarisation signal, hence only the Stokes I data are used in the following analyses. For DS Leo, there are 94 observations between 2006 and 2014. For CN Leo, there are four observations in 2008, which we included only for the FWHM analysis (see Sect. 4).

Similarly to the near-infrared, we computed Stokes I and V profiles with LSD. For EV Lac and DS Leo, we adopted an optical 3500 K VALD mask containing 3240 lines in range 350–1080 nm and with depths larger than 40% the continuum level, similarly to Morin et al. (2008), Bellotti et al. (2022). For CN Leo, we used a 3000 K mask built in a similar manner and containing 3492 absorption lines. The number of lines in both masks already accounts for the removal of the wavelength intervals affected by telluric lines analogously to Bellotti et al. (2022). These are: [627,632], [655,5,657], [686,697], [716,734], [759,770], [813,835], and [895,986] nm. In the next sections, the observations will be phased with the following ephemeris

$$ \text{HJD}_{\text{EVLac}} = 2458738.0805 + \text{P}_{\text{rot,EVLac}} \cdot n_{\text{cyc}} \quad (1) $$

$$ \text{HJD}_{\text{DSLeo}} = 2459158.1132 + \text{P}_{\text{rot,DSLeo}} \cdot n_{\text{cyc}} \quad (2) $$

$$ \text{HJD}_{\text{CNLeo}} = 2458590.0095 + \text{P}_{\text{rot,CNLeo}} \cdot n_{\text{cyc}} \quad (3) $$

where we separately used the first SPIRou observation for each star as heliocentric Julian Date reference, $P_{\text{rot,i}}$ is the stellar rotation period of the ith star (see Sect. 3), and $n_{\text{cyc}}$ represents the rotation cycle (tabulated in Appendix E).

## 3. Longitudinal magnetic field

The disk-integrated, line-of-sight component of the magnetic field ($B_l$) is computed from the first-order moment of a Stokes V profile (Rees & Semel 1979; Folsom et al. 2016). Formally,

$$ B_l [\text{G}] = \frac{-2.14 \times 10^{11} \int vV(v)dv}{\lambda_0 g_{\text{eff}} c \int (I_c - I)dv}, \quad (4) $$

where $\lambda_0$ and $g_{\text{eff}}$ are the normalisation wavelength (in nm) and Landé factor of the LSD profiles, $I_c$ is the continuum level, $v$ is the radial velocity associated to a point in the spectral line profile in the star's rest frame (in km s⁻¹) and c the speed of light in vacuum (in km s⁻¹).

For all stars examined, the adopted normalisation wavelength is 1700 nm and 700 nm for the near-infrared and optical domains, respectively. The normalisation Landé factor is 1.191
<figure><img src="image_2.png" /><figcaption>Fig. 1. Examples of Stokes profiles obtained with SPIRou for EV Lac, DS Leo and CN Leo. In each panel, Stokes I (top), V (middle), and null polarisation profile (bottom; see Bagnulo et al. 2009) are illustrated in units of unpolarised continuum ($I_c$). Clear circular polarisation signatures are present, with a S/N of 5000, 5600, and 3080 for EV Lac, DS Leo, and CN Leo, respectively. The different magnetic activity level of the stars is shown by the different widths of the average line profiles in both unpolarised and polarised light.</figcaption></figure>

between 2020.27 and 2021.55) and 2021b2022a (69 observations between 2021.71 and 2022.46). The S/N ranged between 51 and 151, with a mean of 125, and the mean airmass is 1.3.

CN Leo was observed for 169 nights between April 2019 and June 2022, spanning a period of ~1150 days. The time series was split in four epochs: 2019a (19 observations between 2019.29 and 2019.47), 2019b2020a (37 observations between 2019.83 and 2020.44) and 2020b2021a (63 observations between 2020.83 and 2021.50) 2021b2022a (45 observations between 2021.87 and 2022.44). The S/N ranged between 40 and 146, with a mean of 111, and the average airmass is 1.3. Six observations in February 2019 were excluded due to an opti-

cal component not working nominally at very low temperatures, in a similar way to Bellotti et al. (2023b).

## 2.2. Optical

We used archival optical spectropolarimetric observations collected with ESPaDOnS and Narval. ESPaDOnS is the spectropolarimeter on the 3.6 m CFHT located atop Mauna Kea in Hawaii (Donati et al. 2006a), and Narval is the twin instrument on the 2 m Télescope Bernard Lyot (TBL) at the Pic du Midi Observatory in France (Donati 2003). Data reduction was performed with LIBRE-ESPRIT (Donati et al. 1997), and the continuum normalised spectra were retrieved from PolarBase (Petit et al. 2014).

The optical data sets were already examined in previous studies (Donati et al. 2008; Morin et al. 2008, 2010; Hébrard et al. 2016). For EV Lac, there are 79 observations taken between 2005 and 2016. The 2010 observations were affected by a malfunction of the rhombs which affected the circular polarisation signal, hence only the Stokes I data are used in the following analyses. For DS Leo, there are 94 observations between 2006 and 2014. For CN Leo, there are four observations in 2008, which we included only for the FWHM analysis (see Sect. 4).

Similarly to the near-infrared, we computed Stokes I and V profiles with LSD. For EV Lac and DS Leo, we adopted an optical 3500 K VALD mask containing 3240 lines in range 350–1080 nm and with depths larger than 40% the continuum level, similarly to Morin et al. (2008), Bellotti et al. (2022). For CN Leo, we used a 3000 K mask built in a similar manner and containing 3492 absorption lines. The number of lines in both masks already accounts for the removal of the wavelength intervals affected by telluric lines analogously to Bellotti et al. (2022). These are: [627,632], [655,5,657], [686,697], [716,734], [759,770], [813,835], and [895,986] nm. In the next sections, the observations will be phased with the following ephemeris

$$ \text{HJD}_{\text{EVLac}} = 2458738.0805 + \text{P}_{\text{rot,EVLac}} \cdot n_{\text{cyc}} \quad (1) $$

$$ \text{HJD}_{\text{DSLeo}} = 2459158.1132 + \text{P}_{\text{rot,DSLeo}} \cdot n_{\text{cyc}} \quad (2) $$

$$ \text{HJD}_{\text{CNLeo}} = 2458590.0095 + \text{P}_{\text{rot,CNLeo}} \cdot n_{\text{cyc}} \quad (3) $$

where we separately used the first SPIRou observation for each star as heliocentric Julian Date reference, $P_{\text{rot,i}}$ is the stellar rotation period of the ith star (see Sect. 3), and $n_{\text{cyc}}$ represents the rotation cycle (tabulated in Appendix E).

## 3. Longitudinal magnetic field

The disk-integrated, line-of-sight component of the magnetic field ($B_l$) is computed from the first-order moment of a Stokes V profile (Rees & Semel 1979; Folsom et al. 2016). Formally,

$$ B_l [\text{G}] = \frac{-2.14 \times 10^{11} \int vV(v)dv}{\lambda_0 g_{\text{eff}} c \int (I_c - I)dv}, \quad (4) $$

where $\lambda_0$ and $g_{\text{eff}}$ are the normalisation wavelength (in nm) and Landé factor of the LSD profiles, $I_c$ is the continuum level, $v$ is the radial velocity associated to a point in the spectral line profile in the star's rest frame (in km s⁻¹) and c the speed of light in vacuum (in km s⁻¹).

For all stars examined, the adopted normalisation wavelength is 1700 nm and 700 nm for the near-infrared and optical domains, respectively. The normalisation Landé factor is 1.191
in near-infrared and 1.154 in optical for EV Lac, 1.181 and 1.086 for DS Leo, and 1.223 for CN Leo (only near-infrared). Considering the larger line width in near-infrared than optical, the integration is computed within ±50 km s⁻¹ and ±30 km s⁻¹ from line centre for EV Lac and within ±60 km s⁻¹ and ±30 km s⁻¹ for CN Leo. For DS Leo, the integration intervals are ±25 km s⁻¹ and ±20 km s⁻¹, since it is a slower rotator.

The longitudinal field is a useful diagnostic of surface magnetic features, whose appearance on the visible disk is modulated at the stellar rotation period. For this reason, practical information can be extracted from the time series of $B_l$ values with standard techniques (Folsom et al. 2018; Petit et al. 2021; Carmona et al. 2023; Fouqué et al. 2023; Donati et al. 2023b). We applied a Generalised Lomb-Scargle periodogram (Zechmeister & Kürster 2009) to the entire time series as well as the individual subsets of data for all the examined stars. The results are illustrated in Fig. 2.

For EV Lac and DS Leo, the peak at the expected rotation period is unambiguously higher than a false alarm probability (FAP) of 0.01%, which emphasises the dominant stellar activity signal and confirms $B_l$ as suitable diagnostic. For CN Leo, the interpretation of the periodogram is more elaborate because of strong aliases of the observational cadence (~1 d), the gap between instrument runs (~30 d), and the time span of the time series (~1000 d). By ignoring these features, the periodogram shows a clear and unique peak at the expected stellar rotation period, for which there is no correspondence with the observing window function. If we restrict the periodogram analysis to each of the four subsets of CN Leo (see Sect. 2), only the last one 2021b2022a does not show a spike at the expected rotation period. For all stars, we retrieve $P_{rot}$ in agreement to known values from the literature: $P_{rot,EVLac}$ = 4.36 ± 0.01 d (Morin et al. 2008), $P_{rot,DSLeo}$ = 13.91 ± 0.01 d (Hébrard et al. 2016), and $P_{rot,CNLeo}$ = 2.70 ± 0.01 d (Reinhold & Hekker 2020; Lafarga et al. 2021; Irving et al. 2023).

The longitudinal field computation was carried out for both optical and near-infrared domains for all stars, and is shown in Figs. 3–5. The journal of $B_l$ values is reported in Appendix E.

## 3.1. EV Lac

EV Lac's near-infrared measurements lie between 270 and -470 G, with a mean of -50 G and a median error bar of 19 G. In comparison, the optical data have a larger spread, ranging between 380 and -540 G with a mean of -90 G and a median error bar of 16 G. When a 3000 K, near-infrared mask is used to compute the LSD profiles, the range of $B_l$ spans between -480 and 286 G, with a mean of -52 G and a median error bar of 21 G. The Stokes *I* profiles derived with the colder mask are on average 20% deeper than those associated with the 3500 K mask, and the amplitude of Stokes *V* profiles is on average ~20% larger. These differences in the Stokes profiles cancel out according to Eq. (4), ultimately making the results consistent. Overall, the near-infrared $B_l$ values seem to progressively restore the amplitude with time, after a potential minimum in 2016, making the overall trend look like a beating signal. However, the dearth of data between 2010 and 2019 prevents us from constraining the phenomenon further.

When EV Lac's near-infrared data are phase-folded at $P_{rot,EVLac}$ = 4.36 d, we note a clear rotational modulation for the three epochs (see Fig. 3). The phase variations are modelled with a Levenberg-Marquardt least squares sine fit, and we observe an increasing amplitude of the field variations, indicating a decrease in the field axisymmetry or a field enhancement

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Generalised Lomb-Scargle periodogram of the longitudinal field near-infrared time series, from top to bottom EV Lac, DS Leo, and CN Leo. In each panel the rotation period corresponding to the highest peak is highlighted with a green dotted line, while two FAP levels (0.1% and 0.01%) are shown as grey horizontal lines. The window function of the entire time series is included and mirrored with respect to the x-axis (VanderPlas 2018) to highlight aliases due to the observation cadence and scheduling.</figcaption></figure>

over time. A detailed investigation will be performed with Zeeman-Doppler imaging in Sect. 6.

## 3.2. DS Leo

DS Leo is the least magnetically active star of our sample. The longitudinal field near-infrared measurements fall within ±56 G
in near-infrared and 1.154 in optical for EV Lac, 1.181 and 1.086 for DS Leo, and 1.223 for CN Leo (only near-infrared). Considering the larger line width in near-infrared than optical, the integration is computed within ±50 km s⁻¹ and ±30 km s⁻¹ from line centre for EV Lac and within ±60 km s⁻¹ and ±30 km s⁻¹ for CN Leo. For DS Leo, the integration intervals are ±25 km s⁻¹ and ±20 km s⁻¹, since it is a slower rotator.

The longitudinal field is a useful diagnostic of surface magnetic features, whose appearance on the visible disk is modulated at the stellar rotation period. For this reason, practical information can be extracted from the time series of $B_l$ values with standard techniques (Folsom et al. 2018; Petit et al. 2021; Carmona et al. 2023; Fouqué et al. 2023; Donati et al. 2023b). We applied a Generalised Lomb-Scargle periodogram (Zechmeister & Kürster 2009) to the entire time series as well as the individual subsets of data for all the examined stars. The results are illustrated in Fig. 2.

For EV Lac and DS Leo, the peak at the expected rotation period is unambiguously higher than a false alarm probability (FAP) of 0.01%, which emphasises the dominant stellar activity signal and confirms $B_l$ as suitable diagnostic. For CN Leo, the interpretation of the periodogram is more elaborate because of strong aliases of the observational cadence (~1 d), the gap between instrument runs (~30 d), and the time span of the time series (~1000 d). By ignoring these features, the periodogram shows a clear and unique peak at the expected stellar rotation period, for which there is no correspondence with the observing window function. If we restrict the periodogram analysis to each of the four subsets of CN Leo (see Sect. 2), only the last one 2021b2022a does not show a spike at the expected rotation period. For all stars, we retrieve $P_{rot}$ in agreement to known values from the literature: $P_{rot,EVLac}$ = 4.36 ± 0.01 d (Morin et al. 2008), $P_{rot,DSLeo}$ = 13.91 ± 0.01 d (Hébrard et al. 2016), and $P_{rot,CNLeo}$ = 2.70 ± 0.01 d (Reinhold & Hekker 2020; Lafarga et al. 2021; Irving et al. 2023).

The longitudinal field computation was carried out for both optical and near-infrared domains for all stars, and is shown in Figs. 3–5. The journal of $B_l$ values is reported in Appendix E.

## 3.1. EV Lac

EV Lac's near-infrared measurements lie between 270 and -470 G, with a mean of -50 G and a median error bar of 19 G. In comparison, the optical data have a larger spread, ranging between 380 and -540 G with a mean of -90 G and a median error bar of 16 G. When a 3000 K, near-infrared mask is used to compute the LSD profiles, the range of $B_l$ spans between -480 and 286 G, with a mean of -52 G and a median error bar of 21 G. The Stokes *I* profiles derived with the colder mask are on average 20% deeper than those associated with the 3500 K mask, and the amplitude of Stokes *V* profiles is on average ~20% larger. These differences in the Stokes profiles cancel out according to Eq. (4), ultimately making the results consistent. Overall, the near-infrared $B_l$ values seem to progressively restore the amplitude with time, after a potential minimum in 2016, making the overall trend look like a beating signal. However, the dearth of data between 2010 and 2019 prevents us from constraining the phenomenon further.

When EV Lac's near-infrared data are phase-folded at $P_{rot,EVLac}$ = 4.36 d, we note a clear rotational modulation for the three epochs (see Fig. 3). The phase variations are modelled with a Levenberg-Marquardt least squares sine fit, and we observe an increasing amplitude of the field variations, indicating a decrease in the field axisymmetry or a field enhancement

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Generalised Lomb-Scargle periodogram of the longitudinal field near-infrared time series, from top to bottom EV Lac, DS Leo, and CN Leo. In each panel the rotation period corresponding to the highest peak is highlighted with a green dotted line, while two FAP levels (0.1% and 0.01%) are shown as grey horizontal lines. The window function of the entire time series is included and mirrored with respect to the x-axis (VanderPlas 2018) to highlight aliases due to the observation cadence and scheduling.</figcaption></figure>

over time. A detailed investigation will be performed with Zeeman-Doppler imaging in Sect. 6.

## 3.2. DS Leo

DS Leo is the least magnetically active star of our sample. The longitudinal field near-infrared measurements fall within ±56 G
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Temporal evolution of the longitudinal magnetic field for EV Lac. The shape of the data points distinguishes optical (circles) from near-infrared (squares), and the colour represents the epoch in which the data were collected. Top: full time series of measurements with ESPaDOnS, NARVAL, and SPIRou. Bottom: phase-folded curves of SPIRou data points colour-coded by epoch; the associated least-squares sine fits are shown as dashed lines. The rotation period used is listed in Table 1.

with a median error bar of 7 G, whereas the optical ones range between -25 and 80 G, with a median error of 7 G (excluding a visible outlier at more than 150 G). The dispersion of the near-infrared data is comparable to the optical (20 G against 24 G) and there is a shift of the near-infrared mean towards zero by 16 G relative to the optical one. The phase-folded curves of the two epochs (2020b2021a and 2021b2022a) do not reveal particular changes, as the variations are both characterised by sinusoidal trends and similar amplitudes (see Fig. 4). The results do not change appreciably when a 4000 K synthetic mask is used against the 3500 K. There is a ~40% difference in Stokes *I* depth and amplitude of Stokes *V* between the two masks, which again cancels out in Eq. (4). The range of *B*<sub>l</sub> values with the hotter mask is between -50 and 60 G with a median error bar of 8 G.

### 3.3. CN Leo

For CN Leo, the $B_l$ measurements range from -630 to -240 G with a median error bar of 40 G, implying that we observe only the negative polarity of the large-scale field, similar to AD Leo (Bellotti et al. 2023b). These are compatible with the average $B_l$ of -691 ± 54 G reported in Morin et al. (2010). The field values exhibit a sine-like oscillation of about 1000 d as well as fluctuations within each epoch. Although the scatter of the global and epoch-by-epoch time series is lower or compatible with the

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

**Fig. 4.** Temporal evolution of the longitudinal magnetic field for DS Leo. The format is the same as in Fig. 3.

respective median error bar, which would make the oscillations dubious, the time series is dense, hence we are more sensitive to the epoch-to-epoch variability (see Fig. 5). If we bin the full time series with 50-d intervals, we are still able to capture a sinusoidal variation of 1000 d. When we phase-fold the four epochs at the stellar rotation period, we do not note any evident rotational modulation (the amplitude of the variations is consistent with zero within uncertainties), so we expect the field topology to be a dipole with negative polarity, and with the magnetic axis almost aligned with the rotation axis of the star.

## 4. FWHM of Stokes *I*

Since the Zeeman effect is proportional to Landé factor, modulus of the magnetic field (*B*), and the square of the spectral line wavelength, the width of Stokes *I* profiles computed in the near-infrared domain for stars with intense magnetic fields can be used as a proxy for Zeeman broadening measurements to first order (Donati et al. 2023a; Bellotti et al. 2023b). The value *B* encapsulates the magnetic field at both small and large scales, while the longitudinal field probes the large-scales only. Thus, the temporal analysis of the FWHM of Stokes *I* and its correlation to $B_l$ can shed more light on the link between the two distinct spatial scales. In addition, the FWHM may exhibit a long-term evolution that correlates with the large-scale field evolution, as shown for AD Leo by Bellotti et al. (2023b).

The FWHM analysis can inform the design of activity-filtering techniques for radial velocity searches of exoplanets. Haywood et al. (2022) show that, for the Sun, the rotationally-modulated line broadening correlates with the azimuthal distribution of the small-scale magnetic flux, and Lienhard et al.
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Temporal evolution of the longitudinal magnetic field for EV Lac. The shape of the data points distinguishes optical (circles) from near-infrared (squares), and the colour represents the epoch in which the data were collected. Top: full time series of measurements with ESPaDOnS, NARVAL, and SPIRou. Bottom: phase-folded curves of SPIRou data points colour-coded by epoch; the associated least-squares sine fits are shown as dashed lines. The rotation period used is listed in Table 1.

with a median error bar of 7 G, whereas the optical ones range between −25 and 80 G, with a median error of 7 G (excluding a visible outlier at more than 150 G). The dispersion of the near-infrared data is comparable to the optical (20 G against 24 G) and there is a shift of the near-infrared mean towards zero by 16 G relative to the optical one. The phase-folded curves of the two epochs (2020b2021a and 2021b2022a) do not reveal particular changes, as the variations are both characterised by sinusoidal trends and similar amplitudes (see Fig. 4). The results do not change appreciably when a 4000 K synthetic mask is used against the 3500 K. There is a ~40% difference in Stokes *I* depth and amplitude of Stokes *V* between the two masks, which again cancels out in Eq. (4). The range of *B*<sub>l</sub> values with the hotter mask is between −50 and 60 G with a median error bar of 8 G.

### 3.3. CN Leo

For CN Leo, the $B_l$ measurements range from −630 to −240 G with a median error bar of 40 G, implying that we observe only the negative polarity of the large-scale field, similar to AD Leo (Bellotti et al. 2023b). These are compatible with the average $B_l$ of −691 ± 54 G reported in Morin et al. (2010). The field values exhibit a sine-like oscillation of about 1000 d as well as fluctuations within each epoch. Although the scatter of the global and epoch-by-epoch time series is lower or compatible with the

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

**Fig. 4.** Temporal evolution of the longitudinal magnetic field for DS Leo. The format is the same as in Fig. 3.

respective median error bar, which would make the oscillations dubious, the time series is dense, hence we are more sensitive to the epoch-to-epoch variability (see Fig. 5). If we bin the full time series with 50-d intervals, we are still able to capture a sinusoidal variation of 1000 d. When we phase-fold the four epochs at the stellar rotation period, we do not note any evident rotational modulation (the amplitude of the variations is consistent with zero within uncertainties), so we expect the field topology to be a dipole with negative polarity, and with the magnetic axis almost aligned with the rotation axis of the star.

## 4. FWHM of Stokes *I*

Since the Zeeman effect is proportional to Landé factor, modulus of the magnetic field (*B*), and the square of the spectral line wavelength, the width of Stokes *I* profiles computed in the near-infrared domain for stars with intense magnetic fields can be used as a proxy for Zeeman broadening measurements to first order (Donati et al. 2023a; Bellotti et al. 2023b). The value *B* encapsulates the magnetic field at both small and large scales, while the longitudinal field probes the large-scales only. Thus, the temporal analysis of the FWHM of Stokes *I* and its correlation to $B_l$ can shed more light on the link between the two distinct spatial scales. In addition, the FWHM may exhibit a long-term evolution that correlates with the large-scale field evolution, as shown for AD Leo by Bellotti et al. (2023b).

The FWHM analysis can inform the design of activity-filtering techniques for radial velocity searches of exoplanets. Haywood et al. (2022) show that, for the Sun, the rotationally-modulated line broadening correlates with the azimuthal distribution of the small-scale magnetic flux, and Lienhard et al.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Temporal evolution of the longitudinal magnetic field for CN Leo. The format is the same as in Fig. 3.</figcaption></figure>

(2023) report the efficient activity-filtering in solar radial veloc-
ity time series using a novel proxy based on the total unsigned
magnetic flux. For AU Mic Klein et al. (2021) report the cor-
relation between the radial velocity and FWHM variations,
when the latter were computed with magnetically sensitive lines.
In Appendix B, we describe first-order computations of the
unsigned magnetic field using the FWHM as a proxy.

Following Bellotti et al. (2023b), we applied LSD using
magnetically sensitive ($g_{eff} > 1.2$) and weakly-sensitive ($g_{eff} <
1.2$) line masks. In the near-infrared domain, the masks contain
420 and 400 lines for the 3500 K case (EV Lac and DS Leo),
respectively, and for the 3000 K case (CN Leo) they contain 280
and 300 lines. In optical, the masks contain 1650 and 1500 lines
for the 3500 K case. In the following, the reference data set is the
one obtained using the full mask, for either stellar temperature
or wavelength domain. To compute the FWHM, we modelled
the Stokes *I* LSD profiles with a combination of a Voigt and
a linear model, the latter to account for residuals of continuum
normalisation.

4.1. Rotational modulation and short-term variability

We used Eq. (3) to phase-fold the FWHM time series at each
epoch and for both the full and high-$g_{\text{eff}}$ masks. With the latter
mask, Klein et al. (2021) show enhanced rotationally-modulated
variations for AU Mic, whereas Bellotti et al. (2023b) do not
report a modulation for near-infrared SPIRou observations of
AD Leo, given the large dispersion of the data set. We proceeded
in a similar manner to Bellotti et al. (2023b), and we fitted the
observed FWHM values with a constant and a sine model at the

stellar rotation period. The analysis is summarised in Table 2 and
illustrated in Fig. 6 for DS Leo and in Appendix B for EV Lac
and CN Leo.

For EV Lac, the FWHM of Stokes *I* profiles from near-
infrared observations does not exhibit evident rotational mod-
ulation (see Fig. B.2), since the change in reduced χ² between a
constant and a sine model is only marginal. In some cases, like
2021b with high-*g*<sub>eff</sub> lines, the χ² increases when using a sine
model, but it is not statistically significant (Press et al. 1992).
Only in 2020b2021a we observe a χ² improvement of about two
when using the sine model for either of the chosen masks, and we
correspondingly visualise a hint of rotational modulation, espe-
cially for the full mask. If we colour-code the high-*g*<sub>eff</sub> data by
rotational cycle, we do not observe any specific behaviour, as the
dispersion over different cycles (and within an individual epoch)
is comparable.

For the optical observations, the FWHM of EV Lac’s
Stokes *I* profiles manifests clear rotational modulation in 2006,
2007, and 2010 regardless of the mask employed (see Fig. B.2).
A plausible explanation for observing rotational modulation in
optical but not in near-infrared could be given by the larger con-
tribution of the Zeeman effect in the latter domain. The dis-
tortions induced by the Zeeman effect on the line profile shape
cannot be modelled effectively by a simple Voigt kernel, and
act as a dispersive factor of the FWHM measurement. For this
reason, the FWHM may lose its diagnostic power in the near-
infrared for strong magnetic fields.

For DS Leo, in the SPIRou 2020b2021a epoch, the phase-
folded data points vary with the stellar rotation, quantified by a
decrease in χ² from 2.8 to 0.9 (for the default mask) and from
6.2 to 2.9 (for the high-g<sub>eff</sub> mask) when using a sine model
rather than a constant (see Fig. 6). The sine model for the default
mask fits the data approximately down to noise level. In the
SPIRou 2021b2022a epoch, there is no clear rotational modu-
lation instead and the variations can be equivalently explained
by a constant model. If present, the modulation is enhanced or
quenched when high-g<sub>eff</sub> or low-g<sub>eff</sub> lines are adopted. Colour-
coding the high-g<sub>eff</sub> by cycle number reveals short-term variabil-
ity, as distinct modulations for different rotational cycles can be
extracted. This is more evident for 2020b2021a, as the varia-
tions of the data points appear ‘stacked’, that is, FWHM values
belonging to different rotational cycles exhibit a vertical offset.
This feature could be explained by the presence of differential
rotation, since the shear would have displaced magnetic regions
on the surface during the time span of our observations, mak-
ing the cycle-averaged FWHM oscillate. The quantification of
differential rotation will be investigated in Sect. 6.

For optical observations of DS Leo, we were able to capture
the rotational modulation of the FWHM in most epochs (Fig. 6),
especially when using magnetically sensitive lines. This supports
our weak-field considerations, because the contribution of Zee-
man effect to the line shape is less important for DS Leo both in
near-infrared and in optical, compared to EV Lac and CN Leo.

For CN Leo, the FWHM of near-infrared profiles does
not manifest rotational modulation, similar to EV Lac (see
Fig. B.3), as the χ² associated with a constant and a sine
model is marginally different. Colour-coding the data points
by rotational cycle does not reveal any specific feature that
could be associated to short-term variability. We observe an
increased dispersion at phases larger than 0.5 for almost all
epochs when using high-g<sub>eff</sub> lines, which could bias the mea-
surement of the average FWHM, and which we attributed to the
magnetic field, since low-g<sub>eff</sub> lines are approximately constant at
all phases. For the optical, the limited number of observations
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Temporal evolution of the longitudinal magnetic field for CN Leo. The format is the same as in Fig. 3.</figcaption></figure>

(2023) report the efficient activity-filtering in solar radial veloc-
ity time series using a novel proxy based on the total unsigned
magnetic flux. For AU Mic Klein et al. (2021) report the cor-
relation between the radial velocity and FWHM variations,
when the latter were computed with magnetically sensitive lines.
In Appendix B, we describe first-order computations of the
unsigned magnetic field using the FWHM as a proxy.

Following Bellotti et al. (2023b), we applied LSD using
magnetically sensitive ($g_{eff} > 1.2$) and weakly-sensitive ($g_{eff} <
1.2$) line masks. In the near-infrared domain, the masks contain
420 and 400 lines for the 3500 K case (EV Lac and DS Leo),
respectively, and for the 3000 K case (CN Leo) they contain 280
and 300 lines. In optical, the masks contain 1650 and 1500 lines
for the 3500 K case. In the following, the reference data set is the
one obtained using the full mask, for either stellar temperature
or wavelength domain. To compute the FWHM, we modelled
the Stokes I LSD profiles with a combination of a Voigt and
a linear model, the latter to account for residuals of continuum
normalisation.

4.1. Rotational modulation and short-term variability

We used Eq. (3) to phase-fold the FWHM time series at each
epoch and for both the full and high-$g_{\text{eff}}$ masks. With the latter
mask, Klein et al. (2021) show enhanced rotationally-modulated
variations for AU Mic, whereas Bellotti et al. (2023b) do not
report a modulation for near-infrared SPIRou observations of
AD Leo, given the large dispersion of the data set. We proceeded
in a similar manner to Bellotti et al. (2023b), and we fitted the
observed FWHM values with a constant and a sine model at the

stellar rotation period. The analysis is summarised in Table 2 and
illustrated in Fig. 6 for DS Leo and in Appendix B for EV Lac
and CN Leo.

For EV Lac, the FWHM of Stokes *I* profiles from near-infrared observations does not exhibit evident rotational modulation (see Fig. B.2), since the change in reduced $\chi^2$ between a constant and a sine model is only marginal. In some cases, like 2021b with high-$g_{\text{eff}}$ lines, the $\chi^2$ increases when using a sine model, but it is not statistically significant (Press et al. 1992). Only in 2020b2021a we observe a $\chi^2$ improvement of about two when using the sine model for either of the chosen masks, and we correspondingly visualise a hint of rotational modulation, especially for the full mask. If we colour-code the high-$g_{\text{eff}}$ data by rotational cycle, we do not observe any specific behaviour, as the dispersion over different cycles (and within an individual epoch) is comparable.

For the optical observations, the FWHM of EV Lac’s Stokes *I* profiles manifests clear rotational modulation in 2006, 2007, and 2010 regardless of the mask employed (see Fig. B.2). A plausible explanation for observing rotational modulation in optical but not in near-infrared could be given by the larger contribution of the Zeeman effect in the latter domain. The distortions induced by the Zeeman effect on the line profile shape cannot be modelled effectively by a simple Voigt kernel, and act as a dispersive factor of the FWHM measurement. For this reason, the FWHM may lose its diagnostic power in the near-infrared for strong magnetic fields.

For DS Leo, in the SPIRou 2020b2021a epoch, the phase-
folded data points vary with the stellar rotation, quantified by a
decrease in χ² from 2.8 to 0.9 (for the default mask) and from
6.2 to 2.9 (for the high-gₐeff mask) when using a sine model
rather than a constant (see Fig. 6). The sine model for the default
mask fits the data approximately down to noise level. In the
SPIRou 2021b2022a epoch, there is no clear rotational modu-
lation instead and the variations can be equivalently explained
by a constant model. If present, the modulation is enhanced or
quenched when high-gₐeff or low-gₐeff lines are adopted. Colour-
coding the high-gₐeff by cycle number reveals short-term variabil-
ity, as distinct modulations for different rotational cycles can be
extracted. This is more evident for 2020b2021a, as the varia-
tions of the data points appear ‘stacked’, that is, FWHM values
belonging to different rotational cycles exhibit a vertical offset.
This feature could be explained by the presence of differential
rotation, since the shear would have displaced magnetic regions
on the surface during the time span of our observations, mak-
ing the cycle-averaged FWHM oscillate. The quantification of
differential rotation will be investigated in Sect. 6.

For optical observations of DS Leo, we were able to capture
the rotational modulation of the FWHM in most epochs (Fig. 6),
especially when using magnetically sensitive lines. This supports
our weak-field considerations, because the contribution of Zee-
man effect to the line shape is less important for DS Leo both in
near-infrared and in optical, compared to EV Lac and CN Leo.

For CN Leo, the FWHM of near-infrared profiles does
not manifest rotational modulation, similar to EV Lac (see
Fig. B.3), as the $\chi^2$ associated with a constant and a sine
model is marginally different. Colour-coding the data points
by rotational cycle does not reveal any specific feature that
could be associated to short-term variability. We observe an
increased dispersion at phases larger than 0.5 for almost all
epochs when using high-$g_{\text{eff}}$ lines, which could bias the mea-
surement of the average FWHM, and which we attributed to the
magnetic field, since low-$g_{\text{eff}}$ lines are approximately constant at
all phases. For the optical, the limited number of observations
<table><caption>**Table 2.** Comparison of a constant line against a sine fit for the FWHM phase variations of EV Lac, DS Leo, and CN Leo.</caption><thead><tr><th>Epoch</th><th>Mask</th><th>Mean<br>[km s<sup>-1</sup>]</th><th>Mean error<br>[km s<sup>-1</sup>]</th><th>STD<br>[km s<sup>-1</sup>]</th><th>RMS<sub>const</sub><br>[km s<sup>-1</sup>]</th><th>χ<sup>2</sup><sub>r,const</sub></th><th>RMS<sub>sine</sub><br>[km s<sup>-1</sup>]</th><th>χ<sup>2</sup><sub>r,sine</sub></th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-weight:bold;">EV Lac</td></tr><tr><td rowspan="2">2019b</td><td>Default</td><td>22</td><td>0.65</td><td>1.68</td><td>1.68</td><td>6.9</td><td>1.45</td><td>5.9</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>29</td><td>1.12</td><td>4.39</td><td>4.39</td><td>19.7</td><td>4.17</td><td>19.7</td></tr><tr><td rowspan="2">2020b2021a</td><td>Default</td><td>21</td><td>0.44</td><td>1.00</td><td>1.00</td><td>6.2</td><td>0.82</td><td>4.4</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>28</td><td>0.87</td><td>2.62</td><td>2.62</td><td>12.9</td><td>2.45</td><td>10.7</td></tr><tr><td rowspan="2">2021b</td><td>Default</td><td>21</td><td>0.46</td><td>0.98</td><td>0.98</td><td>4.8</td><td>0.90</td><td>4.5</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>28</td><td>0.74</td><td>2.57</td><td>2.57</td><td>13.0</td><td>2.44</td><td>16.1</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">DS Leo</td></tr><tr><td rowspan="2">2020b2021a</td><td>Default</td><td>13</td><td>0.18</td><td>0.29</td><td>0.29</td><td>2.8</td><td>0.16</td><td>0.9</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>14</td><td>0.21</td><td>0.47</td><td>0.47</td><td>6.2</td><td>0.33</td><td>2.9</td></tr><tr><td rowspan="2">2021b2022a</td><td>Default</td><td>13</td><td>0.18</td><td>0.21</td><td>0.21</td><td>1.5</td><td>0.20</td><td>1.3</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>14</td><td>0.20</td><td>0.32</td><td>0.32</td><td>3.1</td><td>0.30</td><td>2.7</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">CN Leo</td></tr><tr><td rowspan="2">2019a</td><td>Default</td><td>15</td><td>0.34</td><td>1.27</td><td>1.27</td><td>16.7</td><td>1.22</td><td>19.7</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>21</td><td>1.19</td><td>3.73</td><td>3.73</td><td>19.8</td><td>3.44</td><td>22.5</td></tr><tr><td rowspan="2">2019b2020a</td><td>Default</td><td>17</td><td>0.36</td><td>1.41</td><td>1.41</td><td>19.1</td><td>1.35</td><td>16.6</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>25</td><td>1.57</td><td>4.44</td><td>4.44</td><td>15.3</td><td>4.27</td><td>13.6</td></tr><tr><td rowspan="2">2020b2021a</td><td>Default</td><td>16</td><td>0.37</td><td>1.69</td><td>1.69</td><td>45.8</td><td>1.61</td><td>39.8</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>22</td><td>1.51</td><td>5.08</td><td>5.08</td><td>26.8</td><td>4.74</td><td>23.1</td></tr><tr><td rowspan="2">2021b2022a</td><td>Default</td><td>17</td><td>0.34</td><td>1.37</td><td>1.37</td><td>25.7</td><td>1.29</td><td>24.7</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>22</td><td>1.41</td><td>4.94</td><td>4.94</td><td>16.3</td><td>4.65</td><td>15.8</td></tr></tbody></table>

Notes. The columns are: 1) subset of the time series, 2) line list used in LSD computation, 3) mean FWHM, 4) mean error bar, 5) standard deviation of the data sets, 6) RMS (root mean square) residual of a constant line fit equal to the average of the data set, 7) reduced $\chi^2$ of a constant model, 8) RMS residual of a sine fit at the stellar rotation period, and 9) reduced $\chi^2$ of a sine model.

(four) prevents us from performing a rotational modulation analysis.

## 4.2. Correlation between FWHM and $B_l$

The longitudinal magnetic field is a proxy for large-scale magnetic field, and the FWHM is a proxy for the small-scale field, so we inspected their correlations to yield more insights on the link between these two spatial scales. In particular, we computed Pearson correlation coefficients via 5000 bootstrap iterations. The values summarised in Table 3 correspond to the mean and standard deviation of the coefficient distribution. The analysis was carried out for the FWHM computed with the three masks: full, low-$g_{eff}$, and high-$g_{eff}$. The analysis is illustrated in Fig. 7 for optical observations of EV Lac and Fig. B.4 for optical and near-infrared observations of DS Leo. CN Leo was excluded because the FWHM from near-infrared lines is not modulated by rotation.

For EV Lac, there is a positive correlation for 2006 data, which increases when going from low-$g_{eff}$ lines (Pearson $R = 0.19$) to high-$g_{eff}$ lines (Pearson $R = 0.59$). This is the expected behaviour from the principle of Zeeman effect: the line profile is broader when the field is stronger. A positive correlation is also seen in 2007, but its strength does not vary according to the line mask adopted.

For DS Leo, we observe a variety of relations: there is a positive correlation that increases going from low-$g_{eff}$ to high-$g_{eff}$ lines in 2007 and 2011, whereas in 2012 and 2014, there is an anti-correlation that becomes stronger for low-$g_{eff}$ lines. For 2008 and 2010, there is one outlying data point at large $B_l$ values that likely biases the correlation coefficient, as the bulk of the data points does not manifest particular correlations. When looking at the SPIRou 2020b2021a epoch, we do not observe

significant correlation, the Pearson $R$ coefficient being smaller than 0.2 for the three masks.

While not straightforward to interpret, these results altogether suggest an underlying complexity in the link between small- and large-scale field, which may not be encapsulated in a simple scaling. The variety of correlations for different epochs for DS Leo may stem from the fact that, for early M dwarfs, the ratio of magnetic field recovered with polarised light relative to unpolarised light is only a few percent (and it increases to dozens of percents for mid-M, Reiners & Basri 2009; Morin et al. 2010; Kochukhov 2021), hence there is an even less obvious scaling between small- and large-scale fields.

## 4.3. Long-term variations of FWHM

We inspected the long-term behaviour of the epoch-averaged FWHM of Stokes $I$, when computed with the full, high-$g_{eff}$, and low-$g_{eff}$ line masks. For AD Leo, Bellotti et al. (2023b) show that the evolution of the epoch-averaged FWHM correlates with the secular trend of $B_l$ and the magnetic flux measured by modelling the Zeeman broadening and intensification. The results for the stars examined are shown in Fig. 8.

For EV Lac, no clear trend is observed. In near-infrared, the average FWHM decreased from 22 to 21 km s$^{-1}$ for the full mask, while it was stable around 16 km s$^{-1}$ for the low-$g_{eff}$ mask. For high-$g_{eff}$, there is a slight temporal decrease from 29 km s$^{-1}$ to just below 28 km s$^{-1}$. In optical, we observe the same slight oscillation when adopting the three line masks, around 14 km s$^{-1}$, above 12 km s$^{-1}$ and below 12 km s$^{-1}$ for the high-$g_{eff}$, full, and low-$g_{eff}$ masks, respectively.

For DS Leo, the mean FWHM from near-infrared observations is 13 km s$^{-1}$, 14 km s$^{-1}$ and 11 km s$^{-1}$ for the full, high-$g_{eff}$
<table><caption>**Table 2.** Comparison of a constant line against a sine fit for the FWHM phase variations of EV Lac, DS Leo, and CN Leo.</caption><thead><tr><th>Epoch</th><th>Mask</th><th>Mean<br>[km s<sup>-1</sup>]</th><th>Mean error<br>[km s<sup>-1</sup>]</th><th>STD<br>[km s<sup>-1</sup>]</th><th>RMS<sub>const</sub><br>[km s<sup>-1</sup>]</th><th>χ<sup>2</sup><sub>r,const</sub></th><th>RMS<sub>sine</sub><br>[km s<sup>-1</sup>]</th><th>χ<sup>2</sup><sub>r,sine</sub></th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-weight:bold;">EV Lac</td></tr><tr><td rowspan="2">2019b</td><td>Default</td><td>22</td><td>0.65</td><td>1.68</td><td>1.68</td><td>6.9</td><td>1.45</td><td>5.9</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>29</td><td>1.12</td><td>4.39</td><td>4.39</td><td>19.7</td><td>4.17</td><td>19.7</td></tr><tr><td rowspan="2">2020b2021a</td><td>Default</td><td>21</td><td>0.44</td><td>1.00</td><td>1.00</td><td>6.2</td><td>0.82</td><td>4.4</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>28</td><td>0.87</td><td>2.62</td><td>2.62</td><td>12.9</td><td>2.45</td><td>10.7</td></tr><tr><td rowspan="2">2021b</td><td>Default</td><td>21</td><td>0.46</td><td>0.98</td><td>0.98</td><td>4.8</td><td>0.90</td><td>4.5</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>28</td><td>0.74</td><td>2.57</td><td>2.57</td><td>13.0</td><td>2.44</td><td>16.1</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">DS Leo</td></tr><tr><td rowspan="2">2020b2021a</td><td>Default</td><td>13</td><td>0.18</td><td>0.29</td><td>0.29</td><td>2.8</td><td>0.16</td><td>0.9</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>14</td><td>0.21</td><td>0.47</td><td>0.47</td><td>6.2</td><td>0.33</td><td>2.9</td></tr><tr><td rowspan="2">2021b2022a</td><td>Default</td><td>13</td><td>0.18</td><td>0.21</td><td>0.21</td><td>1.5</td><td>0.20</td><td>1.3</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>14</td><td>0.20</td><td>0.32</td><td>0.32</td><td>3.1</td><td>0.30</td><td>2.7</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">CN Leo</td></tr><tr><td rowspan="2">2019a</td><td>Default</td><td>15</td><td>0.34</td><td>1.27</td><td>1.27</td><td>16.7</td><td>1.22</td><td>19.7</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>21</td><td>1.19</td><td>3.73</td><td>3.73</td><td>19.8</td><td>3.44</td><td>22.5</td></tr><tr><td rowspan="2">2019b2020a</td><td>Default</td><td>17</td><td>0.36</td><td>1.41</td><td>1.41</td><td>19.1</td><td>1.35</td><td>16.6</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>25</td><td>1.57</td><td>4.44</td><td>4.44</td><td>15.3</td><td>4.27</td><td>13.6</td></tr><tr><td rowspan="2">2020b2021a</td><td>Default</td><td>16</td><td>0.37</td><td>1.69</td><td>1.69</td><td>45.8</td><td>1.61</td><td>39.8</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>22</td><td>1.51</td><td>5.08</td><td>5.08</td><td>26.8</td><td>4.74</td><td>23.1</td></tr><tr><td rowspan="2">2021b2022a</td><td>Default</td><td>17</td><td>0.34</td><td>1.37</td><td>1.37</td><td>25.7</td><td>1.29</td><td>24.7</td></tr><tr><td><i>g</i><sub>eff</sub> > 1.2</td><td>22</td><td>1.41</td><td>4.94</td><td>4.94</td><td>16.3</td><td>4.65</td><td>15.8</td></tr></tbody></table>

Notes. The columns are: 1) subset of the time series, 2) line list used in LSD computation, 3) mean FWHM, 4) mean error bar, 5) standard deviation of the data sets, 6) RMS (root mean square) residual of a constant line fit equal to the average of the data set, 7) reduced $\chi^2$ of a constant model, 8) RMS residual of a sine fit at the stellar rotation period, and 9) reduced $\chi^2$ of a sine model.

(four) prevents us from performing a rotational modulation analysis.

## 4.2. Correlation between FWHM and $B_l$

The longitudinal magnetic field is a proxy for large-scale magnetic field, and the FWHM is a proxy for the small-scale field, so we inspected their correlations to yield more insights on the link between these two spatial scales. In particular, we computed Pearson correlation coefficients via 5000 bootstrap iterations. The values summarised in Table 3 correspond to the mean and standard deviation of the coefficient distribution. The analysis was carried out for the FWHM computed with the three masks: full, low-$g_{eff}$, and high-$g_{eff}$. The analysis is illustrated in Fig. 7 for optical observations of EV Lac and Fig. B.4 for optical and near-infrared observations of DS Leo. CN Leo was excluded because the FWHM from near-infrared lines is not modulated by rotation.

For EV Lac, there is a positive correlation for 2006 data, which increases when going from low-$g_{eff}$ lines (Pearson $R = 0.19$) to high-$g_{eff}$ lines (Pearson $R = 0.59$). This is the expected behaviour from the principle of Zeeman effect: the line profile is broader when the field is stronger. A positive correlation is also seen in 2007, but its strength does not vary according to the line mask adopted.

For DS Leo, we observe a variety of relations: there is a positive correlation that increases going from low-$g_{eff}$ to high-$g_{eff}$ lines in 2007 and 2011, whereas in 2012 and 2014, there is an anti-correlation that becomes stronger for low-$g_{eff}$ lines. For 2008 and 2010, there is one outlying data point at large $B_l$ values that likely biases the correlation coefficient, as the bulk of the data points does not manifest particular correlations. When looking at the SPIRou 2020b2021a epoch, we do not observe

significant correlation, the Pearson $R$ coefficient being smaller than 0.2 for the three masks.

While not straightforward to interpret, these results altogether suggest an underlying complexity in the link between small- and large-scale field, which may not be encapsulated in a simple scaling. The variety of correlations for different epochs for DS Leo may stem from the fact that, for early M dwarfs, the ratio of magnetic field recovered with polarised light relative to unpolarised light is only a few percent (and it increases to dozens of percents for mid-M, Reiners & Basri 2009; Morin et al. 2010; Kochukhov 2021), hence there is an even less obvious scaling between small- and large-scale fields.

## 4.3. Long-term variations of FWHM

We inspected the long-term behaviour of the epoch-averaged FWHM of Stokes $I$, when computed with the full, high-$g_{eff}$, and low-$g_{eff}$ line masks. For AD Leo, Bellotti et al. (2023b) show that the evolution of the epoch-averaged FWHM correlates with the secular trend of $B_l$ and the magnetic flux measured by modelling the Zeeman broadening and intensification. The results for the stars examined are shown in Fig. 8.

For EV Lac, no clear trend is observed. In near-infrared, the average FWHM decreased from 22 to 21 km s$^{-1}$ for the full mask, while it was stable around 16 km s$^{-1}$ for the low-$g_{eff}$ mask. For high-$g_{eff}$, there is a slight temporal decrease from 29 km s$^{-1}$ to just below 28 km s$^{-1}$. In optical, we observe the same slight oscillation when adopting the three line masks, around 14 km s$^{-1}$, above 12 km s$^{-1}$ and below 12 km s$^{-1}$ for the high-$g_{eff}$, full, and low-$g_{eff}$ masks, respectively.

For DS Leo, the mean FWHM from near-infrared observations is 13 km s$^{-1}$, 14 km s$^{-1}$ and 11 km s$^{-1}$ for the full, high-$g_{eff}$
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Rotational modulation analysis of FWHM for DS Leo. The top two panels correspond to the near-infrared epochs (2020b2021a and 2021b2022a), whereas the remaining panels correspond to the optical epochs (2007, 2008, 2010, 2011, 2012, and 2014). In each panel are shown the phase-folded time series of FWHM computed with the full, low-*g*<sub>eff</sub>, and high-*g*<sub>eff</sub> mask with a horizontal line representing the mean FWHM value (also reported in each legend). For the near-infrared epochs the high-*g*<sub>eff</sub> time series is colour-coded by rotational cycle to inspect signs of short-term variability. Observed stacked data points: observations within the first five cycles tend to fall at high FWHM values, then progressively at lower FWHM between cycle number 5 and 11, and at higher FWHM values again in the latest cycles.</figcaption></figure>

and low-*g*<sub>eff</sub> mask and for both epochs. In optical, the mean FWHM is stable around 8.0 km s<sup>-1</sup>, 9.0 km s<sup>-1</sup> and 8.4 km s<sup>-1</sup> for the full, high-*g*<sub>eff</sub> and low-*g*<sub>eff</sub> masks. Therefore, we do

not observe a significant evolution. The same conclusion is reached when computing near-infrared LSD profiles with a 4000 K mask.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Rotational modulation analysis of FWHM for DS Leo. The top two panels correspond to the near-infrared epochs (2020b2021a and 2021b2022a), whereas the remaining panels correspond to the optical epochs (2007, 2008, 2010, 2011, 2012, and 2014). In each panel are shown the phase-folded time series of FWHM computed with the full, low-*g*<sub>eff</sub>, and high-*g*<sub>eff</sub> mask with a horizontal line representing the mean FWHM value (also reported in each legend). For the near-infrared epochs the high-*g*<sub>eff</sub> time series is colour-coded by rotational cycle to inspect signs of short-term variability. Observed stacked data points: observations within the first five cycles tend to fall at high FWHM values, then progressively at lower FWHM between cycle number 5 and 11, and at higher FWHM values again in the latest cycles.</figcaption></figure>

and low-*g*<sub>eff</sub> mask and for both epochs. In optical, the mean FWHM is stable around 8.0 km s<sup>-1</sup>, 9.0 km s<sup>-1</sup> and 8.4 km s<sup>-1</sup> for the full, high-*g*<sub>eff</sub> and low-*g*<sub>eff</sub> masks. Therefore, we do

not observe a significant evolution. The same conclusion is reached when computing near-infrared LSD profiles with a 4000 K mask.
<table><caption>**Table 3.** Correlation coefficients between $|B_l|$ and FWHM of EV Lac and DS Leo, for the epochs in which we observed rotational modulation at the stellar rotation period.</caption><thead><tr><th>Epoch</th><th>Full</th><th>High-g<sub>eff</sub></th><th>Low-g<sub>eff</sub></th></tr></thead><tbody><tr><td colspan="4">EV Lac</td></tr><tr><td>2006</td><td>0.33 ± 0.39</td><td>0.49 ± 0.36</td><td>0.08 ± 0.48</td></tr><tr><td>2007</td><td>0.27 ± 0.29</td><td>0.30 ± 0.28</td><td>0.29 ± 0.30</td></tr><tr><td colspan="4">DS Leo</td></tr><tr><td>2007</td><td>0.24 ± 0.26</td><td>0.10 ± 0.29</td><td>0.34 ± 0.23</td></tr><tr><td>2008</td><td>-0.35 ± 0.15</td><td>-0.34 ± 0.15</td><td>-0.51 ± 0.12</td></tr><tr><td>2010</td><td>-0.07 ± 0.33</td><td>0.19 ± 0.24</td><td>-0.49 ± 0.42</td></tr><tr><td>2011</td><td>0.84 ± 0.12</td><td>0.86 ± 0.09</td><td>0.20 ± 0.28</td></tr><tr><td>2012</td><td>-0.51 ± 0.25</td><td>0.27 ± 0.29</td><td>-0.66 ± 0.18</td></tr><tr><td>2014</td><td>-0.38 ± 0.20</td><td>-0.15 ± 0.23</td><td>-0.55 ± 0.19</td></tr><tr><td>2020b2021a</td><td>0.09 ± 0.12</td><td>0.20 ± 0.12</td><td>-0.11 ± 0.11</td></tr></tbody></table>

**Notes.** The columns are: 1) subset of the time series, 2) average Pearson correlation coefficient from bootstrap analysis using the full LSD mask, 3) the high-$g_{eff}$ mask, and 4) the low-$g_{eff}$ mask.

For CN Leo, the low-$g_{eff}$ lines are stable around 12 km s$^{-1}$, the full mask values evolve from 15 km s$^{-1}$ to 17 km s$^{-1}$ and then back to 16 km s$^{-1}$, while the high-$g_{eff}$ mask values from 21 km s$^{-1}$ to 25 km s$^{-1}$ and back to 22 km s$^{-1}$. Such temporal variation is anti-correlated with the long-term sine oscillation of $B_l$ (in absolute value), the average FWHM being largest when the $B_l$ is smallest.

## 5. Principal component analysis

Lehmann & Donati (2022) presented a method based on principal component analysis (PCA) that reveals key properties of the stellar magnetic topology and its temporal evolution directly from the LSD Stokes V profiles, without prior assumptions about stellar parameters such as $v_{eq} \sin i$ and inclination. The PCA method provides information about the degree of axisymmetry, the poloidal-to-toroidal fraction of the axisymmetric field, the field complexity (i.e., dipolar or more complex topology), and the temporal evolution of these three parameters. Information about the axisymmetric field are captured by the mean Stokes V profile of the observed time series and information about the non-axisymmetric field in the mean-subtracted Stokes V profiles. The PCA method analyses the mean Stokes V profile to infer the degree of axisymmetry and whether the axisymmetric field is more poloidal or toroidal. Furthermore, by determining the eigenvectors and coefficients of the mean-subtracted Stokes V profiles using PCA, information about the field complexity and its evolution with time can be obtained. For more information about the PCA method see Lehmann & Donati (2022) and Lehmann et al. (2024).

We studied the near-infrared LSD Stokes V time series collected by SPIRou for EV Lac, DS Leo and CN Leo similar to our analysis for AD Leo (Bellotti et al. 2023b). As in Lehmann et al. (2024), we used the S/N-weighted mean profiles and the weighted PCA analysis (Delchambre 2015), as they provide better results for long time series with varying S/N.

### 5.1. EV Lac

The mean Stokes V profile of EV Lac is antisymmetric to the line centre, suggesting that the axisymmetric field is predominantly poloidal (see Fig. 9a), and begins to split at the line centre owing

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Correlation analysis between FWHM and $|B_l|$ for EV Lac. The three panels correspond to the optical epochs in which rotational modulation is present: 2006, 2007, and 2010. In all panels the data points are colour-coded based on the line mask used for LSD: full (black), low-$g_{eff}$ (purple), and high-$g_{eff}$ (yellow).</figcaption></figure>

to the strong magnetic field of EV Lac. The PCA method is based on the weak-field approximation valid for $B \ll 1$ kG. For EV Lac the total unsigned magnetic field values are above this limit with 4.3 kG (Shulyak et al. 2017), which causes the visible splitting at line centre of the mean Stokes V profile. Therefore, we need to be cautious about our conclusions derived with the PCA method. This also applies for CN Leo which exhibits a total field value of ~2.3 kG (Shulyak et al. 2019), although we do not record a clear Stokes V splitting. As a double-check, we simulated the Stokes V profile using synthetic and observed magnetic field maps while artificially increasing or decreasing the magnetic field strength. Our simulations for the three stars of our
<table><caption>**Table 3.** Correlation coefficients between $|B_l|$ and FWHM of EV Lac and DS Leo, for the epochs in which we observed rotational modulation at the stellar rotation period.</caption><thead><tr><th>Epoch</th><th>Full</th><th>High-g<sub>eff</sub></th><th>Low-g<sub>eff</sub></th></tr></thead><tbody><tr><td colspan="4">EV Lac</td></tr><tr><td>2006</td><td>0.33 ± 0.39</td><td>0.49 ± 0.36</td><td>0.08 ± 0.48</td></tr><tr><td>2007</td><td>0.27 ± 0.29</td><td>0.30 ± 0.28</td><td>0.29 ± 0.30</td></tr><tr><td colspan="4">DS Leo</td></tr><tr><td>2007</td><td>0.24 ± 0.26</td><td>0.10 ± 0.29</td><td>0.34 ± 0.23</td></tr><tr><td>2008</td><td>-0.35 ± 0.15</td><td>-0.34 ± 0.15</td><td>-0.51 ± 0.12</td></tr><tr><td>2010</td><td>-0.07 ± 0.33</td><td>0.19 ± 0.24</td><td>-0.49 ± 0.42</td></tr><tr><td>2011</td><td>0.84 ± 0.12</td><td>0.86 ± 0.09</td><td>0.20 ± 0.28</td></tr><tr><td>2012</td><td>-0.51 ± 0.25</td><td>0.27 ± 0.29</td><td>-0.66 ± 0.18</td></tr><tr><td>2014</td><td>-0.38 ± 0.20</td><td>-0.15 ± 0.23</td><td>-0.55 ± 0.19</td></tr><tr><td>2020b2021a</td><td>0.09 ± 0.12</td><td>0.20 ± 0.12</td><td>-0.11 ± 0.11</td></tr></tbody></table>

**Notes.** The columns are: 1) subset of the time series, 2) average Pearson correlation coefficient from bootstrap analysis using the full LSD mask, 3) the high-$g_{eff}$ mask, and 4) the low-$g_{eff}$ mask.

For CN Leo, the low-$g_{eff}$ lines are stable around 12 km s$^{-1}$, the full mask values evolve from 15 km s$^{-1}$ to 17 km s$^{-1}$ and then back to 16 km s$^{-1}$, while the high-$g_{eff}$ mask values from 21 km s$^{-1}$ to 25 km s$^{-1}$ and back to 22 km s$^{-1}$. Such temporal variation is anti-correlated with the long-term sine oscillation of $B_l$ (in absolute value), the average FWHM being largest when the $B_l$ is smallest.

## 5. Principal component analysis

Lehmann & Donati (2022) presented a method based on principal component analysis (PCA) that reveals key properties of the stellar magnetic topology and its temporal evolution directly from the LSD Stokes V profiles, without prior assumptions about stellar parameters such as $v_{eq} \sin i$ and inclination. The PCA method provides information about the degree of axisymmetry, the poloidal-to-toroidal fraction of the axisymmetric field, the field complexity (i.e., dipolar or more complex topology), and the temporal evolution of these three parameters. Information about the axisymmetric field are captured by the mean Stokes V profile of the observed time series and information about the non-axisymmetric field in the mean-subtracted Stokes V profiles. The PCA method analyses the mean Stokes V profile to infer the degree of axisymmetry and whether the axisymmetric field is more poloidal or toroidal. Furthermore, by determining the eigenvectors and coefficients of the mean-subtracted Stokes V profiles using PCA, information about the field complexity and its evolution with time can be obtained. For more information about the PCA method see Lehmann & Donati (2022) and Lehmann et al. (2024).

We studied the near-infrared LSD Stokes V time series collected by SPIRou for EV Lac, DS Leo and CN Leo similar to our analysis for AD Leo (Bellotti et al. 2023b). As in Lehmann et al. (2024), we used the S/N-weighted mean profiles and the weighted PCA analysis (Delchambre 2015), as they provide better results for long time series with varying S/N.

### 5.1. EV Lac

The mean Stokes V profile of EV Lac is antisymmetric to the line centre, suggesting that the axisymmetric field is predominantly poloidal (see Fig. 9a), and begins to split at the line centre owing

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Correlation analysis between FWHM and $|B_l|$ for EV Lac. The three panels correspond to the optical epochs in which rotational modulation is present: 2006, 2007, and 2010. In all panels the data points are colour-coded based on the line mask used for LSD: full (black), low-$g_{eff}$ (purple), and high-$g_{eff}$ (yellow).</figcaption></figure>

to the strong magnetic field of EV Lac. The PCA method is based on the weak-field approximation valid for $B \ll 1$ kG. For EV Lac the total unsigned magnetic field values are above this limit with 4.3 kG (Shulyak et al. 2017), which causes the visible splitting at line centre of the mean Stokes V profile. Therefore, we need to be cautious about our conclusions derived with the PCA method. This also applies for CN Leo which exhibits a total field value of ~2.3 kG (Shulyak et al. 2019), although we do not record a clear Stokes V splitting. As a double-check, we simulated the Stokes V profile using synthetic and observed magnetic field maps while artificially increasing or decreasing the magnetic field strength. Our simulations for the three stars of our
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Long-term evolution of the epoch-averaged FWHM. The panels show the results from optical (circles) and near-infrared (pluses) observations of EV Lac (top), DS Leo (middle), and CN Leo (bottom). In all panels, the data points are colour-coded based on the line mask used for LSD: full (black), low-$g_{\text{eff}}$ (purple), and high-$g_{\text{eff}}$ (yellow). The left y-axis refers to the optical observations, and the right y-axis to the near-infrared observations.</figcaption></figure>

sample show that the conclusions derived with the PCA method remained valid for $B > 1$ kG in the case of EV Lac (and CN Leo), even if it is based on the weak-field approximation.

A signal is present in the first three eigenvectors of the PCA analysis of the mean-subtracted Stokes V profile capturing the information about the non-axisymmetric component of EV Lac’s topology, with the first two displaying an antisymmetric shape and the third a symmetric shape relative to the line centre (see Fig. 9b). The large number of eigenvectors containing a signal symbolises a non-axisymmetric field topology.

Figure 9c shows the per-epoch analysis of EV Lac for the Stokes V mean profile capturing the information of the axisymmetric field and the PCA coefficients of the first three eigenvec-

tors allowing the analysis of the non-axisymmetric field. The coefficients are derived from the weighted PCA of the mean-subtracted Stokes V time series using the weighted mean profile computed over all epochs (e.g., Fig. 9a), and not the weighted mean profile of each individual epoch (e.g., Fig. 9c left column). This because, in the latter case, the mean value of the coefficients would be centred for each epoch, and the amplitudes of the variations could not be compared between epochs.

We observe that the amplitude of the mean Stokes V profile increases from 2019b to 2020b2021a. It remains antisymmetric with respect to the line centre for all epochs, which indicates a poloidal-dominated axisymmetric component for 2019b to 2021b. Since the amplitude of the coefficients also increases from 2019b to 2021a, we conclude that the field of EV Lac strengthens between 2019b and 2021a. The most complex phase variations of the coefficients occur in the first epoch 2019b, suggesting that the topology of EV Lac has the lowest dipole fraction and therefore is most complex in 2019b. Somehow opposite is the second epoch 2020b2021a: the coefficients show sinusoidal curves for all three eigenvectors, reflecting a tilted, dipole-dominated topology. In the last epoch 2021b, the coefficients still show a sinusoidal trend, but become slightly more flat between rotational phase 0.0–0.3.

## 5.2. DS Leo

Figure C.1 shows the PCA analysis for DS Leo. The mean profile reflecting the axisymmetric component of the topology is mainly symmetric to the line centre, indicating a significant toroidal field fraction (>5%) of the axisymmetric field, see Fig. C.1a. The PCA analysis of the mean-subtracted profiles allows the analysis of the non-axisymmetric component. Two eigenvectors emerge from the noise, one antisymmetric and one symmetric, indicating a non-axisymmetric topology. The degree of axisymmetry is however larger than EV Lac’s, as fewer eigenvectors show a signal (see Fig. C.1b).

The per-epoch analysis shows a decrease in amplitude of the mean Stokes V profile, as well as of the phase variations of the coefficients, indicating a weakening of the magnetic field from the first to the second epoch. For both epochs, the coefficient of the first eigenvector exhibits the typical sinusoidal trend of a tilted dipole, and the pointing phase of the dipole appears to be stable. The coefficient of the second eigenvector manifests a more complex trend that evolves from 2020b2021a to 2021b2022a. This reflects the evolution of other magnetic field features in addition to the tilted dipole. We clearly see that the coefficients also vary within the epochs, implying that the magnetic field of DS Leo evolves on the timescale of its rotation period, as already suggested by Lehmann & Donati (2022) in their analysis of the 2008 optical observations of DS Leo.

## 5.3. CN Leo

CN Leo exhibits a strong and antisymmetric Stokes V mean profile, indicating a predominantly poloidal axisymmetric field, as illustrated in Fig. C.2a. The first eigenvector is antisymmetric and there is evidence for a second, symmetric eigenvector, albeit the larger noise (Fig. C.2b). This suggests that CN Leo is the most axisymmetric M dwarf of our sample.

The per-epoch Stokes V mean profiles show a varying amplitude, together with the amplitude of the rotational modulation of the coefficients, which could be due to variability in either field strength or axisymmetry (see Fig. C.2c). The first epoch 2019a appears to be highly axisymmetric, as the coefficients cluster around zero. The second epoch 2019b2020a has a clear
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Long-term evolution of the epoch-averaged FWHM. The panels show the results from optical (circles) and near-infrared (pluses) observations of EV Lac (top), DS Leo (middle), and CN Leo (bottom). In all panels, the data points are colour-coded based on the line mask used for LSD: full (black), low-$g_{\text{eff}}$ (purple), and high-$g_{\text{eff}}$ (yellow). The left y-axis refers to the optical observations, and the right y-axis to the near-infrared observations.</figcaption></figure>

sample show that the conclusions derived with the PCA method remained valid for $B > 1$ kG in the case of EV Lac (and CN Leo), even if it is based on the weak-field approximation.

A signal is present in the first three eigenvectors of the PCA analysis of the mean-subtracted Stokes V profile capturing the information about the non-axisymmetric component of EV Lac's topology, with the first two displaying an antisymmetric shape and the third a symmetric shape relative to the line centre (see Fig. 9b). The large number of eigenvectors containing a signal symbolises a non-axisymmetric field topology.

Figure 9c shows the per-epoch analysis of EV Lac for the Stokes V mean profile capturing the information of the axisymmetric field and the PCA coefficients of the first three eigenvec-

tors allowing the analysis of the non-axisymmetric field. The coefficients are derived from the weighted PCA of the mean-subtracted Stokes V time series using the weighted mean profile computed over all epochs (e.g., Fig. 9a), and not the weighted mean profile of each individual epoch (e.g., Fig. 9c left column). This because, in the latter case, the mean value of the coefficients would be centred for each epoch, and the amplitudes of the variations could not be compared between epochs.

We observe that the amplitude of the mean Stokes V profile increases from 2019b to 2020b2021a. It remains antisymmetric with respect to the line centre for all epochs, which indicates a poloidal-dominated axisymmetric component for 2019b to 2021b. Since the amplitude of the coefficients also increases from 2019b to 2021a, we conclude that the field of EV Lac strengthens between 2019b and 2021a. The most complex phase variations of the coefficients occur in the first epoch 2019b, suggesting that the topology of EV Lac has the lowest dipole fraction and therefore is most complex in 2019b. Somehow opposite is the second epoch 2020b2021a: the coefficients show sinusoidal curves for all three eigenvectors, reflecting a tilted, dipole-dominated topology. In the last epoch 2021b, the coefficients still show a sinusoidal trend, but become slightly more flat between rotational phase 0.0–0.3.

## 5.2. DS Leo

Figure C.1 shows the PCA analysis for DS Leo. The mean profile reflecting the axisymmetric component of the topology is mainly symmetric to the line centre, indicating a significant toroidal field fraction (>5%) of the axisymmetric field, see Fig. C.1a. The PCA analysis of the mean-subtracted profiles allows the analysis of the non-axisymmetric component. Two eigenvectors emerge from the noise, one antisymmetric and one symmetric, indicating a non-axisymmetric topology. The degree of axisymmetry is however larger than EV Lac's, as fewer eigenvectors show a signal (see Fig. C.1b).

The per-epoch analysis shows a decrease in amplitude of the mean Stokes V profile, as well as of the phase variations of the coefficients, indicating a weakening of the magnetic field from the first to the second epoch. For both epochs, the coefficient of the first eigenvector exhibits the typical sinusoidal trend of a tilted dipole, and the pointing phase of the dipole appears to be stable. The coefficient of the second eigenvector manifests a more complex trend that evolves from 2020b2021a to 2021b2022a. This reflects the evolution of other magnetic field features in addition to the tilted dipole. We clearly see that the coefficients also vary within the epochs, implying that the magnetic field of DS Leo evolves on the timescale of its rotation period, as already suggested by Lehmann & Donati (2022) in their analysis of the 2008 optical observations of DS Leo.

## 5.3. CN Leo

CN Leo exhibits a strong and antisymmetric Stokes V mean profile, indicating a predominantly poloidal axisymmetric field, as illustrated in Fig. C.2a. The first eigenvector is antisymmetric and there is evidence for a second, symmetric eigenvector, albeit the larger noise (Fig. C.2b). This suggests that CN Leo is the most axisymmetric M dwarf of our sample.

The per-epoch Stokes V mean profiles show a varying amplitude, together with the amplitude of the rotational modulation of the coefficients, which could be due to variability in either field strength or axisymmetry (see Fig. C.2c). The first epoch 2019a appears to be highly axisymmetric, as the coefficients cluster around zero. The second epoch 2019b2020a has a clear
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** PCA analysis for EV Lac. (a) Mean profile (red) for all observations and its decomposition in the antisymmetric (blue dashed) and symmetric (yellow dotted) component (with respect to the line centre) related to the poloidal and toroidal axisymmetric field, respectively. (b) First three eigenvectors of the mean-subtracted Stokes V profiles. (c) Mean profile (left column), and the coefficients of the first three eigenvectors (three columns to the right) for each season (one season per row). The mean profiles of the individual seasons are plotted in the same format as in panel a. The coefficients are colour-coded by rotation cycle.</figcaption></figure>

sinusoidal variation in the coefficients of the first eigenvector, implying a larger tilt angle of the dipole relative to the other epochs, and the mean of the coefficients is offset towards positive values. The third epoch 2020b2021a is characterised by a less clear sinusoidal variation in the coefficients, with their mean offset towards negative values. The fourth epoch 2021b2022a resemble the first one, with the mean of the coefficients around zero.

**6. Magnetic imaging**

We applied Zeeman-Doppler imaging to reconstruct the large-scale magnetic field at the surface of EV Lac, DS Leo and

CN Leo. The magnetic geometry is described as the sum of a poloidal and a toroidal component, which are both expressed through spherical harmonics decomposition. In particular, we employed the formalism outlined in Lehmann & Donati (2022). The algorithm proceeds iteratively, by synthesising Stokes V profiles and comparing them with the observations, in order to fit the spherical harmonics coefficients $\alpha_{\ell,m}$, $\beta_{\ell,m}$, and $\gamma_{\ell,m}$ (with $\ell$ and $m$ the degree and order of the mode, respectively), until a maximum-entropy solution at a fixed reduced $\chi^2$ is reached (Skilling & Bryan 1984; Semel 1989; Donati & Brown 1997). We employed the zdipy code described in Folsom et al. (2018), and which was
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** PCA analysis for EV Lac. (a) Mean profile (red) for all observations and its decomposition in the antisymmetric (blue dashed) and symmetric (yellow dotted) component (with respect to the line centre) related to the poloidal and toroidal axisymmetric field, respectively. (b) First three eigenvectors of the mean-subtracted Stokes V profiles. (c) Mean profile (left column), and the coefficients of the first three eigenvectors (three columns to the right) for each season (one season per row). The mean profiles of the individual seasons are plotted in the same format as in panel a. The coefficients are colour-coded by rotation cycle.</figcaption></figure>

sinusoidal variation in the coefficients of the first eigenvector, implying a larger tilt angle of the dipole relative to the other epochs, and the mean of the coefficients is offset towards positive values. The third epoch 2020b2021a is characterised by a less clear sinusoidal variation in the coefficients, with their mean offset towards negative values. The fourth epoch 2021b2022a resemble the first one, with the mean of the coefficients around zero.

**6. Magnetic imaging**

We applied Zeeman-Doppler imaging to reconstruct the large-scale magnetic field at the surface of EV Lac, DS Leo and

CN Leo. The magnetic geometry is described as the sum of a poloidal and a toroidal component, which are both expressed through spherical harmonics decomposition. In particular, we employed the formalism outlined in Lehmann & Donati (2022). The algorithm proceeds iteratively, by synthesising Stokes V profiles and comparing them with the observations, in order to fit the spherical harmonics coefficients $\alpha_{\ell,m}$, $\beta_{\ell,m}$, and $\gamma_{\ell,m}$ (with $\ell$ and $m$ the degree and order of the mode, respectively), until a maximum-entropy solution at a fixed reduced $\chi^2$ is reached (Skilling & Bryan 1984; Semel 1989; Donati & Brown 1997). We employed the zdipy code described in Folsom et al. (2018), and which was
implemented by Bellotti et al. (2023b) to incorporate the Unno-Rachkovsky’s solutions to polarised radiative transfer equations in a Milne-Eddington atmosphere (Unno 1956; Rachkovsky 1967; Landi Degl’Innocenti & Landolfi 2004) and the filling factor formalism adopted by Morin et al. (2008).

The filling factors $f_I$ and $f_V$ represent the fraction of the cell of the stellar surface grid covered by magnetic regions and magnetic regions producing net circular polarisation, respectively (Morin et al. 2008; Kochukhov 2021). With the inclusion of $f_V$, we assume that the polarisation signal comes from a multitude of magnetic spots whose local field strength is $B/f_V$ distributed following a certain large-scale structure such that the magnetic field modulus averaged over a grid cell is equal to $B$. In practice, using $f_V$ enables us to reproduce the amplitude (scaling with the magnetic field $B$) and the Zeeman splitting (scaling with $B/f_V$) observed for Stokes V LSD profiles. The values of filling factors are assumed constant throughout the stellar surface grid.

The Unno-Rachkovsky models of the local line profiles (for both Stokes I and V) are described by the following parameters (del Toro Iniesta 2003; Landi Degl’Innocenti & Landolfi 2004): the Gaussian width ($w_G$, related to thermal broadening), the Lorentzian width ($w_L$, related to pressure broadening), the ratio of the line to continuum absorption coefficients ($\eta_0$), and the slope of the source function in the Milne-Eddington atmosphere ($\beta$). Before applying ZDI, and for each epoch of our three stars, we perform a parameter optimisation based on a $\chi^2$ minimisation approach. In practice, we generate a series of synthetic Stokes I profiles for a grid of parameters, and compare them with the median observed Stokes I for a specific epoch, until a minimum $\chi^2$ is found. The parameters corresponding to the $\chi^2$ minimum are then used for the ZDI reconstruction. For EV Lac, we used $w_G = 0.1$ km s$^{-1}$, $w_L = 12.0$ km s$^{-1}$, and $\eta_0 = 9.8$, for DS Leo, we used $w_G = 0.5$ km s$^{-1}$, $w_L = 4.4$ km s$^{-1}$, and $\eta_0 = 17.0$, and for CN Leo, we used $w_G = 0.5$ km s$^{-1}$, $w_L = 10.0$ km s$^{-1}$, and $\eta_0 = 9.5$. The value of $\beta$ is derived considering that, for a Milne-Eddington atmosphere, the local continuum flux relative to the flux at disc centre is

$$I_c/I_c^0 = (1 + \beta \cos \theta)/(1 + \beta), \qquad (5)$$

where $\theta$ is the angle between the line of sight and stellar surface normal (Landi Degl’Innocenti & Landolfi 2004). A standard linear limb darkening law is given, for instance, by

$$I_c/I_c^0 = 1 - \eta + \eta \cos \theta, \qquad (6)$$

where $\eta$ is the limb darkening coefficient (Gray 2005). From Eqs. (5) and (6), we derive

$$\beta = \eta/(1 - \eta). \qquad (7)$$

In our case, we adopted a linear limb darkening coefficient in H band of 0.2 (Claret & Bloemen 2011), hence $\beta$ is consistently fixed to 0.25 in all ZDI reconstructions (for more details, see Erba et al. 2024).

Once the line parameters were fixed, we searched for the optimised value of filling factor $f_V$ for each epoch of our stars. The procedure is similar to an optimisation of stellar parameters or differential rotation search (see Petit et al. 2002). We run ZDI for a grid of $f_V$ values between 1.0% and 100%, and we recorded the $\chi^2$ reached by ZDI each time. The $\chi^2$ distribution is fit with a parabola around the minimum, and the value of $f_V$ corresponding to the minimum $\chi^2$ represents the optimal value. The $1\sigma$ error bars on $f_V$ are determined as the $\Delta\chi^2 = 1$ variation away from the $\chi^2$ minimum (Press et al. 1992).

<figure><img src="image_11.png" /><figcaption>Fig. 10. Reconstructed ZDI maps in flattened polar view of EV Lac for (from left to right) 2019b, 2020b2021a, and 2021b. In each column the radial (top), azimuthal (middle), and meridional (bottom) components of the magnetic field vector are displayed. The radial ticks are located at the rotational phases when the observations were collected, while the concentric circles represent different stellar latitudes: +30° and +60° (dashed lines), and equator (solid line). The colour bar range is set by the maximum (in absolute value) of the magnetic field and illustrates the positive (red) and negative (blue) magnetic polarity for each epoch.</figcaption></figure>

Finally, the maximum degree of the harmonic expansion was set to $l_{max} = 8$ for all stars, consistently with the spatial resolution dictated by $v_{eq} \sin i$ and previous ZDI reconstructions (Morin et al. 2008; Hébrard et al. 2016). We note that only the modes with $l$ up to three have the most significant contribution in the reconstructed field maps. We did not observe appreciable differences in the ZDI reconstructions for filling factor $f_I$ between 0.0 and up to 0.5 in some cases, hence we fixed it to $f_I = 0.0$.

We summarise the properties of all the magnetic field reconstructions in Table 4. The obliquity refers to the colatitude of the maximum of the dipolar component. It is obtained by first computing the poloidal-dipolar component of the field using the $\alpha_{lm}$ coefficient, and then by locating the colatitude associated with the maximum.

## 6.1. EV Lac

The optical maps of EV Lac where reconstructed by Morin et al. (2008) using the 2006 and 2007 data collected with ESPaDOnS and Narval. They found a strong ($B_{mean} = 500$ G), non-axisymmetric, mostly dipolar field, composed mainly by two magnetic spots of distinct polarity at opposite longitudes. They were able to constrain a differential rotation rate of 1.7 mrad d$^{-1}$ which was consistent with solid body rotation within $3\sigma$.

To carry out tomographic inversion for EV Lac near-infrared SPIRou data, we assumed $i = 60^\circ$, $v_{eq} \sin i = 4.0$ km s$^{-1}$, $P_{rot} = 4.36$ d, and solid body rotation (Morin et al. 2008). The Stokes V time series is shown in Fig. D.1. The profiles are fitted down to a $\chi^2$ level of 0.9, 1.2, and 1.4 from an initial value of 1.2, 5.0, and 6.5 for 2019b, 2020b2021a, and 2021b, respectively.

We constrained the filling factor $f_V$ to 9% and 19% for 2020b2021a and 2021b, while we could not for 2019b, hence we set it to 100%. The typical error bar on $f_V$ is at most 1%. Using either 9% or 19% on 2019b has only marginal effects on the map and the magnetic energy repartition, and the same level
implemented by Bellotti et al. (2023b) to incorporate the Unno-Rachkovsky's solutions to polarised radiative transfer equations in a Milne-Eddington atmosphere (Unno 1956; Rachkovsky 1967; Landi Degl'Innocenti & Landolfi 2004) and the filling factor formalism adopted by Morin et al. (2008).

The filling factors $f_I$ and $f_V$ represent the fraction of the cell of the stellar surface grid covered by magnetic regions and magnetic regions producing net circular polarisation, respectively (Morin et al. 2008; Kochukhov 2021). With the inclusion of $f_V$, we assume that the polarisation signal comes from a multitude of magnetic spots whose local field strength is $B/f_V$ distributed following a certain large-scale structure such that the magnetic field modulus averaged over a grid cell is equal to $B$. In practice, using $f_V$ enables us to reproduce the amplitude (scaling with the magnetic field $B$) and the Zeeman splitting (scaling with $B/f_V$) observed for Stokes V LSD profiles. The values of filling factors are assumed constant throughout the stellar surface grid.

The Unno-Rachkovsky models of the local line profiles (for both Stokes I and V) are described by the following parameters (del Toro Iniesta 2003; Landi Degl'Innocenti & Landolfi 2004): the Gaussian width ($w_G$, related to thermal broadening), the Lorentzian width ($w_L$, related to pressure broadening), the ratio of the line to continuum absorption coefficients ($\eta_0$), and the slope of the source function in the Milne-Eddington atmosphere ($\beta$). Before applying ZDI, and for each epoch of our three stars, we perform a parameter optimisation based on a $\chi^2$ minimisation approach. In practice, we generate a series of synthetic Stokes I profiles for a grid of parameters, and compare them with the median observed Stokes I for a specific epoch, until a minimum $\chi^2$ is found. The parameters corresponding to the $\chi^2$ minimum are then used for the ZDI reconstruction. For EV Lac, we used $w_G = 0.1$ km s$^{-1}$, $w_L = 12.0$ km s$^{-1}$, and $\eta_0 = 9.8$, for DS Leo, we used $w_G = 0.5$ km s$^{-1}$, $w_L = 4.4$ km s$^{-1}$, and $\eta_0 = 17.0$, and for CN Leo, we used $w_G = 0.5$ km s$^{-1}$, $w_L = 10.0$ km s$^{-1}$, and $\eta_0 = 9.5$. The value of $\beta$ is derived considering that, for a Milne-Eddington atmosphere, the local continuum flux relative to the flux at disc centre is

$$I_c/I_c^0 = (1 + \beta \cos \theta)/(1 + \beta), \qquad (5)$$

where $\theta$ is the angle between the line of sight and stellar surface normal (Landi Degl'Innocenti & Landolfi 2004). A standard linear limb darkening law is given, for instance, by

$$I_c/I_c^0 = 1 - \eta + \eta \cos \theta, \qquad (6)$$

where $\eta$ is the limb darkening coefficient (Gray 2005). From Eqs. (5) and (6), we derive

$$\beta = \eta/(1 - \eta). \qquad (7)$$

In our case, we adopted a linear limb darkening coefficient in H band of 0.2 (Claret & Bloemen 2011), hence $\beta$ is consistently fixed to 0.25 in all ZDI reconstructions (for more details, see Erba et al. 2024).

Once the line parameters were fixed, we searched for the optimised value of filling factor $f_V$ for each epoch of our stars. The procedure is similar to an optimisation of stellar parameters or differential rotation search (see Petit et al. 2002). We run ZDI for a grid of $f_V$ values between 1.0% and 100%, and we recorded the $\chi^2$ reached by ZDI each time. The $\chi^2$ distribution is fit with a parabola around the minimum, and the value of $f_V$ corresponding to the minimum $\chi^2$ represents the optimal value. The $1\sigma$ error bars on $f_V$ are determined as the $\Delta\chi^2 = 1$ variation away from the $\chi^2$ minimum (Press et al. 1992).

<figure><img src="image_11.png" /><figcaption>Fig. 10. Reconstructed ZDI maps in flattened polar view of EV Lac for (from left to right) 2019b, 2020b2021a, and 2021b. In each column the radial (top), azimuthal (middle), and meridional (bottom) components of the magnetic field vector are displayed. The radial ticks are located at the rotational phases when the observations were collected, while the concentric circles represent different stellar latitudes: +30° and +60° (dashed lines), and equator (solid line). The colour bar range is set by the maximum (in absolute value) of the magnetic field and illustrates the positive (red) and negative (blue) magnetic polarity for each epoch.</figcaption></figure>

Finally, the maximum degree of the harmonic expansion was set to $l_{max} = 8$ for all stars, consistently with the spatial resolution dictated by $v_{eq} \sin i$ and previous ZDI reconstructions (Morin et al. 2008; Hébrard et al. 2016). We note that only the modes with $l$ up to three have the most significant contribution in the reconstructed field maps. We did not observe appreciable differences in the ZDI reconstructions for filling factor $f_I$ between 0.0 and up to 0.5 in some cases, hence we fixed it to $f_I = 0.0$.

We summarise the properties of all the magnetic field reconstructions in Table 4. The obliquity refers to the colatitude of the maximum of the dipolar component. It is obtained by first computing the poloidal-dipolar component of the field using the $\alpha_{lm}$ coefficient, and then by locating the colatitude associated with the maximum.

## 6.1. EV Lac

The optical maps of EV Lac where reconstructed by Morin et al. (2008) using the 2006 and 2007 data collected with ESPaDOnS and Narval. They found a strong ($B_{mean} = 500$ G), non-axisymmetric, mostly dipolar field, composed mainly by two magnetic spots of distinct polarity at opposite longitudes. They were able to constrain a differential rotation rate of 1.7 mrad d$^{-1}$ which was consistent with solid body rotation within $3\sigma$.

To carry out tomographic inversion for EV Lac near-infrared SPIRou data, we assumed $i = 60°$, $v_{eq} \sin i = 4.0$ km s$^{-1}$, $P_{rot} = 4.36$ d, and solid body rotation (Morin et al. 2008). The Stokes V time series is shown in Fig. D.1. The profiles are fitted down to a $\chi^2$ level of 0.9, 1.2, and 1.4 from an initial value of 1.2, 5.0, and 6.5 for 2019b, 2020b2021a, and 2021b, respectively.

We constrained the filling factor $f_V$ to 9% and 19% for 2020b2021a and 2021b, while we could not for 2019b, hence we set it to 100%. The typical error bar on $f_V$ is at most 1%. Using either 9% or 19% on 2019b has only marginal effects on the map and the magnetic energy repartition, and the same level
of $\chi^2$ is achieved. This could result from a lower effective S/N at this epoch, which prevents us from reliably constraining $f_V$. The fact that we can fix $f_V$ = 100% for this epoch means that modelling the horizontal splitting of the Stokes V lobes is not necessary, and weak-field approximation could be equivalently assumed.

For the two remaining epochs, it is not straightforward to claim an evolution of the small-scale features represented by a change in $f_V$ from 9% to 19%. If we set $f_V$ to the median value of the two epochs (i.e., 14%), we observe only slight differences in the reconstructed maps: the field is more axisymmetric (an additional 4%) and in 2020b2021a the field is also more octupolar (additional 4%). While still limited, a variation of $f_V$ has more impact on 2020b2021a than 2021b because of the stronger constraint on $f_V$ in the former epoch.

The maps of the magnetic field and their characteristics are shown in Fig. 10 and Table 4. In all epochs, most of the magnetic energy is stored in the poloidal (>99%) and dipolar (>65%) components, with a substantial contribution from the quadrupolar modes (>10%). The field is non-axisymmetric, and the axisymmetric energy fraction decreases over time from 25% to 2%. Likewise, the obliquity of the field increases from 50° to 87°. The average magnetic field strength for the 2019b epoch is lower than the two other epochs. Considering that the Stokes V models of the 2019b epoch does not fully reproduce the amplitude of the LSD profiles, the magnetic field strength is likely underestimated. A differential rotation search on the SPIRou time series was inconclusive.

## 6.2. DS Leo

The large-scale magnetic field map for DS Leo was reconstructed initially by Donati et al. (2008) using Narval data collected in 2007 and 2008. They obtained a predominantly toroidal field geometry encircling the star, whereas the poloidal component was mainly dipolar. The average field strength was 100 G, and the axisymmetry of the poloidal component decreased between 2007 and 2008. Later, Hébrard et al. (2016) applied ZDI to HARPS-Pol and Narval observations collected in 2014, and consistently found a mostly-toroidal, axisymmetric geometry, with the poloidal component accounting for less than half the magnetic energy and mostly dipolar. Compared to the previous reconstructions, the average field strength decreased down to 80 G and the dipolar component increased from 50% to 88%, already suggesting a rapidly evolving magnetic field. There are two Narval epochs with unpublished data, namely 2010 and 2012, for which we recovered the magnetic field map in this work, as outlined in Appendix A.

We reconstructed the large-scale magnetic field for the 2020b2021a and 2021b2022a SPIRou epochs. The input parameters are $i = 60°$, $v_{eq} \sin i = 2 \text{ km s}^{-1}$ (Hébrard et al. 2016), $P_{\text{rot}} = 13.91 \text{ d}$, and we initially postulated solid body rotation. The time series of Stokes V profiles for the two epochs is shown in Appendix D. The profiles were fitted to a $\chi^2$ level of 1.4 and 1.2 from an initial value of 2.7 and 1.6 for 2020b2021a and 2021b2022a, respectively. We also constrained $f_V$ values of 10% and 6% for the two epochs (with error bars of 1%).

We searched for differential rotation using the method of Donati et al. (2000) and Petit et al. (2002). Basically, we fixed the entropy at a certain value and inspected a dense grid of ($P_{\text{rot,eq}}$, $d\Omega$) pairs to find the combination that minimises the $\chi^2$ between observations and synthetic models, as illustrated in Fig. 11. The best parameters are measured by fitting a 2D paraboloid to the $\chi^2$ distribution, and the error bars are obtained

<figure><img src="image_12.png" /><figcaption>Fig. 11. Differential rotation search for DS Leo. Top: 2020b2021a. Bottom: 2021b2022a. The plots show the $\chi^2$ landscape over a grid of ($P_{\text{rot,eq}}$, $d\Omega$) pairs, with the $1\sigma$ and $3\sigma$ contours. The best values are obtained by fitting a 2D paraboloid around the minimum, while their error bars are estimated from the projection of the $1\sigma$ contour on the respective axis (Press et al. 1992).</figcaption></figure>

from a variation of $\Delta\chi^2 = 1$ away from the minimum (Press et al. 1992; Petit et al. 2002). Differential rotation is implemented in zdipy with the following equation

$$ \Omega(\theta) = \Omega_{\text{eq}} - d\Omega \sin^2(\theta) \quad (8) $$

with $\Omega(\theta)$ the rotation frequency at colatitude $\theta$, $\Omega_{\text{eq}}$ the rotation frequency at equator, and $d\Omega$ the differential rotation rate (Folsom et al. 2018).

We found $P_{\text{rot,eq}} = 13.836 \pm 0.016 \text{ d}$ and $d\Omega = 0.0192 \pm 0.0020 \text{ rad d}^{-1}$ for 2020b2021a, implying a rotation period at the pole of $14.447 \pm 0.069 \text{ d}$ (see Eq. (8)). For 2021b2022a, we obtained $P_{\text{rot,eq}} = 13.469 \pm 0.054 \text{ d}$ and $d\Omega = 0.0448 \pm 0.0045 \text{ rad d}^{-1}$, and the rotation period at the pole was $14.899 \pm 0.172 \text{ d}$. Assuming the ($P_{\text{rot,eq}}$, $d\Omega$) pairs as input parameters for ZDI, the Stokes V fit is improved down to a $\chi^2$ level of 1.25 and 1.10 for 2020b2021a and 2021b2022a, respectively. Finally, the inverse of the differential rotation rate represents the pole-equator lap time ($t_{\text{lap}} = 2\pi/d\Omega$). For 2020b2021a (time span of 256 d) and 2021b2022a (time span of 205 d) respectively, we obtained a $t_{\text{lap}}$ of $327 \pm 34 \text{ d}$ and $140 \pm 14 \text{ d}$. During this time, magnetic regions on the surface are presumably distorted by the shear of differential rotation. The short-term variability emerging from the SPIRou time series of FWHM values (see Fig. 6) can be attributed to differential rotation since the equator-pole lap time is comparable to the time span of the epochs considered, meaning that differential rotation could have shifted the magnetic regions on the surface affecting the mean level of FWHM for that specific cycle.
of $\chi^2$ is achieved. This could result from a lower effective S/N at this epoch, which prevents us from reliably constraining $f_V$. The fact that we can fix $f_V = 100\%$ for this epoch means that modelling the horizontal splitting of the Stokes V lobes is not necessary, and weak-field approximation could be equivalently assumed.

For the two remaining epochs, it is not straightforward to claim an evolution of the small-scale features represented by a change in $f_V$ from 9% to 19%. If we set $f_V$ to the median value of the two epochs (i.e., 14%), we observe only slight differences in the reconstructed maps: the field is more axisymmetric (an additional 4%) and in 2020b2021a the field is also more octupolar (additional 4%). While still limited, a variation of $f_V$ has more impact on 2020b2021a than 2021b because of the stronger constraint on $f_V$ in the former epoch.

The maps of the magnetic field and their characteristics are shown in Fig. 10 and Table 4. In all epochs, most of the magnetic energy is stored in the poloidal (>99%) and dipolar (>65%) components, with a substantial contribution from the quadrupolar modes (>10%). The field is non-axisymmetric, and the axisymmetric energy fraction decreases over time from 25% to 2%. Likewise, the obliquity of the field increases from 50° to 87°. The average magnetic field strength for the 2019b epoch is lower than the two other epochs. Considering that the Stokes V models of the 2019b epoch does not fully reproduce the amplitude of the LSD profiles, the magnetic field strength is likely underestimated. A differential rotation search on the SPIRou time series was inconclusive.

## 6.2. DS Leo

The large-scale magnetic field map for DS Leo was reconstructed initially by Donati et al. (2008) using Narval data collected in 2007 and 2008. They obtained a predominantly toroidal field geometry encircling the star, whereas the poloidal component was mainly dipolar. The average field strength was 100 G, and the axisymmetry of the poloidal component decreased between 2007 and 2008. Later, Hébrard et al. (2016) applied ZDI to HARPS-Pol and Narval observations collected in 2014, and consistently found a mostly-toroidal, axisymmetric geometry, with the poloidal component accounting for less than half the magnetic energy and mostly dipolar. Compared to the previous reconstructions, the average field strength decreased down to 80 G and the dipolar component increased from 50% to 88%, already suggesting a rapidly evolving magnetic field. There are two Narval epochs with unpublished data, namely 2010 and 2012, for which we recovered the magnetic field map in this work, as outlined in Appendix A.

We reconstructed the large-scale magnetic field for the 2020b2021a and 2021b2022a SPIRou epochs. The input parameters are $i = 60°$, $v_{eq} \sin i = 2 \text{ km s}^{-1}$ (Hébrard et al. 2016), $P_{\text{rot}} = 13.91 \text{ d}$, and we initially postulated solid body rotation. The time series of Stokes V profiles for the two epochs is shown in Appendix D. The profiles were fitted to a $\chi^2$ level of 1.4 and 1.2 from an initial value of 2.7 and 1.6 for 2020b2021a and 2021b2022a, respectively. We also constrained $f_V$ values of 10% and 6% for the two epochs (with error bars of 1%).

We searched for differential rotation using the method of Donati et al. (2000) and Petit et al. (2002). Basically, we fixed the entropy at a certain value and inspected a dense grid of ($P_{\text{rot,eq}}$, $d\Omega$) pairs to find the combination that minimises the $\chi^2$ between observations and synthetic models, as illustrated in Fig. 11. The best parameters are measured by fitting a 2D paraboloid to the $\chi^2$ distribution, and the error bars are obtained

<figure><img src="image_12.png" /><figcaption>Fig. 11. Differential rotation search for DS Leo. Top: 2020b2021a. Bottom: 2021b2022a. The plots show the $\chi^2$ landscape over a grid of ($P_{\text{rot,eq}}$, $d\Omega$) pairs, with the $1\sigma$ and $3\sigma$ contours. The best values are obtained by fitting a 2D paraboloid around the minimum, while their error bars are estimated from the projection of the $1\sigma$ contour on the respective axis (Press et al. 1992).</figcaption></figure>

from a variation of $\Delta\chi^2 = 1$ away from the minimum (Press et al. 1992; Petit et al. 2002). Differential rotation is implemented in zdipy with the following equation

$$ \Omega(\theta) = \Omega_{\text{eq}} - d\Omega \sin^2(\theta) \quad (8) $$

with $\Omega(\theta)$ the rotation frequency at colatitude $\theta$, $\Omega_{\text{eq}}$ the rotation frequency at equator, and $d\Omega$ the differential rotation rate (Folsom et al. 2018).

We found $P_{\text{rot,eq}} = 13.836 \pm 0.016 \text{ d}$ and $d\Omega = 0.0192 \pm 0.0020 \text{ rad d}^{-1}$ for 2020b2021a, implying a rotation period at the pole of $14.447 \pm 0.069 \text{ d}$ (see Eq. (8)). For 2021b2022a, we obtained $P_{\text{rot,eq}} = 13.469 \pm 0.054 \text{ d}$ and $d\Omega = 0.0448 \pm 0.0045 \text{ rad d}^{-1}$, and the rotation period at the pole was $14.899 \pm 0.172 \text{ d}$. Assuming the ($P_{\text{rot,eq}}$, $d\Omega$) pairs as input parameters for ZDI, the Stokes V fit is improved down to a $\chi^2$ level of 1.25 and 1.10 for 2020b2021a and 2021b2022a, respectively. Finally, the inverse of the differential rotation rate represents the pole-equator lap time ($t_{\text{lap}} = 2\pi/d\Omega$). For 2020b2021a (time span of 256 d) and 2021b2022a (time span of 205 d) respectively, we obtained a $t_{\text{lap}}$ of $327 \pm 34 \text{ d}$ and $140 \pm 14 \text{ d}$. During this time, magnetic regions on the surface are presumably distorted by the shear of differential rotation. The short-term variability emerging from the SPIRou time series of FWHM values (see Fig. 6) can be attributed to differential rotation since the equator-pole lap time is comparable to the time span of the epochs considered, meaning that differential rotation could have shifted the magnetic regions on the surface affecting the mean level of FWHM for that specific cycle.
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Reconstructed ZDI maps in flattened polar view of DS Leo. Left: 2020b2021a. Right: 2021b2022a. The maps account for the constrained differential rotation. The format is the same as in Fig. 10.</figcaption></figure>

We obtained tighter constraints on $d\Omega$ for 2020b2021a relative to 2021b2022a because the field is less complex in the second epoch, that is to say it is characterised by a lower number of magnetic features tracking differential rotation. Besides, the error bars are smaller than the literature measurements because our time series contained a larger number of observations (Donati et al. 2008; Hébrard et al. 2016). In general, the error bars on the differential rotation parameters are only statistical and do not account for systematics, hence they are likely underestimated. The values of equator/pole period encompass our measurement obtained from the periodogram analysis of $B_l$, and fall within the range of literature values (Donati et al. 2008; Hébrard et al. 2016). Our $d\Omega$ values are generally lower than estimated in the literature, but compatible within $3\sigma$ from Donati et al. (2008), and $1\sigma$ from Hébrard et al. (2016). Such difference could be due to an evolution of the shear at the surface of the star, as studied also by Donati et al. (2023a) for AU Mic, but deciphering the exact mechanism is not a straightforward task. Figure 11 also shows that there is an anti-correlation between $P_{rot,eq}$ and $d\Omega$ (or equivalently an $\Omega_{eq}-d\Omega$ correlation), likely due to the fact that we mainly trace one latitude when searching for differential rotation.

The maps are shown in Fig. 12 and their properties are reported in Table 4. In 2020b2021a, the magnetic energy is distributed almost equally in poloidal and toroidal components, with the dipolar mode accounting for most of the energy (64%). While the toroidal component is mostly axisymmetric (90%), the poloidal component is largely non-axisymmetric (4%). In 2021b2022a, the poloidal component takes over the toroidal one, counting 73% of the magnetic energy against 26%. The dipolar mode remains the dominant one (83%), featuring a moderate increase similarly to the one reported in Hébrard et al. (2016). In terms of axisymmetry, the toroidal component is stable, while the poloidal one increases from 4% to 9%. Between the two epochs, the average field strength decreased from 44 to 18 G.

## 6.3. CN Leo

There is no reconstruction of the magnetic field in the literature for this star. From four spectropolarimetric observations, Morin et al. (2010) detected strong Zeeman signatures corresponding to longitudinal fields of 600 G (as large as the maximum of EV Lac). The shape of Stokes V profiles and the absence of intermittency in the amplitude of the profiles suggested an axisymmetric, poloidal dipole.

From the rotation period $P_{rot} = 2.70 \pm 0.01$ d we measured (see Sect. 3) and the equatorial projected velocity $v_{eq} \sin i =$

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Reconstructed ZDI maps in flattened polar view of CN Leo. From left to right: 2019a, 2019b2020a, 2020b2021a, and 2021b2022a. The format is the same as in Fig. 10.</figcaption></figure>

$2 \pm 1$ km s$^{-1}$ estimated by Cristofari et al. (2023), we obtained $R \sin i = P_{rot}v_{eq} \sin i/50.59 = 0.107 \pm 0.053 R_\odot$, where the denominator accounts for the unit conversion of the variables involved. Considering the radius measurement for CN Leo by Rabus et al. (2019) of $R = 0.151 \pm 0.005 R_\odot$, we infer an inclination of $i = 45 \pm 20^\circ$.

The reconstruction of the large-scale magnetic field at the surface of CN Leo was performed for the 2019a, 2019b2020a, 2020b2021a and 2021b2022a epochs. The input parameters are $i = 45^\circ$, $v_{eq} \sin i = 2.0$ km s$^{-1}$, $P_{rot} = 2.7$ d, and solid body rotation.

The profiles were fitted to a $\chi^2$ level of 0.9, 1.0, 1.2, and 1.2 from an initial value of 2.9, 3.5, 6.8, and 6.7 for 2019a, 2019b2020a, 2020b2021a, and 2021b2022a, respectively. The $f_V$ values were constrained to be 16%, 13%, 19% and 16% for the same epochs, with an error bar of 1%. The Stokes V time series is shown in Appendix D.

The maps of the magnetic field and their characteristics are shown in Fig. 13 and Table 4. In all epochs, most of the magnetic energy is stored in the poloidal (>99%), dipolar (>99%) and axisymmetric (>99%) components, without any significant variation over time. The field strength follows the sine-like long-term oscillation exhibited by $B_l$, meaning that it starts from an average field of 350 G, it decreases to 320 G, then increases to 490 G and finally decreases to 420 G.

## 6.4. Summary of the ZDI reconstructions

The evolution of the large-scale magnetic field for EV Lac, DS Leo, and CN Leo is summarised in Fig. 14. For completeness, we also included the evolution of AD Leo as reconstructed by Bellotti et al. (2023b), which sees a decrease of field strength accompanied by an increase of the obliquity of the magnetic field axis.

For EV Lac, the magnetic topology was found to be predominantly non-axisymmetric using optical data
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Reconstructed ZDI maps in flattened polar view of DS Leo. Left: 2020b2021a. Right: 2021b2022a. The maps account for the constrained differential rotation. The format is the same as in Fig. 10.</figcaption></figure>

We obtained tighter constraints on $d\Omega$ for 2020b2021a relative to 2021b2022a because the field is less complex in the second epoch, that is to say it is characterised by a lower number of magnetic features tracking differential rotation. Besides, the error bars are smaller than the literature measurements because our time series contained a larger number of observations (Donati et al. 2008; Hébrard et al. 2016). In general, the error bars on the differential rotation parameters are only statistical and do not account for systematics, hence they are likely underestimated. The values of equator/pole period encompass our measurement obtained from the periodogram analysis of $B_l$, and fall within the range of literature values (Donati et al. 2008; Hébrard et al. 2016). Our $d\Omega$ values are generally lower than estimated in the literature, but compatible within $3\sigma$ from Donati et al. (2008), and $1\sigma$ from Hébrard et al. (2016). Such difference could be due to an evolution of the shear at the surface of the star, as studied also by Donati et al. (2023a) for AU Mic, but deciphering the exact mechanism is not a straightforward task. Figure 11 also shows that there is an anti-correlation between $P_{rot,eq}$ and $d\Omega$ (or equivalently an $\Omega_{eq}-d\Omega$ correlation), likely due to the fact that we mainly trace one latitude when searching for differential rotation.

The maps are shown in Fig. 12 and their properties are reported in Table 4. In 2020b2021a, the magnetic energy is distributed almost equally in poloidal and toroidal components, with the dipolar mode accounting for most of the energy (64%). While the toroidal component is mostly axisymmetric (90%), the poloidal component is largely non-axisymmetric (4%). In 2021b2022a, the poloidal component takes over the toroidal one, counting 73% of the magnetic energy against 26%. The dipolar mode remains the dominant one (83%), featuring a moderate increase similarly to the one reported in Hébrard et al. (2016). In terms of axisymmetry, the toroidal component is stable, while the poloidal one increases from 4% to 9%. Between the two epochs, the average field strength decreased from 44 to 18 G.

## 6.3. CN Leo

There is no reconstruction of the magnetic field in the literature for this star. From four spectropolarimetric observations, Morin et al. (2010) detected strong Zeeman signatures corresponding to longitudinal fields of 600 G (as large as the maximum of EV Lac). The shape of Stokes V profiles and the absence of intermittency in the amplitude of the profiles suggested an axisymmetric, poloidal dipole.

From the rotation period $P_{rot} = 2.70 \pm 0.01$ d we measured (see Sect. 3) and the equatorial projected velocity $v_{eq} \sin i =$

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Reconstructed ZDI maps in flattened polar view of CN Leo. From left to right: 2019a, 2019b2020a, 2020b2021a, and 2021b2022a. The format is the same as in Fig. 10.</figcaption></figure>

$2 \pm 1$ km s$^{-1}$ estimated by Cristofari et al. (2023), we obtained $R \sin i = P_{rot}v_{eq} \sin i/50.59 = 0.107 \pm 0.053 R_\odot$, where the denominator accounts for the unit conversion of the variables involved. Considering the radius measurement for CN Leo by Rabus et al. (2019) of $R = 0.151 \pm 0.005 R_\odot$, we infer an inclination of $i = 45 \pm 20^\circ$.

The reconstruction of the large-scale magnetic field at the surface of CN Leo was performed for the 2019a, 2019b2020a, 2020b2021a and 2021b2022a epochs. The input parameters are $i = 45^\circ$, $v_{eq} \sin i = 2.0$ km s$^{-1}$, $P_{rot} = 2.7$ d, and solid body rotation.

The profiles were fitted to a $\chi^2$ level of 0.9, 1.0, 1.2, and 1.2 from an initial value of 2.9, 3.5, 6.8, and 6.7 for 2019a, 2019b2020a, 2020b2021a, and 2021b2022a, respectively. The $f_V$ values were constrained to be 16%, 13%, 19% and 16% for the same epochs, with an error bar of 1%. The Stokes V time series is shown in Appendix D.

The maps of the magnetic field and their characteristics are shown in Fig. 13 and Table 4. In all epochs, most of the magnetic energy is stored in the poloidal (>99%), dipolar (>99%) and axisymmetric (>99%) components, without any significant variation over time. The field strength follows the sine-like long-term oscillation exhibited by $B_l$, meaning that it starts from an average field of 350 G, it decreases to 320 G, then increases to 490 G and finally decreases to 420 G.

## 6.4. Summary of the ZDI reconstructions

The evolution of the large-scale magnetic field for EV Lac, DS Leo, and CN Leo is summarised in Fig. 14. For completeness, we also included the evolution of AD Leo as reconstructed by Bellotti et al. (2023b), which sees a decrease of field strength accompanied by an increase of the obliquity of the magnetic field axis.

For EV Lac, the magnetic topology was found to be predominantly non-axisymmetric using optical data
<table><caption>**Table 4.** Properties of the magnetic map for EV Lac, DS Leo, and CN Leo, reconstructed from different epochs of SPIRou near-infrared time series.</caption><thead><tr><th rowspan="2"></th><th colspan="3">EV Lac</th><th colspan="2">DS Leo</th><th colspan="3">CN Leo</th></tr><tr><th>2019b</th><th>2020b2021a</th><th>2021b</th><th>2020b2021a</th><th>2021b2022a</th><th>2019a</th><th>2019b2020a</th><th>2020b2021a</th><th>2021b2022a</th></tr></thead><tbody><tr><td>f<sub>V</sub> [%]</td><td>...</td><td>9</td><td>19</td><td>10</td><td>6</td><td>15</td><td>12</td><td>18</td><td>15</td></tr><tr><td>B<sub>mean</sub> [G]</td><td>72.5</td><td>297.7</td><td>265.1</td><td>44.3</td><td>17.9</td><td>350.2</td><td>321.8</td><td>486.9</td><td>420.1</td></tr><tr><td>B<sub>max</sub> [G]</td><td>175.8</td><td>700.0</td><td>690.5</td><td>115.2</td><td>34.9</td><td>650.4</td><td>607.8</td><td>937.6</td><td>820.9</td></tr><tr><td>B<sub>pol</sub> [%]</td><td>99.1</td><td>99.7</td><td>99.2</td><td>55.9</td><td>72.9</td><td>99.9</td><td>99.9</td><td>99.9</td><td>99.9</td></tr><tr><td>B<sub>tor</sub> [%]</td><td>0.9</td><td>0.3</td><td>0.8</td><td>44.1</td><td>27.0</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>B<sub>dip</sub> [%]</td><td>63.1</td><td>88.7</td><td>82.8</td><td>64.0</td><td>83.2</td><td>99.4</td><td>99.4</td><td>99.3</td><td>99.2</td></tr><tr><td>B<sub>quad</sub> [%]</td><td>25.4</td><td>9.4</td><td>14.6</td><td>28.5</td><td>13.1</td><td>0.5</td><td>0.6</td><td>0.6</td><td>0.7</td></tr><tr><td>B<sub>oct</sub> [%]</td><td>10.2</td><td>1.3</td><td>2.4</td><td>3.9</td><td>2.6</td><td>0.0</td><td>0.0</td><td>0.1</td><td>0.1</td></tr><tr><td>B<sub>axisym</sub> [%]</td><td>25.1</td><td>12.1</td><td>1.8</td><td>39.4</td><td>30.6</td><td>99.5</td><td>99.6</td><td>99.9</td><td>99.9</td></tr><tr><td>B<sub>axisym.pol</sub> [%]</td><td>24.6</td><td>12.0</td><td>1.2</td><td>3.9</td><td>9.2</td><td>99.5</td><td>99.6</td><td>99.9</td><td>99.9</td></tr><tr><td>Obliquity [°]</td><td>50.5</td><td>68.5</td><td>86.5</td><td>75.5</td><td>69.5</td><td>4.5</td><td>3.5</td><td>1.5</td><td>1.5</td></tr></tbody></table>

**Notes.** The following quantities are listed: filling factor on Stokes V; mean magnetic strength; maximum magnetic strength; poloidal and toroidal magnetic energy as a fraction of the total energy; dipolar, quadrupolar and octupolar magnetic energy as a fraction of the poloidal energy; axisymmetric magnetic energy as a fraction of the total energy; and obliquity of the dipolar component relative to the rotation axis. To compute the local magnetic field (i.e., within a grid cell of the ZDI stellar model) the field strength should be divided by the associated filling factor $f_V$.

(Morin et al. 2008), and we observed an even less axisymmetric field from SPIRou data, with the negative pole of the dipolar component lying approximately at equator. From the SPIRou observations, we also noted a rise in field strength from 57 G in 2019b to around 200 G in 2020b2021a and 2021b, following a similar evolution with respect to $B_t$.

For DS Leo, the poloidal component saw a progressive increase with respect to the toroidal one (from 20 to 70% of the total magnetic energy), and it remained highly non-axisymmetric throughout (Donati et al. 2008; Hébrard et al. 2016). The complexity of the poloidal component of the field has also decreased in the same timescale, since the field is predominantly dipolar. Moreover, the average reconstructed field strength has reduced, from an initial 100 G in 2007 to 15 G in the 2021b2022a SPIRou epoch. The panel of DS Leo in Fig. 14 includes two field topologies, in 2010 and 2012, which we reconstructed using unpublished Narval data (see Appendix A).

For CN Leo, there is no substantial evolution in four years, the topology being mainly poloidal, dipolar and axisymmetric. The only feature is a varying field strength, whose oscillations correlate with the longitudinal field evolution (in absolute value).

Similarly to what reported by Bellotti et al. (2023b), we caution that the field strength reported may be underestimated, owing to the limitation of the ZDI model at reproducing the shape of the Stokes V lobes (see e.g., Fig. D.1). The near-infrared Stokes V profiles manifest evident noise that is not rotationally modulated, and which prevents the ZDI model to capture all the information present in the profile.

To describe the magnetic field vector, we employed the formalism of Lehmann & Donati (2022), which sees the substitution of $\beta_{lm}$ with $\alpha_{lm} + \beta_{lm}$ in the spherical harmonics equations. The main effect of using the latter is to change the appearance of the meridional and azimuthal magnetic field maps, making them resemble the radial field map, and avoid unnecessary large values of $\beta_{lm}$. If we use the standard formalism (see e.g., Donati et al. 2006b) with $\alpha_{lm}$ and $\beta_{lm}$ disjointed, the reconstructed maps of EV Lac vary less than 3% in poloidal energy fraction, less than 5% in dipolar, quadrupolar, and octupolar components, and less than 2% in axisymmetric fraction. The average field strength decreases by at most 100 G. For CN Leo, the variations are similar or lower. This is in agreement with what was previously

found also for AD Leo (see appendix of Bellotti 2023). Overall, the properties of the magnetic topology and its long-term evolution are reconstructed consistently adopting either of the two spherical harmonics formalisms.

## 7. Discussion and conclusions

In this work we presented the long-term spectropolarimetric monitoring of three well-studied, active M dwarfs: EV Lac, DS Leo, and CN Leo. We used archival optical data collected with ESPaDOnS and Narval, as well as near-infrared SPIRou data obtained as part of the Legacy Survey between 2019 and 2022. We carried out distinct analyses to capture the evolution of the magnetic field, employing the longitudinal field as a proxy of the large-scale field component and the FWHM of Stokes I LSD profiles as a proxy for the small-scale component. We analysed qualitatively the secular evolution of the large-scale magnetic field via principal component analysis and reconstructed its topology by means of Zeeman-Doppler imaging.

We found different trends in the magnetic field behaviour, potentially hinting at a variety of magnetic cycles. Our conclusions are the following:

1. The longitudinal magnetic field analysis showed a pulsating trend for EV Lac with variations between ±400 G around 2007, then between ±200 G in 2016, and finally between ±300 G in 2020. However, the dearth of observations between 2010 and 2019 prevent us from drawing firm conclusions on this behaviour. DS Leo did not manifest any specific trend as the mean field remained reasonably stable around 0 G, and the measurements within ±50 G. CN Leo exhibited a sine-like trend, characterised by a period of 2.7 yr, an amplitude of 100 G, and an average of −500 G.

2. We observed rotational modulation of the FWHM of Stokes I for optical observations of EV Lac and DS Leo, and near-infrared observations of DS Leo, likely stemming from distinct contributions of the Zeeman effect in the two wavelength domains. The FWHM of near-infrared line profiles for stars with intense magnetic fields is expected to be more impacted than at optical wavelengths, resulting in larger scatter in the time series.
<table><caption>**Table 4.** Properties of the magnetic map for EV Lac, DS Leo, and CN Leo, reconstructed from different epochs of SPIRou near-infrared time series.</caption><thead><tr><th rowspan="2"></th><th colspan="3">EV Lac</th><th colspan="2">DS Leo</th><th colspan="3">CN Leo</th></tr><tr><th>2019b</th><th>2020b2021a</th><th>2021b</th><th>2020b2021a</th><th>2021b2022a</th><th>2019a</th><th>2019b2020a</th><th>2020b2021a</th><th>2021b2022a</th></tr></thead><tbody><tr><td>f<sub>V</sub> [%]</td><td>...</td><td>9</td><td>19</td><td>10</td><td>6</td><td>15</td><td>12</td><td>18</td><td>15</td></tr><tr><td>B<sub>mean</sub> [G]</td><td>72.5</td><td>297.7</td><td>265.1</td><td>44.3</td><td>17.9</td><td>350.2</td><td>321.8</td><td>486.9</td><td>420.1</td></tr><tr><td>B<sub>max</sub> [G]</td><td>175.8</td><td>700.0</td><td>690.5</td><td>115.2</td><td>34.9</td><td>650.4</td><td>607.8</td><td>937.6</td><td>820.9</td></tr><tr><td>B<sub>pol</sub> [%]</td><td>99.1</td><td>99.7</td><td>99.2</td><td>55.9</td><td>72.9</td><td>99.9</td><td>99.9</td><td>99.9</td><td>99.9</td></tr><tr><td>B<sub>tor</sub> [%]</td><td>0.9</td><td>0.3</td><td>0.8</td><td>44.1</td><td>27.0</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>B<sub>dip</sub> [%]</td><td>63.1</td><td>88.7</td><td>82.8</td><td>64.0</td><td>83.2</td><td>99.4</td><td>99.4</td><td>99.3</td><td>99.2</td></tr><tr><td>B<sub>quad</sub> [%]</td><td>25.4</td><td>9.4</td><td>14.6</td><td>28.5</td><td>13.1</td><td>0.5</td><td>0.6</td><td>0.6</td><td>0.7</td></tr><tr><td>B<sub>oct</sub> [%]</td><td>10.2</td><td>1.3</td><td>2.4</td><td>3.9</td><td>2.6</td><td>0.0</td><td>0.0</td><td>0.1</td><td>0.1</td></tr><tr><td>B<sub>axisym</sub> [%]</td><td>25.1</td><td>12.1</td><td>1.8</td><td>39.4</td><td>30.6</td><td>99.5</td><td>99.6</td><td>99.9</td><td>99.9</td></tr><tr><td>B<sub>axisym.pol</sub> [%]</td><td>24.6</td><td>12.0</td><td>1.2</td><td>3.9</td><td>9.2</td><td>99.5</td><td>99.6</td><td>99.9</td><td>99.9</td></tr><tr><td>Obliquity [°]</td><td>50.5</td><td>68.5</td><td>86.5</td><td>75.5</td><td>69.5</td><td>4.5</td><td>3.5</td><td>1.5</td><td>1.5</td></tr></tbody></table>

**Notes.** The following quantities are listed: filling factor on Stokes V; mean magnetic strength; maximum magnetic strength; poloidal and toroidal magnetic energy as a fraction of the total energy; dipolar, quadrupolar and octupolar magnetic energy as a fraction of the poloidal energy; axisymmetric magnetic energy as a fraction of the total energy; and obliquity of the dipolar component relative to the rotation axis. To compute the local magnetic field (i.e., within a grid cell of the ZDI stellar model) the field strength should be divided by the associated filling factor $f_V$.

(Morin et al. 2008), and we observed an even less axisymmetric field from SPIRou data, with the negative pole of the dipolar component lying approximately at equator. From the SPIRou observations, we also noted a rise in field strength from 57 G in 2019b to around 200 G in 2020b2021a and 2021b, following a similar evolution with respect to $B_t$.

For DS Leo, the poloidal component saw a progressive increase with respect to the toroidal one (from 20 to 70% of the total magnetic energy), and it remained highly non-axisymmetric throughout (Donati et al. 2008; Hébrard et al. 2016). The complexity of the poloidal component of the field has also decreased in the same timescale, since the field is predominantly dipolar. Moreover, the average reconstructed field strength has reduced, from an initial 100 G in 2007 to 15 G in the 2021b2022a SPIRou epoch. The panel of DS Leo in Fig. 14 includes two field topologies, in 2010 and 2012, which we reconstructed using unpublished Narval data (see Appendix A).

For CN Leo, there is no substantial evolution in four years, the topology being mainly poloidal, dipolar and axisymmetric. The only feature is a varying field strength, whose oscillations correlate with the longitudinal field evolution (in absolute value).

Similarly to what reported by Bellotti et al. (2023b), we caution that the field strength reported may be underestimated, owing to the limitation of the ZDI model at reproducing the shape of the Stokes V lobes (see e.g., Fig. D.1). The near-infrared Stokes V profiles manifest evident noise that is not rotationally modulated, and which prevents the ZDI model to capture all the information present in the profile.

To describe the magnetic field vector, we employed the formalism of Lehmann & Donati (2022), which sees the substitution of $\beta_{lm}$ with $\alpha_{lm} + \beta_{lm}$ in the spherical harmonics equations. The main effect of using the latter is to change the appearance of the meridional and azimuthal magnetic field maps, making them resemble the radial field map, and avoid unnecessary large values of $\beta_{lm}$. If we use the standard formalism (see e.g., Donati et al. 2006b) with $\alpha_{lm}$ and $\beta_{lm}$ disjointed, the reconstructed maps of EV Lac vary less than 3% in poloidal energy fraction, less than 5% in dipolar, quadrupolar, and octupolar components, and less than 2% in axisymmetric fraction. The average field strength decreases by at most 100 G. For CN Leo, the variations are similar or lower. This is in agreement with what was previously

found also for AD Leo (see appendix of Bellotti 2023). Overall, the properties of the magnetic topology and its long-term evolution are reconstructed consistently adopting either of the two spherical harmonics formalisms.

## 7. Discussion and conclusions

In this work we presented the long-term spectropolarimetric monitoring of three well-studied, active M dwarfs: EV Lac, DS Leo, and CN Leo. We used archival optical data collected with ESPaDOnS and Narval, as well as near-infrared SPIRou data obtained as part of the Legacy Survey between 2019 and 2022. We carried out distinct analyses to capture the evolution of the magnetic field, employing the longitudinal field as a proxy of the large-scale field component and the FWHM of Stokes I LSD profiles as a proxy for the small-scale component. We analysed qualitatively the secular evolution of the large-scale magnetic field via principal component analysis and reconstructed its topology by means of Zeeman-Doppler imaging.

We found different trends in the magnetic field behaviour, potentially hinting at a variety of magnetic cycles. Our conclusions are the following:

1. The longitudinal magnetic field analysis showed a pulsating trend for EV Lac with variations between ±400 G around 2007, then between ±200 G in 2016, and finally between ±300 G in 2020. However, the dearth of observations between 2010 and 2019 prevent us from drawing firm conclusions on this behaviour. DS Leo did not manifest any specific trend as the mean field remained reasonably stable around 0 G, and the measurements within ±50 G. CN Leo exhibited a sine-like trend, characterised by a period of 2.7 yr, an amplitude of 100 G, and an average of −500 G.

2. We observed rotational modulation of the FWHM of Stokes I for optical observations of EV Lac and DS Leo, and near-infrared observations of DS Leo, likely stemming from distinct contributions of the Zeeman effect in the two wavelength domains. The FWHM of near-infrared line profiles for stars with intense magnetic fields is expected to be more impacted than at optical wavelengths, resulting in larger scatter in the time series.
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Evolution of the magnetic topology of AD Leo, EV Lac, CN Leo, and DS Leo. The temporal baseline covers 15 yr, including the archival optical maps (Donati et al. 2008; Morin et al. 2008; Hébrard et al. 2016) and the near-infrared maps reconstructed with SPIRou data in this work and in Bellotti et al. (2023b) for AD Leo. The field strength, level of axisymmetry, and dominant field component (poloidal or toroidal) are encoded with the symbol size, shape, and colour, respectively. The red and blue data points represent poloidal- and toroidal-dominated field geometries, whereas round and star-like shapes indicate mostly axisymmetric and non-axisymmetric configurations.</figcaption></figure>

3. For the optical epochs of EV Lac and DS Leo where the rotational modulation of the FWHM is evident, we observed both negative and positive correlations with $|B_l|$, whereas in some cases this correlation was absent. The interpretation of these results is not straightforward. They indicate an underlying complexity in the link between the small-scale field (traced by FWHM) and the large-scale field (traced by $B_l$), which may not be encapsulated in a simple scaling.

4. The long-term behaviour of the epoch-averaged FWHM of Stokes *I* reflects the secular evolution of the longitudinal field. This is enhanced in particular when high-$g_{\text{eff}}$ lines are used, whereas low-$g_{\text{eff}}$ lines tend to yield stable FWHM values throughout the time series.

5. Using PCA, we obtained the variations in axisymmetry, poloidal-to-toroidal fraction of the axisymmetric component, and complexity for EV Lac, DS Leo, and CN Leo directly from the LSD Stokes *V* time series, without prior assumptions on stellar parameters or the need of elaborate field topology models.

6. The ZDI reconstructions also revealed distinct secular changes in the magnetic topology. For EV Lac the field remained mainly poloidal, dipolar, and with a significant contribution of quadrupolar field. The axisymmetry has decreased almost to 0% in the most recent epoch. For

DS Leo, the poloidal component increased over time with a simultaneous decrease in axisymmetry. The toroidal component has a substantial contribution to the energy budget still, and is mainly axisymmetric, similarly to the optical reconstructions (Donati et al. 2008; Hébrard et al. 2016). For CN Leo, the field maintained a poloidal, dipolar, and axisymmetric configuration, with the field strength following a fluctuating trend similar to $B_l$.

If we group the long-term evolution of the large-scale field geometry according to spectral type, it seems that the cyclic variability decreases towards later spectral types, as also pointed out by Fuhrmeister et al. (2023) from a long-term monitoring of chromospheric activity indicators. From our work we see that the field of DS Leo (early M type) undergoes rapid variations on short timescales; the field of EV Lac still manifests an evident topological change, but on longer timescales (similar to AD Leo, Bellotti et al. 2023b); and for CN Leo (late M type), there is no substantial evolution. This phenomenon could be correlated to the magnetic field intensity, with stronger fields (i.e. later spectral types) being able to quench surface shear and induce intrinsic stability, as already pointed out by Donati et al. (2008) and Morin et al. (2008).

With the exception of DS Leo, activity cycles were reported for EV Lac and CN Leo from photometric or chromospheric activity monitoring. For EV Lac, Mavridis & Avgoloupis (1986) constrained a 5 yr activity cycle based on flare activity and B-band photometric oscillations over a baseline of 9 yr. Given the gap in our $B_l$ time series between 2010 and 2016, we cannot easily compare the reported timescale with our $B_l$ variations. For CN Leo, photometric analysis of the All Sky Automated Survey (ASAS) light curves revealed an activity cycle of 8.9 ± 0.2 yr (Suárez Mascareño et al. 2016) or 12.0 ± 4.4 yr (Irving et al. 2023). These are 3.3 and 4.4 times longer than the 2.7 yr sine-like variations we recorded for $B_l$. In contrast, Fuhrmeister et al. (2023) and Mignon et al. (2023) did not report any significant long-term periodicity from chromospheric activity monitoring of CN Leo. These result may resemble what is described for Sun-like stars because there are cases in which magnetic cycles and polarity flips are phased with chromospheric cycles (Boro Saikia et al. 2016; Jeffers et al. 2017), whereas in other cases regular chromospheric oscillations are not reflected in polarimetric variations (Boro Saikia et al. 2022). At the same time, different activity indexes may be sensitive to different regions in the stellar atmospheres, and therefore timescales, which may explain why they did not capture long-term variability for CN Leo (e.g., Cortés-Zuleta et al. 2023; Mignon et al. 2023).

In conclusion, the peculiar evolution of the large-scale magnetic field for the stars examined here provides us with a glimpse of a potential variety in cyclic trends for M dwarfs. This advocates for additional spectropolarimetric monitoring of active M dwarfs for a better constraint on the timescales involved and the extent of such variety over distinct stellar masses and rotation periods.

**Acknowledgements.** We thank the anonymous referee for the thorough review of the manuscript and the suggestions. We acknowledge funding from the French National Research Agency (ANR) under contract number ANR-18-CE31-0019 (SPlaSH). SB acknowledges funding from the European Space Agency (ESA), under the visiting researcher programme. XD acknowledges funding within the framework of the Investissements d'Avenir programme (ANR-15-IDEX-02), through the funding of the 'Origin of Life' project of the Univ. Grenoble-Alpes. EM acknowledges funding from FAPEMIG under project number APQ-02493-22 and research productivity grant number 309829/2022-4 awarded by the CNPq, Brazil. BK acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Evolution of the magnetic topology of AD Leo, EV Lac, CN Leo, and DS Leo. The temporal baseline covers 15 yr, including the archival optical maps (Donati et al. 2008; Morin et al. 2008; Hébrard et al. 2016) and the near-infrared maps reconstructed with SPIRou data in this work and in Bellotti et al. (2023b) for AD Leo. The field strength, level of axisymmetry, and dominant field component (poloidal or toroidal) are encoded with the symbol size, shape, and colour, respectively. The red and blue data points represent poloidal- and toroidal-dominated field geometries, whereas round and star-like shapes indicate mostly axisymmetric and non-axisymmetric configurations.</figcaption></figure>

3. For the optical epochs of EV Lac and DS Leo where the rotational modulation of the FWHM is evident, we observed both negative and positive correlations with $|B_l|$, whereas in some cases this correlation was absent. The interpretation of these results is not straightforward. They indicate an underlying complexity in the link between the small-scale field (traced by FWHM) and the large-scale field (traced by $B_l$), which may not be encapsulated in a simple scaling.

4. The long-term behaviour of the epoch-averaged FWHM of Stokes *I* reflects the secular evolution of the longitudinal field. This is enhanced in particular when high-$g_{\text{eff}}$ lines are used, whereas low-$g_{\text{eff}}$ lines tend to yield stable FWHM values throughout the time series.

5. Using PCA, we obtained the variations in axisymmetry, poloidal-to-toroidal fraction of the axisymmetric component, and complexity for EV Lac, DS Leo, and CN Leo directly from the LSD Stokes *V* time series, without prior assumptions on stellar parameters or the need of elaborate field topology models.

6. The ZDI reconstructions also revealed distinct secular changes in the magnetic topology. For EV Lac the field remained mainly poloidal, dipolar, and with a significant contribution of quadrupolar field. The axisymmetry has decreased almost to 0% in the most recent epoch. For

DS Leo, the poloidal component increased over time with a simultaneous decrease in axisymmetry. The toroidal component has a substantial contribution to the energy budget still, and is mainly axisymmetric, similarly to the optical reconstructions (Donati et al. 2008; Hébrard et al. 2016). For CN Leo, the field maintained a poloidal, dipolar, and axisymmetric configuration, with the field strength following a fluctuating trend similar to $B_l$.

If we group the long-term evolution of the large-scale field geometry according to spectral type, it seems that the cyclic variability decreases towards later spectral types, as also pointed out by Fuhrmeister et al. (2023) from a long-term monitoring of chromospheric activity indicators. From our work we see that the field of DS Leo (early M type) undergoes rapid variations on short timescales; the field of EV Lac still manifests an evident topological change, but on longer timescales (similar to AD Leo, Bellotti et al. 2023b); and for CN Leo (late M type), there is no substantial evolution. This phenomenon could be correlated to the magnetic field intensity, with stronger fields (i.e. later spectral types) being able to quench surface shear and induce intrinsic stability, as already pointed out by Donati et al. (2008) and Morin et al. (2008).

With the exception of DS Leo, activity cycles were reported for EV Lac and CN Leo from photometric or chromospheric activity monitoring. For EV Lac, Mavridis & Avgoloupis (1986) constrained a 5 yr activity cycle based on flare activity and B-band photometric oscillations over a baseline of 9 yr. Given the gap in our $B_l$ time series between 2010 and 2016, we cannot easily compare the reported timescale with our $B_l$ variations. For CN Leo, photometric analysis of the All Sky Automated Survey (ASAS) light curves revealed an activity cycle of 8.9 ± 0.2 yr (Suárez Mascareño et al. 2016) or 12.0 ± 4.4 yr (Irving et al. 2023). These are 3.3 and 4.4 times longer than the 2.7 yr sine-like variations we recorded for $B_l$. In contrast, Fuhrmeister et al. (2023) and Mignon et al. (2023) did not report any significant long-term periodicity from chromospheric activity monitoring of CN Leo. These result may resemble what is described for Sun-like stars because there are cases in which magnetic cycles and polarity flips are phased with chromospheric cycles (Boro Saikia et al. 2016; Jeffers et al. 2017), whereas in other cases regular chromospheric oscillations are not reflected in polarimetric variations (Boro Saikia et al. 2022). At the same time, different activity indexes may be sensitive to different regions in the stellar atmospheres, and therefore timescales, which may explain why they did not capture long-term variability for CN Leo (e.g., Cortés-Zuleta et al. 2023; Mignon et al. 2023).

In conclusion, the peculiar evolution of the large-scale magnetic field for the stars examined here provides us with a glimpse of a potential variety in cyclic trends for M dwarfs. This advocates for additional spectropolarimetric monitoring of active M dwarfs for a better constraint on the timescales involved and the extent of such variety over distinct stellar masses and rotation periods.

*Acknowledgements.* We thank the anonymous referee for the thorough review of the manuscript and the suggestions. We acknowledge funding from the French National Research Agency (ANR) under contract number ANR-18-CE31-0019 (SPlaSH). SB acknowledges funding from the European Space Agency (ESA), under the visiting researcher programme. XD acknowledges funding within the framework of the Investissements d'Avenir programme (ANR-15-IDEX-02), through the funding of the 'Origin of Life' project of the Univ. Grenoble-Alpes. EM acknowledges funding from FAPEMIG under project number APQ-02493-22 and research productivity grant number 309829/2022-4 awarded by the CNPq, Brazil. BK acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme
(Grant agreement No. 865624). Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. The observations at the CFHT were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. We gratefully acknowledge the CFHT QSO observers who made this project possible. This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna; Astropy, 12 a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018); NumPy (van der Walt et al. 2011); Matplotlib: Visualization with Python (Hunter 2007); SciPy (Virtanen et al. 2020).

## References

Artigau, É., Cadieux, C., Cook, N. J., et al. 2022, AJ, 164, 84

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Bagnulo, S., Landolfi, M., Landstreet, J. D., et al. 2009, PASP, 121, 993

Bellotti, S. 2023, PhD Thesis, IRAP, France

Bellotti, S., Petit, P., Morin, J., et al. 2022, A&A, 657, A107

Bellotti, S., Fares, R., Vidotto, A. A., et al. 2023a, A&A, 676, A139

Bellotti, S., Morin, J., Lehmann, L. T., et al. 2023b, A&A, 676, A56

Bertaux, J. L., Lallement, R., Ferron, S., Boonne, C., & Bodichon, R. 2014, A&A, 564, A46

Boro Saikia, S., Jeffers, S. V., Morin, J., et al. 2016, A&A, 594, A29

Boro Saikia, S., Marvin, C. J., Jeffers, S. V., et al. 2018, A&A, 616, A108

Boro Saikia, S., Lüftinger, T., Folsom, C. P., et al. 2022, A&A, 658, A16

Browning, M. K. 2008, ApJ, 676, 1262

Carmona, A., Delfosse, X., Bellotti, S., et al. 2023, A&A, 674, A110

Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039

Chabrier, G., & Küker, M. 2006, A&A, 446, 1027

Charbonneau, P. 2010, Liv. Rev. Sol. Phys., 7, 3

Charbonneau, P. 2020, Liv. Rev. Sol. Phys., 17, 4

Claret, A., & Bloemen, S. 2011, A&A, 529, A75

Cook, N. J., Artigau, É., Doyon, R., et al. 2022, PASP, 134, 114509

Cortés-Zuleta, P., Boisse, I., Klein, B., et al. 2023, A&A, 673, A14

Cristofari, P. I., Donati, J. F., Folsom, C. P., et al. 2023, MNRAS, 522, 1342

da Gomes Silva, J., Santos, N. C., Bonfils, X., et al. 2012, A&A, 541, A9

Delchambre, L. 2015, MNRAS, 446, 3545

del Toro Iniesta, J. C. 2003, Introduction to Spectropolarimetry (Cambridge: Cambridge University Press)

Donati, J. F. 2003, in Solar Polarization, eds. J. Trujillo-Bueno, & J. Sanchez Almeida, ASP Conf. Ser., 307, 41

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

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

Donati, J. F., Mengel, M., Carter, B. D., et al. 2000, MNRAS, 316, 699

Donati, J. F., Catala, C., Landstreet, J. D., & Petit, P. 2006a, in Solar Polarization 4, eds. R. Casini, & B. W. Lites, ASP Conf. Ser., 358, 362

Donati, J. F., Howarth, I. D., Jardine, M. M., et al. 2006b, MNRAS, 370, 629

Donati, J. F., Morin, J., Petit, P., et al. 2008, MNRAS, 390, 545

Donati, J. F., Kouach, D., Moutou, C., et al. 2020, MNRAS, 498, 5684

Donati, J. F., Cristofari, P. I., Finociety, B., et al. 2023a, MNRAS, 525, 455

Donati, J. F., Lehmann, L. T., Cristofari, P. I., et al. 2023b, MNRAS, 525, 2015

Dorman, B., Nelson, L. A., & Chau, W. Y. 1989, ApJ, 342, 1003

Dungee, R., van Saders, J., Gaidos, E., et al. 2022, ApJ, 938, 118

Durney, B. R., De Young, D. S., & Roxburgh, I. W. 1993, Sol. Phys., 145, 207

Erba, C., Folsom, C. P., David-Uraz, A., et al. 2024, ApJ, submitted [arXiv:2401.09722]

Fares, R., Moutou, C., Donati, J. F., et al. 2013, MNRAS, 435, 1451

Finley, A. J., & Matt, S. P. 2017, ApJ, 845, 46

Folsom, C. P., Petit, P., Bouvier, J., et al. 2016, MNRAS, 457, 580

Folsom, C. P., Bouvier, J., Petit, P., et al. 2018, MNRAS, 474, 4956

Folsom, C. P., Ó Fionnagáin, D., Fossati, L., et al. 2020, A&A, 633, A48

Fouqué, P., Martioli, E., Donati, J. F., et al. 2023, A&A, 672, A52

Fuhrmeister, B., Czesla, S., Perdelwitz, V., et al. 2023, A&A, 670, A71

Gastine, T., Morin, J., Duarte, L., et al. 2013, A&A, 549, L5

Gray, D. F. 2005, The Observation and Analysis of Stellar Photospheres (Cambridge: Cambridge University Press)

Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951

Haywood, R. D., Milbourne, T. W., Saar, S. H., et al. 2022, ApJ, 935, 6

Hébrard, É. M., Donati, J. F., Delfosse, X., et al. 2016, MNRAS, 461, 1465

Huerta, M., Johns-Krull, C. M., Prato, L., Hartigan, P., & Jaffe, D. T. 2008, ApJ, 678, 472

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

Irving, Z. A., Saar, S. H., Wargelin, B. J., & do Nascimento, J.-D. 2023, ApJ, 949, 51

Jeffers, S. V., Boro Saikia, S., Barnes, J. R., et al. 2017, MNRAS, 471, L96

Kitchatinov, L. L., Moss, D., & Sokoloff, D. 2014, MNRAS, 442, L1

Klein, B., Donati, J.-F., Moutou, C., et al. 2021, MNRAS, 502, 188

Klein, B., Zicher, N., Kavanagh, R. D., et al. 2022, MNRAS, 512, 5067

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

Kochukhov, O., & Lavail, A. 2017, ApJ, 835, L4

Kochukhov, O., Makaganiuk, V., & Piskunov, N. 2010, A&A, 524, A5

Lafarga, M., Ribas, I., Reiners, A., et al. 2021, A&A, 652, A28

Landi Degl'Innocenti, E. 1992, in Magnetic Field Measurements, eds. F. Sanchez, M. Collados, & M. Vazquez, 71

Landi Degl'Innocenti, E., & Landolfi, M. 2004, Polarization in Spectral Lines (Dordrecht: Kluwer Academic Publishers), 307

Lavail, A., Kochukhov, O., & Wade, G. A. 2018, MNRAS, 479, 4836

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

Lehmann, L. T., Donati, J. F., Fouqué, P., et al. 2024, MNRAS, 527, 4330

Lienhard, F., Mortier, A., Cegla, H. M., et al. 2023, MNRAS, 522, 5862

Lopez-Santiago, J., Martino, L., Míguez, J., & Vázquez, M. A. 2020, AJ, 160, 273

Lorenzo-Oliveira, D., Freitas, F. C., Meléndez, J., et al. 2018, A&A, 619, A73

Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119

Maeder, A., & Meynet, G. 2000, ARA&A, 38, 143

Mavridis, L. N., & Avgoloupis, S. 1986, A&A, 154, 171

Meunier, N. 2021, arXiv e-prints [arXiv:2104.06072]

Mignon, L., Meunier, N., Delfosse, X., et al. 2023, A&A, 675, A168

Morin, J. 2012, in EAS Pub. Ser., eds. C. Reylé, C. Charbonnel, & M. Schultheis, 57, 165

Morin, J., Donati, J. F., Petit, P., et al. 2008, MNRAS, 390, 567

Morin, J., Donati, J. F., Petit, P., et al. 2010, MNRAS, 407, 2269

Morin, J., Delfosse, X., Donati, J. F., et al. 2011, in SF2A-2011: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. G. Alecian, K. Belkacem, R. Samadi, & D. Valls-Gabaud, 503

Mullan, D. J., & MacDonald, J. 2001, ApJ, 559, 353

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

Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9

Petit, P., Donati, J. F., & Collier Cameron, A. 2002, MNRAS, 334, 374

Petit, P., Louge, T., Théado, S., et al. 2014, PASP, 126, 469

Petit, P., Folsom, C. P., Donati, J. F., et al. 2021, A&A, 648, A55

Phan-Bao, N., Lim, J., Donati, J.-F., Johns-Krull, C. M., & Martín, E. L. 2009, ApJ, 704, 1721

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in FORTRAN. The Art of Scientific Computing (Cambridge: Cambridge University Press)

Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279

Rabus, M., Lachaume, R., Jordán, A., et al. 2019, MNRAS, 484, 2674

Rachkovsky, D. N. 1967, Izvestiya Ordena Trudovogo Krasnogo Znamenii Krymskoj Astrofizicheskoj Observatorii, 37, 56

Rees, D. E., & Semel, M. D. 1979, A&A, 74, 1

Reid, I. N., Cruz, K. L., Allen, P., et al. 2004, AJ, 128, 463

Reiners, A. 2012, Liv. Rev. Sol. Phys., 9, 1

Reiners, A., & Basri, G. 2009, A&A, 496, 787

Reiners, A., Shulyak, D., Käpylä, P. J., et al. 2022, A&A, 662, A41

Reinhold, T., & Hekker, S. 2020, A&A, 635, A43

Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr., 90, 054005

Saar, S. H. 1994, in Infrared Solar Physics, eds. D. M. Rabin, J. T. Jefferies, & C. Lindsey, 154, 493

Sanderson, T. R., Appourchaux, T., Hoeksema, J. T., & Harvey, K. L. 2003, J. Geophys. Res. (Space Physics), 108, 1035

Schrijver, C. J., & Zwaan, C. 2000, Solar and Stellar Magnetic Activity (Cambridge: Cambridge University Press)

See, V., Matt, S. P., Finley, A. J., et al. 2019, ApJ, 886, 120

Segura, A., Walkowicz, L. M., Meadows, V., Kasting, J., & Hawley, S. 2010, Astrobiology, 10, 751

Semel, M. 1989, A&A, 225, 456

Shulyak, D., Reiners, A., Wende, S., et al. 2011, in 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, eds. C. Johns-Krull, M. K. Browning, & A. A. West, ASP Conf. Ser., 448, 1263

Shulyak, D., Reiners, A., Engeln, A., et al. 2017, Nature, 5, 0184

Shulyak, D., Reiners, A., Nagel, E., et al. 2019, A&A, 626, A86

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

Skumanich, A. 1972, ApJ, 171, 565

Stelzer, B., Marino, A., Micela, G., López-Santiago, J., & Liefke, C. 2013, MNRAS, 431, 2063

Suárez Mascareño, A., Rebolo, R., & González Hernández, J. I. 2016, A&A, 595, A12