© The Authors 2023

1 Introduction

Stellar activity is driven by an underlying magnetic dynamo that is predicted to change periodically according to the theory of an α-Ω dynamo (Roberts & Stix 1972). Dynamo simulations for the fully convective M dwarf Proxima Centauri also directly led to magnetic cycles (Yadav et al. 2016). This periodic behaviour manifests itself as stellar activity cycles. These cycles have also been observed on stars other than the Sun by time series of X-ray and chromospheric activity indicators, or photometry.

While the Sun varies in X-rays as a coronal indicator by more than a factor of ten during its cycle (Peres et al. 2000), only a few stars have been found to exhibit long-term cyclic behaviour in X-rays. These cycles all had a much lower amplitude (Robrade et al. 2012; Wargelin et al. 2017; Coffaro et al. 2020).

Chromospheric indicators, on the other hand, have been used widely for a cycle search in solar-type stars. The monitoring program conducted at the Mt. Wilson Observatory (Wilson 1978) using the line index of the Ca II H & Klines revealed that stellar activity cycles are a ubiquitous phenomenon (Baliunas et al. 1995). Short and long cycles have been found for several stars (Brandenburg et al. 2017), starting with the Sun, which exhibits the approximately 80-yr-long Gleissberg cycle (Gleissberg 1945) and the 11-yr-long Schwabe cycle (Richards et al. 2009), although it also shows a transient periodic behaviour of about one-half to 2 yr (Ballester et al. 2002; Mursula et al. 2003). Multiple cycles have also been found in other stars, such as in ε Eridani, which displays two distinct cycles at about 3 and 12 yr (Brandenburg et al. 2017; Jeffers et al. 2022b). In addition to these decade-long cycles, F-type stars with activity cycles shorter than 1 yr have been revealed (Jeffers et al. 2018; Mittag et al. 2019). Although some authors sort different cycles lengths into an inactive and active branch, the existence of these two branches has been questioned by Boro Saikia et al. (2018) and Brown et al. (2022), for example.

While extensive work about activity cycles has been conducted on FGK stars (see e.g. Lovis et al. 2011), only a few M dwarfs have been under investigation, mostly using photometric time series. For example, Suárez Mascareño et al. (2016) derived magnetic cycles and rotation periods in 125 late-type stars, 70 of which were M dwarfs. About half of them exhibited long-term photometric variability. Another study that used photometric data to derive rotation periods and cycles was conducted by Díez Alonso et al. (2019). They found 12 stars with long-term cyclic variations with an observation time baseline of more than 9 yr. Moreover, Küker et al. (2019) found that in a sample of 31 M dwarfs, 19 exhibited a cycle. Using chromospheric indicators, fewer cycles were found for M dwarfs by Robertson et al. (2013) by employing time series of Hα indices. They also identified linear and quadratic trends as possible parts of activity cycles, which are too long to be covered by the available data. While they searched 93 K and M dwarfs, only six periods and seven long-term trends were identified by Robertson et al. (2013). Using Hα and Ca II H&K data, Suárez Mascareño et al. (2018) found that in a sample of 71 M dwarfs, 13 exhibited a cycle. Moreover, Gomes da Silva et al. (2012) found a correlation between long-term activity variations measured by a Na I D index and radial velocity (RV) in about one-third of the 27 stars they analysed.

This long-term correlation between an activity indicator and RV may impose problems for exoplanet detections with longer orbital periods, as is already known from short time-activity variations caused by the rotation period of the star, for instance, which may hide or mimic a planetary RV signal (Queloz et al. 2001; Dumusque et al. 2011; Jeffers et al. 2022a). This is caused by the asymmetry that is imprinted on the spectrum by surface heterogeneities. The asymmetry in turn leads to a shift of the line centre depending on the size and location of the active region on the stellar surface. Accordingly, Stock et al. (2020) found a long-term periodic behaviour in their RV measurements, which were mainly conducted with the CARMENES spectrograph. In the framework of an exoplanet search, they found tentative activity cycles of 600 and 2900 days for the two M dwarfs GJ 251 and Lalande 21185, respectively.

This study searches for more activity cycles in M dwarfs using data from the CARMENES spectrograph for different chromospheric line measurements. The paper is structured as follows: in Sect. 2 we present the data, their reduction, and the method for the measurement of our chromospheric indicators. In Sect. 3, we detail our search methods for activity cycles and trends. In Sect. 4, we present and discuss our results. We give a summary and concluding remarks in Sect. 5.

Fig. 1 Distribution of the observation time baseline for our sample of 211 stars with more than 19 observations and a minimum observation time baseline of 200 days.

2 Observations and measurements of chromospheric indicators

All spectra used for the present analysis were taken with the CARMENES spectrograph, installed at the 3.5 m Calar Alto telescope (Quirrenbach et al. 2020). CARMENES is a spectrograph covering the wavelength range from 5200 to 9600 Å in the visual channel (VIS) and from 9600 to 17 100 Å in the near-infrared channel (NIR). The instrument provides a spectral resolution of ~ 94600 in VIS and ~80400 in NIR. The CARMENES consortium is conducting a survey of ~350 M dwarfs to find low-mass exoplanets (Alonso-Floriano et al. 2015; Reiners et al. 2018). Since the cadence of the spectra is optimised for the planet search, usually no continuous time series are obtained.

2.1 Data reduction and sample construction

In our analysis, we considered a sample of 362 M dwarfs observed by CARMENES, resulting in more than 19000 spectra taken before February 2022. From this sample, we excluded known binaries (Baroch et al. 2018; Schweitzer et al. 2019; Baroch et al. 2021). Binaries may hamper our analysis because the line shifts induced by binarity may be large enough to shift parts of the line out of the integration range for line indices or pseudo-equivalent width (pEW). Additionally, we had to constrain the stellar sample further because the CARMENES data are optimised for exoplanetary search, and therefore not all stars are observed with the same frequency or time baseline (i.e. the elapsed time between the first and last observation). We wish to search for long-term systematic variations, therefore we restricted the sample to all stars with a time coverage of more than 200 days with at least 20 observations to allow for a statistically meaningful analysis. Even though it is hard to detect periods with so few observations, linear and quadratic trends may still be detected. In this fashion, we obtained a sub-sample of 211 stars for our analysis. We show the distribution of the observation time baseline for these 211 stars in Fig. 1.

The stellar spectra were reduced using the CARMENES reduction pipeline (Zechmeister et al. 2014; Caballero et al. 2016). Subsequently, we corrected them for barycentric and RV motions and carried out a correction for telluric absorption lines (Nagel et al. 2019) using the molecfit package¹. No correction for airglow emission lines was attempted because they do not play a role for the activity indicators we used.

2.2 Measurement of activity indicators

To assess the activity state of the stars in each spectrum, we employed pEW measurements because the spectra of M dwarfs do not show an identifiable continuum (because molecular absorption lines are present nearly everywhere). We measured the pEW of the blue Ca II infrared triplet (IRT) line at 8500.35 Å (hereafter CaIIIRT₁) and of the middle Ca II IRT line at 8544.44 Å(hereafter Ca II IRT₂), and the Hα line (all wave-lengths are given for vacuum conditions). The three Ca II IRT lines act very similarly regarding activity (Lafarga et al. 2021), therefore we used only the two bluer lines because they are located on the same order of the spectrograph.

The pEW was calculated using the expression

$$\text{pEW}={\displaystyle \int \left(1-F_{\text{core}}/F_0\right)\text{d}\lambda ,}$$(1)

where F_core is the flux density in the line band, F₀ is the mean flux density in the two reference intervals (representing the pseudo-continuum), and λ denotes the wavelength. Thus, time series of pEW measurements were obtained for each star. The integration ranges for the individual lines are given in Table 1.

The uncertainty of individual pEW measurements depends on the S/N of the spectra, which can be particularly problematic for the Hα line, where the S/N is typically better for bright, early M dwarfs than for fainter mid-type M dwarf. We did not compute pEWs when the S/N in the pseudo-continuum was lower than 15. Therefore, no pEW(Hα) measurements could be obtained for some of the latest and faintest M dwarfs, and we omitted these stars from further analysis when fewer than 20 measurements were available. When the S/N was good enough to compute pEWs, we calculated their error by a bootstrap analysis using the flux density errors from the pipeline to re-compute the pEWs 1000 times. This yielded quite low statistical errors; in the case of pEW(Hα) and pEW(Ca II IRT), the relative statistical error is about 10⁻⁷−10⁻⁵ for most of the spectra, while the time series of the chromospheric indicators of inactive stars exhibits a larger scatter of about 10% of the median value of the pEW in Hα and of 3% in both pEW(CaIIIRT). This is apparently caused by intrinsic variations of the activity level of the stars, and we therefore applied this as error to all chromospheric indicator measurements. We also display this scatter as error bars in all the figures showing time series.

In addition to these classical chromospheric indicators, we also considered more recently introduced activity indicators, specifically, a TiO index defined as the ratio of the average flux density in a band red- and blue-wards of the band head, as introduced by Schöfer et al. (2019). The wavelength bands, which slightly differ from Schöfer et al. (2019), are given in Table 1. An error was again estimated by the variation in inactive stars to be 0.3% (which roughly agrees with a bootstrap analysis in this case). Moreover, we used RV, chromatic index (CRX; as a measure of change in RV shift as a function of wavelength), and the differential line width (dLw; as a measure of changes in line width compared to a constructed average). All three indicators and their errors were obtained for the VIS arm of CARMENES from the reduction pipeline serval (Zechmeister et al. 2018). They have been used for example in a search for rotation periods by Lafarga et al. (2021).

For a sub-sample of 186 stars with Admeasurements compiled by Perdelwitz et al. (2021) from seven different instruments (most spectra were obtained with the HARPS², HIRES³, and Narval spectrographs), we included these data in our cycle search. Perdelwitz et al. (2021) used the pipeline-reduced spectra from the respective archives and computed the chromospheric line fluxes in the Ca IIH & K lines after rectification and subtraction of the photospheric flux using PHOENIX (Hauschildt et al. 1999) photospheric model spectra from the Husser et al. (2013) database. While some of the obtained ${R^{\prime }}_{\text{HK}}$ time series temporally overlap with the CARMENES data, most of them do not. More-over, the observational cadence is usually coarser, while the time baseline is much longer than for the CARMENES data. This longer baseline makes the data more suitable for a cycle search.

2.3 Compilation of cycles from the literature

From the literature, we collected a list of 57 proposed cycles in 41 stars as linear or quadratic trends for our sample stars. Most of these cycles originate from an analysis of photometric data. For example, the studies by Suárez Mascareño et al. (2016), which were dedicated to a cycle search, and of Díez Alonso et al. (2019), where some cycles emerged as by-products of a rotation period search, used only photometric data. We also included cycles found in photometric data by Küker et al. (2019). In contrast, Robertson et al. (2013) used an Hα index for an activity cycle search, and Perdelwitz et al. (2021) also performed a period analysis on their ${R^{\prime }}_{\text{HK}}$ data with three published cycles. We also have some stars in common with the sample of Suárez Mascareño et al. (2018), who used the Ca II H&K and Hα index for their cycle search.

Other activity cycles have been published as part of RV exo-planet searches that also used the CARMENES data, for example by Stock et al. (2020) for Lalande 21185 (including RV data from the SOPHIE instrument) and GJ 251 (including data from the HARPS instrument), for which cycle lengths of 2900 days and 600 days were found. Moreover, Lopez-Santiago et al. (2020) found a probable cycle length of 14 yr for the star GJ 3512 by performing a combined fit to the rotation period and activity cycle.

3 Searching for systematic long-term variations

We used several methods to search for systematic long-term variations in the available data to account for the different time scales involved. While we searched for cycles that were potentially longer than a decade, our longest time baseline is only 6.1 a (= years). This implies that we may cover only parts of activity cycles. Therefore, we also searched for linear and quadratic trends, which may arise when the data cover only the gradual decay or increase in activity or include an activity maximum or minimum. We caution, nevertheless, that these systematic long-term variations need not to lead to periodic behaviour. Therefore, we distinguish between periodic behaviour, where we do cover at least one full period, and linear and quadratic trends, which are indicative of systematic, but not necessarily cyclic behaviour.

Moreover, in all our search algorithms, we decided to tie a detection to at least three indicators and used in turn not as strict criteria for the individual indicators, as we would have applied if we had accepted detections for single indicators. We opted for this alternative since correlated (or anti-correlated) behaviour between the indicators is expected in the case of systematic long-term variability since they (optimally) all should react to the changing activity level. We decided in this context to use at least three indicators, because for only two indicators, the two used CallIRT lines may trigger a detection alone. With more than three indicators, this criterion is very strict because it may lead to non-detections when the different indicators have different sensitivities.

Finally, before searching for cycles or parts of it, we applied a 3σ clipping to each of the indicator time series using the python routine scipy.sigmaclip. Thus we avoided outliers caused by flaring, bad weather, or instrumental effects, which may confound the search. Moreover, after conducting the automatic variability search, we excluded all time series that had data gaps spanning more than 40% of the temporal baseline, were longer than 550 days, or both, because we noted visually that long data gaps like this were often problematic.

3.1 Searching for linear trends

The linear trend search was conducted as follows: we used each of the pEW of Hα, Ca II IRT₁, Ca II IRT₂, and the I_TiO indices as well as CRX and dLw to compute a linear regression with respect to time, employing scipy. stats. linregress in python, which outputs also Pearson's correlation coefficient r and the corresponding p-value. This additional output allowed us to identify trends in time. These trends were searched for using the whole data set for each star because we are interested in finding stars for which the whole time series covers the activity increase or decrease phase of an activity cycle. The identified trends are linear by the definition of the Pearson's correlation coefficient, but a subsequent visual inspection showed that in a few cases, the linear description only holds as a first-order approximation, but a parabolic description appears to be more appropriate. These time series were also found in our quadratic trend search and were marked accordingly. Therefore, we claim a linear trend when the Pearson correlation coefficients are abs(r) > 0.55 and p < 0.002 in three or more of the indicators. We determined these thresholds empirically using a test sample of 20 stars that included 3 stars whose linear trends were found by eye. The relatively low threshold of abs(r) > 0.55 was justified a posteriori because we did not reject any of the automatically found trends by eye (but only found some to be better described by quadratic trends).

An instructive example of the trends we found is shown in Fig. 2 for the M1.0V star J23245+578/BD+57 2735. The pEW(Hα) and pEW(Ca II IRT₁) both show about the same slope. This is the case for most of the trends we found. In one of the time series, we nevertheless found one convincing example in which the trends of Hα and the pEW(Ca II RT₁) were of opposite sign. In the M0.0 V star J13450+176/BD+18 2776, we found a highly significant Pearson correlation coefficient of 0.95 (p-value 10⁻¹⁵) between values of pEW(Ca II IRT₁) and pEW(Ca II IRT₂), while the correlation of the values of pEW(Hα) and the two pEW(Ca II RT) are −0.73 and −0.69, respectively (p-values < 3.1 × 10⁻⁵). We show the corresponding time series in Fig. 3 with the pEW(Ca II IRT₂) instead of I_Ti0 as an exception. This behaviour may be explained by the relative inactivity of the star. With an increase in activity level in M dwarfs, we first expect a deepening of the Hα line, which is followed by a fill-in, until the line finally enters emission at even higher activity levels (Cram & Mullan 1979). For the CallIRT lines, only a fill-in is expected instead. This behaviour can lead to an anti-correlation of the two lines for the lowest activity stars, where Hα still deepens with increasing activity.

Fig. 2 Time series of J23245+578/BD+57 2735. We show the pEW(Hα) (top, black dots), pEW(Ca II IRT1 ) (middle, blue dots), and ITiO (bottom, magenta dots). Additionally, the best linear fit is indicated as the solid line. The error bars are not statistical, but determined by the activity jitter of inactive stars.

3.2 Searching for quadratic trends

We used a parabolic fit to identify another type of systematic long-term variability. We applied a collection of ten criteria, nine of which had to be fulfilled to accept a quadratic trend. We used a reduced χ² like quantity to access the quality of the fit. However, the correspondence is not exact owing to our use of a generic jitter estimate and the absence of the degrees of freedom in the definition: (i)${\chi }^2=\frac{1}{n_{\text{spec}}}\frac{{\displaystyle \sum {\left({\text{pEW}}_i-\text{poly}\right)}^2}}{{\displaystyle \sum {\text{err}}_i^2}}<1.05$ for pEW(Hα) (ii-iii) The same holds for Ca II IRT₁ and Ca II IRT₂ with a threshold of 1.0 (iv) The same holds for the TiO index with a threshold of 1.03. To fine-tune these first four criteria, we visually inspected a sub-sample of 20 time series, which included five time series for which a quadratic trend was found by eye. (v)–(viii) All maxima or minima of the parabolic fit had to lie within the observed time baseline, (ix) The maximum difference between the corresponding dates of the maxima or minima of the quadratic fits of the different activity indicators must be lower than 15% of the observation time baseline and shorter than 150 days, with the exception of one indicator, (x) Same as (ix), but without any exception. These last two criteria make use of the expected correlated behaviour of the different indicators and at least (ix) must be triggered for an automatic detection.

As an instructive example of our quadratic trend search, we show the M0.0V star J14257+236W/BD+24 2733A in Fig. 4. This clearly demonstrates the need of longer time baselines to clarify whether these quadratic trends lead to periodic behaviour or are just systematic episodes embedded in chaotic variations.

Fig. 3 Time series of J13450+176/BD+18 2776. We show the pEW(Hα) (top, black dots), pEW(Ca II IRT1) (middle, blue dots), and pEW(Ca II IRT2) (bottom, green dots). Additionally, the best-fit linear trend is indicated as the solid line.

Fig. 4 Time series of J14257+236W/BD+24 2733A. We show pEW(Hα) (top) and pEW(Ca II IRT1) (bottom). Additionally, the best-fit quadratic trend is indicated as the solid line.

3.3 Searching for cycles

To search for cyclic behaviour, we applied the generalised Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) as implemented in PyAstronomy⁴ (Czesla et al. 2019) to each of our chromospheric indicator time series, namely pEW(Hα), pEW(Ca II IRT₁), pEW(Ca II IRT₂), I_Ti0, ${R^{\prime }}_{\text{HK}}$ (whenever available), RV, CRX, and dLw. In our further analysis, we searched for periods longer than 1 yr and shorter than 90% of the time baseline of the individual data set to cover whole cycles. We accepted periods with an FAP lower than 0.1% as significant. When we found significant periods in three or more indicators, we verified whether the difference between the maximum peak periods of at least three indicators (not necessarily the significant ones) was less than 20%, to ensure that the indicators showed (about) the same period. In this case, we accepted the period as an automatic detection. We also included ${R^{\prime }}_{\text{HK}}$ data in this procedure, but we caution that due to the usually much longer time baseline, we often found longer periods in these data as maximum peak. The search was therefore basically conducted on the CARMENES data alone. Since we allow for a 20% difference of the periods in the individual chromospheric indicators we assign a 20% error to all our measured cycles.

In an extended search, we also accepted periods shorter than lyr but longer than 150 days (to avoid rotation periods). This led to an automatic detection of nine additional stars, which were all rejected by visual inspection. They are usually caused by the observing pattern or instead are weak quadratic trends. The only example of a star with a 0.7 a period, where the short period may be present in the TiO index alone, is the M 5.0 star J20260+585/Wolf 1069, which we show in Fig. A.1. Nevertheless, this GLS also shows a small peak in the window function that is caused by the observing scheme at about the period we found. We therefore also rejected this star by visual inspection.

All automatically found cycle candidates were inspected visually. This sometimes revealed that the period we found can also be explained by a quadratic or linear trend, especially when the period we found is near the length of the time baseline of the observations.

An example of an formerly unknown cycle candidate that fulfilled our automatic detection criteria and was not excluded by visual inspection is the M0.1 V star J22330+093/BD+08 4887, for which we show the GLS and time series in Fig. 5. While pEW(Hα) and pEW(Call IRT) do show about the same period, the TiO index does not show the cycle period, but a shorter period is preferred. Since the broad peak in the GLS reveals that the period is not well constrained and no ${R^{\prime }}_{\text{HK}}$ data are available, the activity cycle period needs to be confirmed by photometry or other spectroscopic data.

Fig. 5 GLS and time series of J22330+093/BD+08 4887. We show the GLS (left, top) with pEW(Hα) in black, pEW(Ca II IRT1) in blue, pEW(Ca II IRT2) in green, the TiO index in magenta, RV in yellow, dLw in grey, CRX in orange, and the window function in red. The dashed blue line indicates FAP=0.001 for the pEW(Ca II IRT1) data. Additionally, we show the time series for pEW(Hα) (right, top), for pEW(Ca II IRT1) (left, bottom), and the TiO index (right, bottom). The dashed black lines indicate the best sine fit for each indicator.

4 Results and discussion

With the methods described in the previous section, we found 26 automatic detections of long-term variability in our sample of 211 stars. Seven of these are linear trends, 14 are quadratic trends, and 12 are cycles (some stars exhibit more than one type of variability). From these numbers, we already excluded one linear trend, one quadratic trend, and one cycle, whose time series showed data gaps longer than 550 days and longer than 40% of the observation time baseline. Our results are summarised in Table 2, where we list the Karmn number of the star, a common name, the spectral type, the number of spectra we used, and the time baseline of the CARMENES observations. Furthermore, in the case of automatically found linear trends, we list the highest Pearson correlation coefficient (which does not necessarily correspond to the lowest or even a significant p-value), and in case of a quadratic trend, the number of flags that were triggered. In the case of significant GLS periods, we list the period length, and when available, the period length found in the ${R^{\prime }}_{\text{HK}}$ data. Additionally, we cite known literature values.

After the automatic detection, we inspected all variable time series by eye to determine the most appropriate detection in cases of multiple detections. For example, in all three cases when we rejected the linear trend, a quadratic trend was also found (twice by automatic detection, once by eye), and the quadratic trend was always preferred by visual evaluation. Six of the seven rejected cycles also showed a quadratic trend (four by automatic detection, two by visual inspection). Two of the stars with automatically detected quadratic trends (J04290+219/BD+21 652, J05314-036/HD 36395) also have a longer cycle detected in the ${R^{\prime }}_{\text{HK}}$ data, and the cycle length we found is about the length of the observation time baseline.

Finally, four linear trends, 11 quadratic trends, and five cycle periods passed the visual inspection. Therefore, 20 of our 26 automatically found long-term variable stars pass the visual inspection, which additionally leads to an appropriate categorisation. Since our adopted errors are relatively large, we refrained from applying an F-test and performed the sorting by eye. The stars that passed visual inspection are marked with an asterisk in Table 2, where we also list some remarks from the visual inspection for many stars.

4.1 Comparison to literature

For seven of the stars showing systematic variability, we could find published cycle lengths. Unfortunately, the situation is complicated and we find only few agreements. For J02222+478/BD+47 612, we find a quadratic trend that agrees with the linear trend found by Robertson et al. (2013). The same is true for J16254+543/GJ 625, where we find a quadratic trend that agrees with the cycle period found by Suárez Mascareño et al. (2018). For J05314-036/HD 36395, a long and a short cycle is known. Since the time baseline of the CARMENES data is about 1.5 a shorter than the shorter cycle, the quadratic trend agrees with this short cycle, while the period found in the ${R^{\prime }}_{\text{HK}}$ data agrees with the long cycle.

For J19169+051N/V 1428 Aql, a short and a long cycle are known from photometry. As the observation time baseline of the CARMENES data is rather short for this star, we find a quadratic trend, while in the ${R^{\prime }}_{\text{HK}}$ data, we find a period that roughly agrees with the longer literature period. For J22565+165/HD 216899, we find a period of 5.4 a without visual confirmation that is much shorter than the period found in the ${R^{\prime }}_{\text{HK}}$ data. J06105-218/HD 42581 A and J07274+052/Luyten's star are discussed in Sect. 4.2 because we found a visually convincing cycle period.

4.2 Promising activity cycle candidates

In the following, we discuss the five stars that are the most promising candidates for activity cycles for which more than a full cycle is covered by the CARMENES observations and that passed the visual inspection. The first of these stars is J22330+093/BD+08 4887. This star was discussed in Sect. 3.3, see Fig. 5.

In the M1.0 dwarf J05415+534/HD 233153 a rather short period of only 1.4 a was found that we show in the left panel of Fig. 6 in Hα and Ca II IRT₁. The trend that was also detected automatically is overlaid. It may indicate two different cycle periods. However, the TiO index does not show any variation, not even the trend.

For the M3.5 V star J18346+401/LP 229-017, we show the GLS and time series in the right panel of Fig. 6. Although the period has about the same length as the observation time baseline, we opted for the cycle and not for a quadratic trend by visial inspection in this special case. The period can be seen quite clearly in Hα and Ca II IRT, but not in TiO, where a shorter period is found.

The case for the M0.5 dwarf J06105-218/HD 42581 A is slightly more complicated. A cycle of about 8.3 a is known from the literature (Suárez Mascareño et al. 2016; Díez Alonso et al. 2019). In the ${R^{\prime }}_{\text{HK}}$ data, we find a peak at about twice this period, which may be a sub-harmonic because there is also a (non-significant) peak at 7.9 a, which would agree with the literature value. In the CARMENES data, we also find a cycle of 2.7 a length, which is visually confirmed. Either this is a quasiperiodic episode, or this is a second, shorter cycle. We show the data of this star in Fig. A.2.

For the M3.5 dwarf J07274+052/Luyten's star, we find a period of 2.8 a, in agreement with the findings from the ${R^{\prime }}_{\text{HK}}$ data, which has been published previously (Perdelwitz et al. 2021). There is a second published period of 6.6±1.3a (Suárez Mascareño et al. 2016) that may be about twice the shorter cycle. A second significant peak lies at 6.1 a in the GLS of the ${R^{\prime }}_{\text{HK}}$ data. The correct period cannot be decided from the available data, which are shown in Fig. A.3.

4.3 Known cycles without detected long-term variation in the chromospheric indicators

We fail to find systematic long-term activity variations with our automatic search algorithms for 34 stars with proposed cycles from the literature. These non-detections fall into three categories: First, the star was excluded from our data analysis due to binarity (1 star) or fewer than 20 observations (4 stars); or it showed possibly problematic long data gaps (1 stars). Second, visual inspections revealed indications of linear or quadratic trends or even cycles (18 stars) in individual chromospheric indicators. Third, we do not find any systematic long-term variations (10 stars). We list all these stars in Table 3 and also remark on the visual inspection of the CARMENES data.

An example that demonstrates the difficulties well is the star J06371+175/HD260655. We were able to visually identify a linear trend in the pEW(Ca II IRT) time series, while we find a quadratic trend in pEW(Hα) and I_TiO. Since the baseline of the CARMENES observations is 6.1 a, this neither agrees with the 7.3 a cycle from the ${R^{\prime }}_{\text{HK}}$ data nor with literature values of 2–3 a from Díez Alonso et al. (2019) or 4.9 a from Suárez Mascareño et al. (2016). Since the trend we found in the data indicates a cycle that is longer than 10 yr, the star may exhibit multiple cycles. We caution, however, that the two cycle length values from the literature (both photometry) also disagree, and claims of more than two cycles per star would remain highly speculative given the available data. Alternatively, this may indicate quasi-periodic episodes or changes in the cycle length, such as those known from the Sun.

Another instructive example is the star J06548+332/Wolf 294/GJ 251. Stock et al. (2020) found a significant 600 d cycle in RV data using HIRES data alone, while the signal is not present in their CARMENES RV data. Their inspection of chromospheric indicators reveals a 660 days cycle in the Hα indicator and a 300 d cycle length in CRX and NalD. While we find a significant peak at about 630 days in the TiO index and pEW(Ca ll IRT₂), the GLS of the pEW(Hα) shows an even higher peak at about 1680 days (including newer data), while the ${R^{\prime }}_{\text{HK}}$ data favour about twice that period with a 3300 days period, which may indicate that the period in pEW(Hα) is a harmonic, the longer period of which cannot be found because the observation time baseline is too short. Further investigation is required to determine whether there are indeed a short (600 days = 1.6 a) and a long cycle (9 a).

The star J11033+359/Lalande 21185 is a peculiar case, for which we find a visually convincing period in Hα, but not in any other chromospheric indicator. This period was also reported with a similar length by Stock et al. (2020) for their (shorter time baseline) CARMENES data alone. However, they replaced that period for a period that was more than twice longer based on their combined CARMENES and SOPHIE data analysis. This value in turn is about half the value we found from the ${R^{\prime }}_{\text{HK}}$ data.

On the other hand, there are examples such as J07446+035/YZ CMi, J10196+198/AD Leo, J17303+055/BD+05 3409, J20525–169/LP816-060, or J22096–046/BD–05 5715, where we find a linear or quadratic trend that agrees well with the cycles proposed in literature by visual inspection in single activity indicators. Moreover, the star J18580+059 /BD+05 3993 was discussed in detail by Suárez Mascareño et al. (2018). Our observation started just after the time baseline discussed by them, but a downward trend was expected to occur afterwards, which is what we observe, but only by visual inspection.

Although the existence of multiple cycles in a single star may explain some of the discrepancies we found, the star may also undergo phases of quasi-systematic episodes and therefore show a different behaviour at different times. Another problem is that different activity indicators reveal different behaviour. Trends or periods are often only found in one activity indicator, while the others show a constant, chaotic, or even contradicting systematic behaviour. While the first two possibilities can be explained by the different sensitivity of the indicators, the latter is not expected. It may indicate time lags between the indicators in the (quasi-)systematic development. Possible examples for different sensitivity are the stars J08161+013/GJ 2066 and J11033+359/Lalande 21185, for which we find a quadratic trend or even a cycle in a single activity indicator, while the others are rather constant. An example in which activity indicators show a different behaviour is J06548+332/Wolf 294 (discussed above, linear trend in pEW(Ca II IRT), quadratic trend in pEW(Hα)). Since the discrepancy is seen only in the last observing block, this indicates a time lag between the indicators.

Another problem may be imposed by the observation scheme. An example is J08298+297/DX Cnc, for which the observation pattern shows some tight blocks with larger gaps in between. Together with the rather high activity level of J08298+297/DX Cnc, this may be the cause for our finding only insignificant peaks in the GLS at about the known cycle period in pEW(Hα) and pEW(Ca II IRT).

Fig. 6 Same as in Fig. 5. GLS and time series of J05415+534/HD 233153 (left) and J18346+401/LP 229-017 (right). We note that for J05415+534/HD 233153, a linear trend was found as well, which we include in the sine fit for pEW(Hα) and pEW(Ca II IRT1). We further note that for J18346+401/LP 229-017, the highest peak of pEW(Hα) in the GLS is at about twice the frequency as for the other chromospheric indicators that triggered the cycle detection. We therefore used half the frequency for the sine fit of pEW(Hα) as well.

4.4 Promising cycle candidates from ${R^{\prime }}_{HK}$ alone

Since the ${R^{\prime }}_{\text{HK}}$ data usually show a much longer observation time baseline and ${R^{\prime }}_{\text{HK}}$ is known to be suitable for activity cycle search in more solar-type stars (Baliunas et al. 1995; Brandenburg et al. 2017), we decided to perform an additional search on the ${R^{\prime }}_{\text{HK}}$ data alone. This can be achieved even though the data originate from seven different telescopes, because instrument offsets were minimised by extracting the ${R^{\prime }}_{\text{HK}}$ values using a direct comparison to model spectra (for further discussion, see Perdelwitz et al. 2021). Since only a single chromospheric indicator was used this time and a cross-check with other indicators is therefore missing, we opted in this search for a much better FAP < 0.0005 and excluded all stars with fewer than 20 spectra or automatically detected periods shorter than 1 yr. Since the ${R^{\prime }}_{\text{HK}}$ data often have a much coarser sampling, we also excluded stars with data gaps longer than a quarter of the observation time baseline and longer than 1000 days.

We repeated the GLS analysis and found that 22 cycle candidates fulfilled these criteria. Out of these, 7 correspond to stars that also show a variability detection in the CARMENES data and are therefore listed in Table 2. Another 6 correspond to stars with a proposed cycle from the literature and are therefore listed in Table 3. A cycle period detection based on the ${R^{\prime }}_{\text{HK}}$ data alone is indicated in both of the tables. The remaining 9 cycle candidates are detections based on ${R^{\prime }}_{\text{HK}}$ alone, and none of these was rejected by visual inspection. We list them in Table 4 with their spectral type, number of observations, the ${R^{\prime }}_{\text{HK}}$ observation time baseline, and remarks on the visual inspection of the CARMENES data. We show the GLS and time series of ${R^{\prime }}_{\text{HK}}$ of the most promising 6 candidates in Figs. A.4, A.5, and A.6.

This better performance of the ${R^{\prime }}_{\text{HK}}$ data in the cycle search is partly caused by the longer time baseline, which allows us to find cycles that are not yet covered by the CARMENES data. We therefore state the CARMENES baselines in Table 4, where appropriate. Nevertheless, in some of these cases, no trends are found, either. While the reason may be that the CARMENES observations incidentally cover the maximum or minimum of the cycle, the different sensitivity of the lines may also play a role. We discuss this in Sect. 4.6 and reveal that the relative amplitude in ${R^{\prime }}_{\text{HK}}$ data is higher than of the other line pEWs.

Fig. 7 Number of stars per spectral sub-type (blue bars, left y-axis), and number of stars with an automatic variability detection (red bars, left y-axis). We show the fraction of stars with variability per spectral sub-type as black dots corresponding to the right y-axis.

4.5 Dependence on spectral type

To identify a possible dependence on spectral sub-type, we show in Fig. 7 the number of stars and the number of automatic long-term variability detections (including the nine detections based on ${R^{\prime }}_{\text{HK}}$ alone). The percentage of systematically variable stars drops toward later spectral types from about 15–30% in early-M dwarf stars to lower than 5% for mid-M dwarfs. We do not find any sign of a systematic long-term variability for any star with spectral type later than M4.0 V. This cannot be caused by observational biases because the mean observation time baseline of the included M4 stars is even longer than for M0 stars. Although the mean number of observations for M4 stars is some-what lower than for M0 stars, the number of stars with more than 50 usable observations is higher for M4 stars than for M0 stars.

Noise should not play a role here either because we excluded spectra with a too low level of the signal-to-noise ratio, and this is usually only an issue for stars with spectral type M5 and later because the exposure time never exceeds 1800 s. Either these latest type, fully convective stars do not have activity cycles, or they may not found due to frequent flaring.

Flaring was also identified as the main reason for variability in fully convective M dwarfs in a study by Medina et al. (2022), who found only for one very young star a correlation of the EW(Hα) to the rotational phase in a sample of 13 stars. The short timescales of steallar variability in their study, which was observed with high cadence, suggest an origin of the variability in flares. Although we excluded large flares, frequent low-level flaring would lead to a broadening of the distribution of the individual indicators and would thus not be handled by σ clipping. This can lead to a veiling of possible activity cycles.

The stars with cycle values proposed in the literature also include stars with spectral type M4 and later. We therefore note that all these stellar cycles have been detected by photometry and not by chromospheric indicators. At least for J07446+035/YZ CMi and J20525-169/LP 816-060, visual inspection also revealed weak linear or parabolic trends in these stars, while for J08298+297/DX Cnc, we find an insignificant period of the cycle length from the literature. This may indicate that in these relatively late-M dwarfs, the sensitivity of the chro-mospheric indicators to cyclic activity is lower than in early-type M dwarfs. Next to flaring, this may be caused by generally high chromospheric filling factors, for example, which would lead to low amplitudes in the cycle modulation and might therefore hamper the detection of a cycle.

4.6 Sensitivity of the different indicators

We analysed the chromospheric indicator that was most often involved in the trend or cycle detection for the different search methods. For the linear trend search, the CaIIIRT lines, Hα, and dLw most often led to a trend with between 5 and 7 detections each, while the TiO index and CRX led to one detection, and RV to zero detection. For the parabolic search, we find that pEW(Hα), pEW(Ca II IRT), and I_TiO each were triggered for all detected stars. Finally, for the cycle search, again the Ca II IRT lines were involved most often in a cycle detection (12 detections each), followed by Hα with 8 detections, dLw with 6 detections, TiO and CRX with 4 detections, and RV with 3 detections. We conclude that long-term variability is most easily found in the Ca II IRT lines, followed by Hα and dLw, while TiO, RV, and CRX are less well suited.

For a better comparison, we also considered the amplitudes of the different indicators. We computed the amplitude as the relative difference of the median of the ten highest and ten lowest data points in the time series with a detected (linearorquadratic) trend or period (regardless of whether the indicator triggered the detection). For the time series showing a trend, the amplitude is only a lower threshold. This leads to median amplitudes of 14% for pEW(Hα), of 4% for pEW(Ca II IRT), and of 0.08% for I_TiO. This underlines the higher sensitivity of Hα compared to Ca II IRT and especially I_TiO. The highest sensitivity is found for ${R^{\prime }}_{\text{HK}}$, with a relative median amplitude of 46%. The relatively low amplitudes of pEW(Ca II IRT) come as a surprise because Ca II IRT triggered the most automatic detections. This better performance compared to pEW(Hα) may be caused by a larger influence of flares on the pEW(Hα), because the line is formed at larger heights in the chromosphere than Ca II IRT.

4.7 Comparison of cycles to rotation periods

While Böhm-Vitense (2007) found two distinct relations between rotation period and cycle length for active and inactive G and K stars, Küker et al. (2019) found virtually no dependence of the cycle length on rotation period for early-M dwarfs. They found two groups with characteristic cycle times, however, first the fast rotators with rotation periods shorter than one day and cycle lengths of about 1 yr, and second, stars with longer rotation periods and a cycle length between 3 and 5 yr. Moreover, Suárez Mascareño et al. (2016) did not find a correlation between rotation period and cycle length for the M dwarfs included in their study. Nevertheless, they find a weak correlation for earlier-type stars. Suárez Mascareño et al. (2018) found a weak correlation of rotation period with cycle length for M dwarfs using literature values in addition to their measurements, but also stated that a correlation test showed that there is a 37% probability of the correlation being spurious.

Although we have 25 stars with known rotation periods shorter than one day in our sample, these stars are all of spectral type M3.5 or later, and no long-term variation is found for any of these stars. If the finding of Küker et al. (2019) is correct that such short cycles are found only for fast rotators, our visual evaluation that none of the found periods shorter than 1 yr is also consistent because none of these stars is a fast rotators.

We searched for a correlation of rotation period P_rot and cycle length P_cyc using our cycle detections from CARMENES and ${R^{\prime }}_{\text{HK}}$ data and all literature values from Table 3. This resulted in 34 stars (13 from our own measurements and 21 from the literature) for which rotation periods are also available. We performed an analysis similar to that of Suárez Mascareño et al. (2018) and computed the slope of log(P_cyc/P_rot) versus log(1/P_rot), which gives 1.1 ± 0.1 and therefore does not deviate significantly from 1.0. This implies that there is no correlation between rotation period and cycle length.

5 Conclusions

We searched in a sample of 211 M dwarfs for long-term systematic variability in chromospheric indicators. The observations were taken from the more than 19 000 CARMENES observations, considering only stars with more than 19 observations and an observation time baseline of more than 200 days. We not only considered the pEW of Hα and the Ca II IRT lines as chromospheric indicators, but also more recently defined activity indicators: the TiO-index, RV, CRX, and dLw. In a sub-sample of 186 stars, we also used already published ${R^{\prime }}_{\text{HK}}$data with much longer observation time baselines, originating from different instruments. In these we performed an automatic search for periods and other long-term systematics, namely linear trends and quadratic trends. The automatic search led to 26 detections of this systematic long-term variability. By visual inspection, we rejected 6 of these detections and sorted the remaining stars uniquely into our different variability categories, leaving us with four linear trends, 11 quadratic trends, and five periods. We show all the five stars with proposed periods in the Figs. 5, 6, A.2, and A.3. Analysing the ${R^{\prime }}_{\text{HK}}$ data alone and applying a stricter FAP threshold than for the combined data set, we found an additional nine promising cycle candidates, which are shown in Figs. A.4, A.5, and A.6. Furthermore, we found that from the analysed indicators that are available from the CARMENES observations, the pEWs of Ca II IRT and Hα lines are best suited for a study of long-term variability. Comparison of the relative median amplitudes of the systematic variation shows that ${R^{\prime }}_{\text{HK}}$ has the highest amplitude, followed by pEW(Hα), pEW(Ca II IRT), and ITiO. Therefore, we propose that future activity cycle searches performed in chromospheric indicators should concentrate on ${R^{\prime }}_{\text{HK}}$, pEW(Hα), or pEW(CaII IRT). On the other hand, we were unable to confirm many cycles proposed in the literature on the basis of photometric data, especially for later-type stars, which may suggest that photometry is even better suited for an activity cycle search than chromospheric indicators in these stars. Most important for RV exoplanet searches, we found that RV variations are not a sensitive indicator of systematic long-term chromospheric variability in our study.

The detection rate of systematic long-term variability in all our sample stars is about 12% from the CARMENES data alone. This is low in comparison to Suárez Mascareño et al. (2018), for instance, who found cycles in about 18% of their M0–M3 stars using chromospheric indicators as well. When we only consider our M0–M3 stars, we have an even slightly higher detection rate of about 25%, which is to be expected because we additionally searched for linear and quadratic trends. The detection rate of the long-term variability dramatically decreases along the M spectral sequence to later-type stars. For stars M4.0 and later, we do not find any activity cycle. These fully convective stars may either exhibit no cycles caused by the different dynamo underlying their activity phenomena, or the existence of cycles is veiled in the chromospheric indicators by frequent flaring. This latter possibility is underlined by the study of Suárez Mascareño et al. (2016), who found cycles in about 50% of their M dwarfs for fully convective stars as well, using photometry. Again, this suggests that photometry is better suited for cycle searches in M dwarfs.

Although the CARMENES observation span little more than 6 yr at most, we found some short-cycle candidates in agreement with literature values, or promising parts of cycles where sometimes longer (mostly photometrically determined) cycles were already known. We also identified a few cases for which modulation periods reported in the literature could not be confirmed. The most promising linear and quadratic trends stress that the observations need to be continued in the future to accumulate data that cover at least a whole cycle or reveal the current data to be a quasi-systematic episode.

Acknowledgements

B.F. acknowledges funding by the DFG under Schm 1032/69-1. We thank our referee for the careful reading and the suggestions for improvement. CARMENES is an instrument for the Centro Astronómico Hispanoen Andalucía de Calar Alto (CAHA, Almería, Spain). We acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG) through the priority programme SPP 1992 “Exploring the Diversity of Extra-solar Planets” (JE 701/5-1), and from the Agencia Estatal de Investigación 10.13039/501100011033 of the Ministerio de Ciencia e Innovación and the ERDF “A way of making Europe” through projects PID2019-109522GB- C5[1:4], PGC2018-098153-B-C33, and the Centre of Excellence “Severo Ochoa” and “María de Maeztu” awards to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), Instituto de Astrofísica de Andalucía (SEV-2017-0709), and Centro de Astrobiología (MDM-2017-0737), and the Generalitat de Catalunya/CERCA programme. CARMENES is funded by the German Max-Planck-Gesellschaft (MPG), the Spanish Consejo Superior de Investigaciones Científicas (CSIC), the European Union through FEDER/ERF FICTS-2011-02 funds, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Ciències de l'Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the Spanish Ministry of Economy, the German Science Foundation through the Major Research Instrumentation Programme and DFG Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía.

Appendix A Further examples of long-term variability

Fig. A.1 Same as in Fig. 5. GLS and time series for J20260+585/Wolf 1069. In pEW(Hα), frequent flaring and the spacing of the observation leads to the period we found, while the same period may be present in the ITiO data.

Fig. A.2 Same as in Fig. 5, but for I06105–218/HD 42581 A.

Fig. A.3 Same as in Fig. 5, but for I07274+052/Luyten's star. We omit the error bars for clarity.

Fig. A.4 GLS and time series of J10251–102/BD–09 3070 (M1.0 V, left) and J10289+008/BD+01 2447 (M2.0 V, right). In the GLS (top), we show the power of the ${R^{\prime }}_{\text{HK}}$ data in cyan, the dashed cyan line marks the 0.0005 FAP level, and the red line is the power of the window function. In the time series (bottom), we show the ${R^{\prime }}_{\text{HK}}$ values as cyan dots and the best-fit sine function as the dashed black line.

Fig. A.5 Same as Fig. A.4, but for J11421+267/Ross 905 (M2.5 V, left) and J14294+155/Rossl30 (M2.0 V, right).

Fig. A.6 Same as Fig. A.4, but for J20450+444/BD+44 3567 (M1.5 V, left) and J20533+621/HD 199305 (Ml.O V, right).

References

All Tables

All Figures

	Fig. 1 Distribution of the observation time baseline for our sample of 211 stars with more than 19 observations and a minimum observation time baseline of 200 days.
In the text

	Fig. 2 Time series of J23245+578/BD+57 2735. We show the pEW(Hα) (top, black dots), pEW(Ca II IRT1 ) (middle, blue dots), and ITiO (bottom, magenta dots). Additionally, the best linear fit is indicated as the solid line. The error bars are not statistical, but determined by the activity jitter of inactive stars.
In the text

	Fig. 3 Time series of J13450+176/BD+18 2776. We show the pEW(Hα) (top, black dots), pEW(Ca II IRT1) (middle, blue dots), and pEW(Ca II IRT2) (bottom, green dots). Additionally, the best-fit linear trend is indicated as the solid line.
In the text

	Fig. 4 Time series of J14257+236W/BD+24 2733A. We show pEW(Hα) (top) and pEW(Ca II IRT1) (bottom). Additionally, the best-fit quadratic trend is indicated as the solid line.
In the text

Fig. 5 GLS and time series of J22330+093/BD+08 4887. We show the GLS (left, top) with pEW(Hα) in black, pEW(Ca II IRT1) in blue, pEW(Ca II IRT2) in green, the TiO index in magenta, RV in yellow, dLw in grey, CRX in orange, and the window function in red. The dashed blue line indicates FAP=0.001 for the pEW(Ca II IRT1) data. Additionally, we show the time series for pEW(Hα) (right, top), for pEW(Ca II IRT1) (left, bottom), and the TiO index (right, bottom). The dashed black lines indicate the best sine fit for each indicator.

In the text

Fig. 6 Same as in Fig. 5. GLS and time series of J05415+534/HD 233153 (left) and J18346+401/LP 229-017 (right). We note that for J05415+534/HD 233153, a linear trend was found as well, which we include in the sine fit for pEW(Hα) and pEW(Ca II IRT1). We further note that for J18346+401/LP 229-017, the highest peak of pEW(Hα) in the GLS is at about twice the frequency as for the other chromospheric indicators that triggered the cycle detection. We therefore used half the frequency for the sine fit of pEW(Hα) as well.

In the text

	Fig. 7 Number of stars per spectral sub-type (blue bars, left y-axis), and number of stars with an automatic variability detection (red bars, left y-axis). We show the fraction of stars with variability per spectral sub-type as black dots corresponding to the right y-axis.
In the text

	Fig. A.1 Same as in Fig. 5. GLS and time series for J20260+585/Wolf 1069. In pEW(Hα), frequent flaring and the spacing of the observation leads to the period we found, while the same period may be present in the ITiO data.
In the text

	Fig. A.2 Same as in Fig. 5, but for I06105–218/HD 42581 A.
In the text

	Fig. A.3 Same as in Fig. 5, but for I07274+052/Luyten's star. We omit the error bars for clarity.
In the text

Fig. A.4 GLS and time series of J10251–102/BD–09 3070 (M1.0 V, left) and J10289+008/BD+01 2447 (M2.0 V, right). In the GLS (top), we show the power of the ${R^{\prime }}_{\text{HK}}$ data in cyan, the dashed cyan line marks the 0.0005 FAP level, and the red line is the power of the window function. In the time series (bottom), we show the ${R^{\prime }}_{\text{HK}}$ values as cyan dots and the best-fit sine function as the dashed black line.

In the text

	Fig. A.5 Same as Fig. A.4, but for J11421+267/Ross 905 (M2.5 V, left) and J14294+155/Rossl30 (M2.0 V, right).
In the text

	Fig. A.6 Same as Fig. A.4, but for J20450+444/BD+44 3567 (M1.5 V, left) and J20533+621/HD 199305 (Ml.O V, right).
In the text