Stellar Cycles in Fully Convective Stars and a New Interpretation of Dynamo Evolution

ASAS-3 optical photometry data were downloaded from the ASAS All-Star Catalogue. ⁴ All ASAS-3 data use V-band filters, and individual measurements are ascribed one of four grades: A—best data; B—mean data; C—A and B without a measured indication; D—worst data. We opted to use only measurements with an A or B grade. For most stars in our sample, the MAG_2 (four pixel) aperture data had the least scatter and therefore the lowest empirical uncertainty; to be consistent, we therefore used this aperture for all our stars. We note, however, that the choice of aperture can sometimes have a significant effect on the inferred cycle period, and could explain (at least partly) why our results sometimes differ from those of Suárez Mascareño et al. (2016). ⁵ ASAS-3 has a FWHM point-spread function (PSF) of typically 23'', corresponding to approximately 1.5 pixels.

ASAS-SN optical photometry data were collected using the light-curve server. ⁶ Early ASAS-SN data were gathered using V-band filters, while the most recent data use g band. The ASAS-SN light-curve server does not perform proper-motion corrections, unlike the ASAS All-Star Catalogue, so proper-motion-corrected coordinates had to be specified for each query. ASAS-SN has a 16'' FWHM PSF, so we allowed for a proper motion of up to 1'' either side of the input coordinates; extractions use a two pixel (16'') aperture. We then compiled these proper-motion-corrected light curves into single light curves for each star once we had downloaded all available data.

TESS optical photometry data were downloaded using the Lightkurve Python package (Astropy Collaboration et al. 2018; Lightkurve Collaboration et al. 2018; Brasseur et al. 2019; Ginsburg et al. 2019), which accesses publicly available data directly from the Mikulski Archive for Space Telescopes (MAST). Only data products produced via their Science Processing Operations Center Science Analysis Pipeline (Jenkins et al. 2016) were used. For each star, we also used the pipeline aperture.

2.2.1. ASAS-SN

Initially, the early V-band and recent g-band data were disjointed in the compiled ASAS-SN light curves. We therefore needed to color-correct these data. Using the empirical UBVR-uvgr relation found by Kent (1985) and improved by Windhorst et al. (1991), we can convert from V- into g-band magnitudes using

$$\begin{eqnarray}&&{\rm{V}}={\rm{g}}-0.03-0.42({\rm{g}}-{\rm{r}}),\end{eqnarray} \tag{ 1 }$$

where V is the V-band magnitude, g is the g-band magnitude, and g − r is the g–r color index, given by

$$\begin{eqnarray}&&g-r=1.02(B-V)-0.22,\end{eqnarray} \tag{ 2 }$$

where B − V is the B − V color index. However, M stars are relatively cool objects and emit mostly toward the red end of the optical spectrum; as such, the B − V color index is not the best diagnostic for these stars. In addition, the relation given in Equation (1) was derived using predominantly O–K stars. After color-correcting V-band data using these formulae, residual offsets were therefore required. These offsets were computed by comparing the means of overlapping g- and V-band data.

Further internal cross-calibration is then required because ASAS-SN data are collected using multiple telescopes positioned around the world, with small but sometimes significant differences between nominally identical instruments. The ASAS-SN pipeline automatically adjusts for these differences, but color-dependent residuals may remain, as we found for our reddish M stars.

To correct for differences between different ASAS-SN telescopes, we first selected the telescope with the largest number of measurements, thus defining our reference data set. We then calibrated the remaining telescopes sequentially in descending order of relative overlap. To calibrate a telescope, we applied an offset such that the averages of overlapping data points between the telescope and the reference data set were equal. Once a telescope had been calibrated, it then became part of the reference data set for subsequent calibrations. Fortunately, we did not find any instances of telescopes with no overlapping data. However, if a particular telescope was found to show unusual behavior within a given light curve (relative to the other telescopes) then this telescope's data were removed. For example, telescope "bA" showed an unusually large scatter for GJ 447, obscuring its cyclic modulations, but appeared perfectly normal for GJ 581, GJ 729, and GJ 849; this telescope's measurements were therefore removed in the former case, but kept in the latter three.

2.2.2. TESS

2.2.3. Outlier Removal

When assessing our results, it is important to keep in mind that some target-star intensity extractions may suffer from contamination by binary companions or other stars that appear nearby in the sky. Winters et al. (2019) list data from an all-sky volume-limited (25pc) survey for stellar companions to 1120 M-dwarf primaries, and serves as our principal reference regarding binarity because it includes all our stars of interest except for Proxima (GJ 551). Two systems are listed as binaries: GJ 234 and GJ 896.

The GJ 234 pair has a semimajor axis of 11 (e.g., Gatewood et al. 2003; Kervella et al. 2019), but contamination by the secondary is small, with estimated intensity ratios of 100:5.3 in the B band (Gatewood et al. 2003), 100:5.9 in V (Henry et al. 1999), and 100:25 in K (Coppenbarger et al. 1994). GJ 896 A and B had a separation of 5.35'' in 2004 (Winters et al. 2019) with an estimated V magnitude of 10.29 for component A and 2.12 mag dimmer for B (100:14 ratio). Davison et al. (2015) measured a separation of 7'' and listed the spectral types as M3.5+M4.0. These stars are mostly blended in ASAS and ASAS-SN data, and it is difficult to determine if the lack of an obvious cycle in GJ 896A is intrinsic, from contamination by the secondary star, or variable extraction efficiencies as the pair's center of emission moves around its center of mass and proper motion.

The multiplicity of GJ 273 is somewhat unclear, with Vrijmoet et al. (2020) identifying it as a suspected double based on Gaia data, while Winters et al. (2019) conclude that the putative companion is probably an unassociated background object. Astudillo-Defru et al. (2017) analyzed HARPS RV data and report that GJ 273 has four planets, without any mention of spectral evidence for a companion; in our analysis, we treat this star as a single.

Although GJ 273, along with GJ 317 and GJ 406, were not included in their study of binarity in nearby stars using data from HIPPARCOS and Gaia DR2, Kervella et al. (2019) calculated minimum (inclination-dependent) masses for companions to the rest of the stars in our sample. Except for the two known binary systems, minimum masses were no higher than 0.25M_J assuming zero inclination and orbital radii of r = 1 au, with mass scaling as $\sqrt{r}$. All our stars have clear rotational modulation, indicating that inclination cannot be too small (assuming rotational and orbital coplanarity), so it seems unlikely that any of the stars under study here (except for GJ 234 and GJ 896) have stellar mass companions. Many of these stars are also the subjects of radial velocity (RV) studies, which would readily reveal binarity unless the orbital inclination were very small. We particularly mention GJ 358 because of its remarkable changing cycle periods: Bonfils et al. (2013) explicitly treat this star as a single after excluding known and discovered spectroscopic binaries and visual pairs from their analysis of HARPS data.

Intensity extractions can also be affected by stars that appear nearby on the sky even if they are not physically associated. GJ 406 has a very high proper motion (472 yr⁻¹) and has recently passed quite close to two objects that may provide contaminating flux. Likewise, Proxima lies in a very crowded part of the sky and is currently moving at 386 yr⁻¹ within a group of stars that, although dimmer, may be bright enough to cause significant contamination within the ASAS-3 and ASAS-SN extraction regions. These effects should be slowly varying, but they complicate interpretation of observations and may be at least partly responsible for the general brightening seen in Proxima's ASAS and ASAS-SN data.

A further consideration to make when assessing our results is which stars are fully convective. As mentioned in Section 1, spectral type is often used to distinguish partially from fully convective M stars. However, spectral type determinations are not always a robust metric, and the boundary between subtypes is often blurred. Indeed, some of these difficulties were mentioned in Section 1 and are discussed in detail in Section 2.3. Fortunately, recent studies have alleviated some of the ambiguity in identifying fully convective M stars.

Using data from Gaia DR2, Jao et al. (2018) identified a gap in the color–magnitude diagram for the lower main sequence, which they proposed could be used to distinguish partially from fully convective M stars. van Saders & Pinsonneault (2012) first predicted this gap as the result of nonequilibrium ³He fusion prior to stars becoming fully convective; using numerical simulations, Feiden et al. (2021) were able to reproduce this gap. Using this distinction, specifically in the V–K_S versus M_V plane, we identify 11 of our 15 M stars as being fully convective (see Table 1). We note, however, that for a few high proper-motion stars not included in Gaia DR2, we had to use DR3 parallaxes.

Herein, we predominately use traditional time-series analysis techniques, specifically, Lomb–Scargle (L–S) periodograms (Lomb 1976; Scargle 1982), for two key reasons. First, more advanced techniques, such as Gaussian processes (GPs), come with a significant computational expense. Depending on the number of covariance function hyperparameters, the number of mean function parameters, the priors, and the size of the data set itself, fitting a GP model can easily take several orders of magnitude times the time it takes to obtain an L–S periodogram. Second, as we discuss below, the literature shows that stellar cycles often appear to be well approximated by sines. We note two exceptions, however: GJ 358 and GJ 551 (Proxima), whose cyclic modulations cannot adequately be described by one or two single-period sinusoids; these stars are therefore better suited to a GP analysis.

To compute the false-alarm probabilities (FAPs) on our detections, we use the approximation proposed by Baluev (2008), often referred to as the Baluev estimate. We compared this method with a more robust bootstrap simulation method for six ASAS-3 and ASAS-SN light curves, and found that results were usually within 10%, and many times within 5%; the largest difference we found was FAP_Baluev: FAP_bootstrap ≈ 0.79. ⁸ Due to the computational expense of bootstrap simulations, without much improvement to performance, we prefer the Baluev estimate. In our work, we consider an FAP of ≤0.1% to constitute a significant detection; the same threshold was used in Suárez Mascareño et al. (2016).

In theory, typical time-series analysis techniques, such as fast Fourier transforms (FFTs) and L–S periodograms, may not necessarily be the best tools for the analysis of stellar cycles. Both are limited in that they model signals with sinusoids; while this may be a good model for idealized systems, it often fails in reality. Indeed, stellar cycles need not be sinusoidal in nature, nor strictly periodic.

To identify stellar cycles in unevenly sampled light curves, Olspert et al. (2018b) introduced a Bayesian generalized Lomb–Scargle periodogram with trend (BGLST); using synthetic quasiperiodic data, they showed that their BGLST performed better than the standard L–S periodogram. Moreover, in their follow-up paper, Olspert et al. (2018a) further showed that GPs can be used effectively in the analysis of stellar cycles - in their case, using Ca ii H and K data.

The greatest benefit of GPs over traditional time-series analysis techniques, such as the L–S periodogram, is their flexibility. Much like L–S periodograms, GPs can be used with unevenly sampled data. However, GPs need not assume the data contain a periodic, sinusoidal signal. The assumptions of a GP are encoded in its mean and covariance functions, which, in theory, allows for any number of arbitrary assumptions to be made.

Olspert et al. (2018a), in their study of stellar cycles, compared the performance of a Bayesian harmonic (i.e., sinusoidal) regression model, and GP models with periodic and quasiperiodic covariance functions. They found that where traditional methods suggest double cycles, their quasiperiodic GP model often found only a single cycle; as the assumption of a quasiperiodic signal has more physical justifications, Olspert et al. (2018a) concluded that double cycles are rarer than initial results (e.g., Wilson 1978; Baliunas et al. 1995; Brandenburg et al. 1998; Saar & Brandenburg 1999) suggest.

However, while a quasiperiodic model is more physically justified, its results are often not much of an improvement over a simple harmonic model. Olspert et al. (2018a) identified harmonic signals in the light curves in 36 stars; quasiperiodic signals were also identified in 26 of these 36 stars. Comparing these 26 quasiperiodic periods with their harmonic equivalents, 22 were comfortably within 3σ. This suggests that while stellar cycles need not be strictly periodic nor sinusoidal, they can often be well approximated by sines.

2.6.1. Rotation Periods from Photometric Time Series

To determine rotation periods, we used L–S periodograms, as well as data from both ASAS-SN and TESS. ASAS-3 data are noisier and lower cadence than ASAS-SN data, and therefore provide no additional benefit for rotational studies; for these reasons, ASAS-3 data are not included. The first step in our analysis was to separate the ASAS-SN data by observing season. We then computed L–S periodograms for each of these seasons, and, provided the period of the peak was greater than 15 days, used the tallest peak in each season's L–S periodogram to fit a sine model using the Levenberg–Marquardt (L–M) algorithm (a sophisticated nonlinear least-squares fitting algorithm; see Moré 1978)—as shown in Figure 1 for GJ 54.1. Note that there were often no significant peaks at P_rot/2.

Zoom In Zoom Out Reset image size

We discarded periodograms computed from seasons with low data density or incomplete seasons, ⁹ and identified significant peaks in each periodogram. After disregarding anomalous peaks, such as the broad peak at approximately 180 days in the third row of Figure 1, we inferred the rotation period by computing the mean period of the remaining peaks, and estimated the uncertainty as the standard deviation on this mean.

However, if the rotation period of a star was sufficiently short (shorter than 15 days), we used the higher cadence but shorter duration TESS data (TESS observation windows typically span ∼27 days) to better constrain this period. If a star had multiple TESS observations, we found that in every case, the inferred periods were identical to within the (tightly constrained) uncertainties.

Our analysis using TESS data proceeded in much the same way as our analysis using ASAS-SN data; an example is shown in Figure 2 for GJ 234. Due to our currently limited understanding of TESS systematics, however, the identification of rotation periods is restricted to below approximately half a TESS observing window (hence our 15 days cutoff)—see, for example, Anthony et al. (2022) and Claytor et al. (2022).

Zoom In Zoom Out Reset image size

To infer the rotation period of a star from its TESS light curve, we fit a Gaussian to the tallest peak in its periodogram. Our estimate for the rotation period was thus the mean of this Gaussian, with the uncertainty given by its standard deviation. Since we were often limited to one or two TESS observations, following the same procedure as described above for ASAS-SN resulted in unrealistically small errors in some cases and no errors in others; we therefore adopted this alternative method. We note that when we tested this approach on ASAS-SN data, it resulted in larger uncertainties due to the broad peaks in these periodograms. However, TESS observations only span ∼27 days, and so these broadening effects (e.g., resulting from differential rotation and multiple starspot generations) are minimal. TESS is also a space-based observatory, with a much higher signal-to-noise ratio than ASAS-SN. Given these differences, we deem this alternative method to be appropriate. Results are presented in Table 1 and are discussed in Section 3.

2.6.2. Cycle Periods from Photometric Time Series

To determine cycle periods, we typically use L–S periodograms incorporating data from both ASAS-3 and ASAS-SN; exceptions to this are Proxima (for which we also have ASAS-4 data and use GPs), and GJ 358 (for which we also use GPs)—see Section 2.5 for more details.

For each star, we computed a L–S periodogram from its light curve. We then identified significant peaks in periodograms with periods of more than 1.5 yr, and used the periods of these peaks to fit sine functions to the light curve; see Figures 3 and 4 for examples. All ASAS-3/-4 light curves and associated periodograms (16 images) are provided in Figure Set 3; likewise for ASAS-SN (15 images) in Figure Set 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (16 panels)

Tarball of figure set images

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (15 panels)

Tarball of figure set images

If a periodogram contained multiple significant peaks, we fit sine functions using each peak and then combined these functions to create a superposition. If the fit of the superposition was worse than the fit of any single component, then the worst fitting components were removed until the fit of the superposition was better than the fit of any single component, or only the best-fitting component was remaining. If the period(s) of the remaining component(s) was within the period of observation, then we deem this cycle to be "well defined" - and poorly constrained otherwise.

Using the analysis procedures described in Section 2.6.1, high-confidence rotation periods were determined for all but two of our stars. For one exception, GJ 729, P_rot is shorter than 15 days, but TESS data were not available for our work. We therefore used the same method as for stars with longer rotation periods, computing the mean and standard deviation of similar significant peaks in each ASAS-SN observing season's periodogram; these periodograms are shown in Figure 5. This resulted in a rotation period of 2.9 ± 0.1 days, in agreement with Suárez Mascareño et al. (2016) and Ibañez Bustos et al. (2020), who found a rotation periods of 2.9 ± 0.1 days (from photometric monitoring) and 2.848 ± 0.001 days (from chromospheric indicators), respectively.

Zoom In Zoom Out Reset image size

For the second exception, GJ 273, we could not determine a rotation period from our data sources. Breaking ASAS-SN data down by observing season, as described in Section 2.6.1, and shown in Figure 6, produced only one significant peak at ∼42 days. However, no other observing seasons showed any significant peaks at all. Computing periodograms from the entire ASAS-SN light curve, as shown in Figure 7, also produced no significant peaks with a period shorter than one year. Figures 6 and 7 show that this star has fewer ASAS-SN observations than many of the other stars in our sample, especially in the early seasons (which use V-band filters). The sparsity of these observations undoubtedly contributed to our inability to identify rotational modulations within this light curve. For these reasons, we adopted the P_rot of 115.9 ± 19.4 days measured by Suárez Mascareño et al. (2015) using HARPS Ca ii H and K and Hα data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given analysis uncertainties and the inherent variation in rotational measurements caused by starspot evolution, migration in latitude and longitude, and differential rotation, our results are generally in good agreement with previous measurements. Where there are significant disagreements (GJ 234, GJ 628, and LP 816-60), we believe our measurements are more trustworthy given the higher quality and cadence of TESS and ASAS-SN data compared to the ASAS-3 data used by Suárez Mascareño et al. (2016). Note that two of our stars, GJ 317 and GJ 406, have P_rot values reported for the first time.

L–S analysis yields FAPs shorter than 0.1% using ASAS-3/-4 and/or ASAS-SN data for 12 stars, not including a few cases where the significant peak corresponds to a cycle period exceeding the period of observation. Below, we comment on three stars of particular note. For stars not mentioned below, light curves and associated L–S periodograms are provided in Figure Sets 1 (ASAS-3/-4) and 2 (ASAS-SN).

3.2.1. LP 816-60

3.2.2. GJ 358

The GJ 358 ASAS-SN light curve (Figure 8) shows that the cyclic modulations appear to be reasonably well described by two sinusoidal functions, following the P_cyc and P_cyc/2 rule suggested by do J.-D. Nascimento et al. (2023, in preparation). Interestingly, both the ASAS-3 and ASAS-SN light curves (Figures 9 and 8) show long-term linear trends whose gradients are of opposite sign; this may suggest that GJ 358 exhibits a longer-term Gleissberg-esque cycle, although many more years of data will be needed to confirm this.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Further study of the light curves from ASAS-3 and ASAS-SN, both of which are shown in a single plot in Figure 10, reveals that the GJ 358 cycle period has been steadily decreasing for at least the last 20 yr. ASAS-3 data show cyclic modulations with periods greater than four years, while the most recent modulation seen in ASAS-SN data has a period of just 1.5 yr. This star therefore appears to be quite unusual. Unfortunately, there is a ∼6 yr gap between ASAS-3 and ASAS-SN data on this object. This gap could be filled with ASAS-4 data, but these data are not publicly available. Analysis of these data would prove very interesting and might reveal, for example, whether the cycle period has changed monotonically, and when and how the longer-term trends seen in ASAS-3 and ASAS-SN data intersected.

Zoom In Zoom Out Reset image size

At present, it is also unclear whether GJ 358's cycle is approaching some kind of (quasi-)steady state, or if this cycle period will reach a local minimum before increasing again. It is also unclear if this changing period is related to the long-term cycle suggested by the linear trends seen in the GJ 358 ASAS-3 and ASAS-SN data. Future observations of this object are therefore needed to better understand this seemingly exceptional cyclic activity.

3.2.3. GJ 551

GJ 551 (Proxima Centauri) is an interesting case for several reasons. It is the most well-studied example of a fully convective M star exhibiting cyclic activity, and it is the only fully convective star (in fact, the only M star) to have undergone long-term X-ray monitoring. Suárez Mascareño et al. (2016) found a 7 yr cycle in ASAS-3 data, and Wargelin et al. (2017) included an additional ∼5 yr of ASAS-4 data, plus UV and X-ray observations to further support this conclusion.

More recent optical monitoring data, however, are less clear, as reflected in the ASAS-SN results listed in Table 1 and shown in Figure 11, which includes ASAS-4 data extending into 2019 that were presented in Damasso et al. (2020). ASAS-3 and ASAS-4 data were cross calibrated (private communication, G. Pojmański), and then calibrated versus ASAS-SN data following the "overlapping data" procedures described in Section 2.2.1. As seen in Figure 11, periodic behavior is quite apparent at the beginning, but steadily weakens while average brightness increases. We suspect that the increasing stellar contamination noted in Section 2.3 is significantly affecting the data, but attempting to model and correct the contamination is beyond the scope of this paper. In contrast, X-ray data from Swift, now spanning more than 12 yr (with some gaps), indicate well-behaved cyclic behavior, currently with a period of ∼9 yr (B. J. Wargelin et al. 2023, in preparation).

Zoom In Zoom Out Reset image size

It should be noted that some of the cycles presented in Table 1 may not be related to a star's global magnetic field, but could be the result of other physical phenomena (or, if they are, may not be Schwabe cycles; do J.-D. Nascimento et al. 2023, in preparation). For example, Rossby waves (horizontally flowing eddies in the convective zone/photosphere), as seen in the Sun, can act as dynamos (e.g., Gilman 1969), leading to perceived cyclic activity that is not necessarily related to the star's global magnetic field.

In the case of our Sun, there is also the Rieger cycle (Rieger et al. 1984), where hard solar flares are observed to occur in groups with a mean spacing of ∼154 days. Similar cycles surely exist in other stars; while flares are a result of magnetic activity, they are local effects, and it is not well understood how they relate to a star's global magnetic field - even in the case of our Sun (e.g., Toriumi & Park 2022).

One way to assess the validity of a star's apparent cycles is to see if they are consistent through ASAS-3 and ASAS-SN data. If the best-fitting cyclic models from both data sets show a reasonable level of agreement (in terms of both phase and amplitude) where they overlap, then these cycles are more likely to be real. Note, however, that this method assumes a roughly constant cycle period, which, as discussed in Section 1, may not always be the case.

To investigate the amplitude and phase matches, we extend the cycles identified in Table 1 such that the ASAS-3 and ASAS-SN cyclic models for each star overlap by one year. In addition to this, we also offset the ASAS-3 data such that the mean magnitude of this light curve is the same as the mean magnitude of the ASAS-SN light curve, thus allowing easy comparison of the models' phases and amplitudes. Example plots are shown in Figure 12, showing the results for GJ 273, GJ 234, and GJ 447. All phase-match plots (10 images) are provided in Figure Set 12. Note, however, that we place little emphasis on these plots; indeed, they are primarily included so that the inferred cycles from both data sets can be more easily compared.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (10 panels)

Tarball of figure set images

As can be seen in the top panel of Figure 12, both the ASAS-3 and ASAS-SN cyclic models on GJ 273 agree very well, in terms of both phase and amplitude, where they overlap. We interpret this as enhanced evidence for these cycles being real. On the other hand, GJ 234 (middle panel of Figure 12) shows good amplitude agreement, but poorer phase agreement, while GJ 447 (bottom panel of Figure 12) shows poor agreement in both cases. In the final column of Table 1, we assign each star's ASAS-3 and ASAS-SN cycles a subjective phase-match grade of either good, fair, or poor. For reference, we judge that GJ 273, GJ 234, GJ 447 show good, fair, and poor phase agreement, respectively. The remaining plots are provided in Figure Set 12.

The convective turnover time is a measure of how long it takes for convection to move material from the bottom of the convective zone to the top, or vice versa, and it is related to the ratio of the thickness of the convective envelope (thin in F stars, fully convective by ∼M3.5) and the mean convective velocity (which increases with temperature). To compare stars with differently sized convective zones, temperatures, and rotation periods, a dimensionless quantity known as the Rossby number, Ro, is commonly used. Another dimensionless quantity is the linear αΩ dynamo number, D, a measure of the dynamo's strength, which is proportional to Ro⁻². The Rossby number of a star is defined as

$$\begin{eqnarray}&&{Ro}=\displaystyle \frac{{P}_{\mathrm{rot}}}{\tau },\end{eqnarray} \tag{ 3 }$$

where P_rot is the rotation period, and τ is the convective turnover time. Note, however, that there are two definitions for the convective turnover time: a local and a global definition. The local convective turnover time, τ_L , applies to the bottom of the convective zone (in the tachocline, where the main dynamo amplification takes place in some models), while the global convective turnover time, τ_G , is an average over the entire convective zone. The global convective turnover time is therefore a better quantity for M stars, where the tachocline is negligible or absent.

The choice of a convective turnover time is important. We need both a local (for tachocline-based dynamos) and a global (for full convection zone dynamos) τ, which we denote τ_L and τ_G, respectively. Since we are interested in the operation of dynamos, a τ derived from any indirect activity diagnostic (e.g., Ca ii HK or X-ray emission) is of dubious utility. This is because these empirical τs implicitly assume that the function $f(M,\mathrm{log}g,[\mathrm{Fe}/{\rm{H}}])$ that best connects the given emission to rotation can be defined to be τ, ignoring the complex physics connecting magnetic flux to the heating that produces the emission in question (e.g., Cuntz et al. 1999). The hidden inclusion of this additional heating physics makes these empirical τs inappropriate for understanding how dynamos themselves operate. Note that a τ derived by best-fitting rotation with unsigned magnetic flux measurements would be more appropriate. This has not yet been attempted to our knowledge, likely because data on unsigned fluxes are sparse and generally come with large errors (e.g., Saar et al. 1994; Reiners 2012). Calculating such a τ is beyond the scope of the current work; however, new flux measurement methods may make it feasible in the near future (e.g., Lehmann et al. 2015; Mortier 2016).

Unfortunately, theoretically based τs are also problematic. Due to uncertainties concerning the M dwarf internal structure near the fully convective limit (Baraffe & Chabrier 2018; Jao et al. 2022), any purely theoretical τ is at best approximate at low masses (Jao et al. 2022). There are also problems defining a local τ in the center of fully convective stars. Indeed, with no convection zone bottom, the tachocline concept driving the need for a local τ requires modification itself. Fortunately, a recent paper has worked to improve τ models in low-mass stars (Landin et al. 2023). We adopt their values here at age ≈1Gyr. This is post-zero-age main sequence (ZAMS) for all masses here except 0.1 M_⊙, but we note from their Figure 1 that τ_G(0.1M_⊙) evolves negligibly for ages >0.3 Gyr. Indeed, over the mass range considered here, τ varies little along the main sequence, except at the highest masses, where differences of ∼30% may accrue at age extremes. We use the V–K_S versus T_eff relation of Pecaut & Mamajek (2013; and later improvements; Mamajek, E.E. 2021, private communication) to match to our stars. As a cross check, we compare the calculations at low mass with a different kind of empirical τ calculation based purely on stellar properties - derived equating the star's output bolometric flux with convective flux (Corsaro et al. 2021). They find ${\tau }_{{\rm{G}}}\propto R{[M/({LR})]}^{1/3}\propto {M}^{1/3}{T}_{\mathrm{eff}}^{-4/3}$ for fully convective stars, where M, L, and R are the stellar mass, luminosity, and radius, respectively (taken from Pecaut & Mamajek 2013). We find that a scaled version of these physics-based empirical τ_G for the most part compares well with the new τ_G of Landin et al. (2023; Figure 13; where we also display their τ_L). The exception is in the narrow temperature range 3250 K < T_eff < 3400 K, very near the fully convective boundary, where τ_G (Landin) shows a spike that τ_G (Corsaro) does not. (We also note that T_eff(0.1M_⊙) is cooler in Landin et al. 2023). We suggest that τ_G calculations near the stellar center may still have problems (Jao et al. 2022), and in keeping with the idea that τ_L and τ_G should roughly scale with each other (e.g., Montesinos et al. 2001), we adopt the Corsaro values for T_eff ≤ 3400 K. We use τ_G for fully convective stars, and either τ_G or τ_L otherwise, depending on the circumstances.

Zoom In Zoom Out Reset image size

With estimated convective turnover times and Rossby numbers for our 15 M-type stars, we then used results from Saar & Brandenburg (1999), Lehtinen et al. (2016), and Olspert et al. (2018a) to gather a comparison sample of 40 F-, G-, and K-type stars with robust, well-defined cycles. Of these 40 stars, listed in Table 2, 10 have double cycles. Note, however, Olspert et al. (2018a) used multiple different methods to identify stellar cycles; for consistency with our approach, we use results obtained via their harmonic model (see Section 2.5 for more details). To estimate the convective turnover times of these stars, we use the method described above. ¹⁰ We note, however, that we exclude HD 18256 (F6) as it is too hot for our τ estimations; interestingly, in their study of Ca ii H and K cycles in 15 stars, Baliunas et al. (2004) found that this star had by the far the largest anharmonicity, which is a measure of the deviation from single-period sinusoidal behavior, in their sample.

4.1.1. Data Selection

Figure 14 shows a plot of cycle amplitude against period for the well-defined cycles in Table 1, along with a least-squares fit following ${A}_{\mathrm{cyc}}\propto {P}_{\mathrm{cyc}}^{0.94\pm 0.11}$. When fitting this power law, only stars with a single well-defined cycle per data set were used (i.e., when a star had a single cycle in both ASAS-3 and ASAS-SN data, then both cycles were used). The suspect 5.1 yr ASAS-4 cycle of Proxima was ignored during the fitting (see Section 2.3).

Zoom In Zoom Out Reset image size

As can be seen in Figure 14, cycles with longer periods generally also have larger amplitudes (ignoring stars with more than one concurrent cycle, whose cycles are connected with dotted lines and that are shown with faded markers). We speculate that this may be because longer periods allow a greter buildup of magnetic energy, resulting in a more pronounced cycle; this process is interrupted in stars with more than one concurrent cycle, and therefore they appear to contradict this trend. This relation may also explain why short-period cycles are detected less frequently in ASAS-3 data, which are noisier than their ASAS-SN successor. It may be the case that short-period cycles are present in ASAS-3 data, but cannot easily be resolved as their amplitudes are comparable to the ASAS-3 noise level.

Since stars with more than one well-defined cycle per data set were excluded from the power-law fit shown in Figure 14, it is interesting to note where these excluded cycles reside in relation to this fit. Indeed, it can be seen that for stars with two concurrent cycles, one circle is usually closer to this power law than the other; these other cycles may therefore be "false cycles" (Section 3.3) or indications of a multibranched relation, as seen by some in the P_cyc–P_rot relation (e.g., Saar & Brandenburg 1999; Böhm-Vitense 2007; Lehtinen et al. 2016). It should also be noted that A_cyc can be quite variable (more so in the Sun than in P_cyc), and so some outliers may be expected. That said, the (somewhat tentative) agreement among these excluded cycles is interesting, and we interpret this as further evidence for this relation between cycle amplitude and period.

Multiple branches in the A_cyc–P_cyc relation may explain why, for example, the ASAS-3 cycle of GJ 358 and ASAS-SN cycle of GJ 447 are outliers from the power law shown in Figure 14. As discussed in Section 3.2.2, however, GJ 358 already appears to be a highly unusual star. GJ 447, on the other hand, does have a second cycle period in Table 1, although it is not well constrained. This second cycle has an estimated period of 10.6 ± 6.9 yr and an estimated amplitude of 21.5 ± 1.0 mmag, which would put it very close to the power law shown in Figure 14, although with a large period uncertainty.

Using chromospheric Ca ii H and K data on FGK stars, Saar & Brandenburg (2002) found stronger evidence for multiple branches in the A_cyc–P_cyc relation owing to their larger sample of stars. However, as the amplitudes from Saar & Brandenburg (2002) are from chromospheric indicators, they are not directly comparable to our photometric amplitudes. It would be very interesting to see how the photometric cycle amplitudes of F–M type stars compare, but this is left to future work.

In the case of our Sun, the solar cycle does not have a strict 11 yr periodicity; indeed, the ∼4 year range of the solar P_cyc (Donahue & Baliunas 1992) is fully ±18% of the average value. Moreover, if a given solar cycle has a shorter than average period, then it will usually have a larger than average amplitude, and vice versa (related to the Waldmeier effect; Waldmeier 1935). For the M-type stars analyzed in this work, the observation intervals cover only a handful of cycles in the best case, so it is difficult to determine whether these cycles reflect a natural range in P_cyc as seen in the Sun. However, from Figure 14, it is clear that for M-type stars that exhibit two apparent cycles, the shorter cycle usually has a greater amplitude; in cases where the shorter P_cyc is approximately half the longer P_cyc, this is consistent with one polarity being stronger than the other (i.e., a Hale, rather than a Schawbe, cycle, as suggested by do J.-D. Nascimento et al. 2023, in preparation). Alternatively, if, as Olspert et al. (2018a) suggest, these double cycles are the result of a single quasiperiodic cycle, then perhaps this suggests that M dwarfs behave similarly to the Sun - with shorter than average cycles having greater than average amplitudes, and vice versa.

The plot of the cycle period against rotation period for the well-defined cycles in Table 1, along with previously detected cycles in FGK stars presented in Table 2, shows high scatter, without strong correlations. However, convective zones are believed to be a key ingredient of stellar cycles, and these are not accounted for in this plot. We therefore place little emphasis on this plot (and choose to not include it here) as it ignores a major physical parameter of current stellar cycle models: the Rossby number.

Figure 15 plots the cycle period against the Rossby number for the well-defined cycles presented in Table 1, along with the previously detected cycles in FGK stars presented in Table 2. This figure shows that when the convective turnover time is accounted for in the form of the Rossby number, M-type stars appear to have similar cycles to those of FGK stars at equivalent Ro. However, it is clear that our M dwarfs favor shorter cycle periods, extending the lower cluster of FGK stars to lower P_cyc. M dwarfs also isolate themselves at low Rossby numbers, although this may be a sampling effect; perhaps M dwarfs extend the upper cluster of FGK stars to lower Ro as well.

Zoom In Zoom Out Reset image size

In Figure 15, the longer cycle periods seen in FGK stars likely reflect the fact that these stars have generally been monitored over longer periods than M dwarfs. The evidence that M dwarfs exhibit shorter P_cyc (relative to FGK stars) is strong, but we also found indications of longer, currently unresolvable or poorly constrained cycles in 11 of our M dwarfs. Reproducing these figures with another 10 yr of ASAS-SN data may therefore prove interesting.

Figure 16 shows a plot of the cycle period divided by the rotation period against inverse Rossby number for the well-defined cycles presented in Table 1, along with the previously detected cycles in FGK stars presented in Table 2.

Zoom In Zoom Out Reset image size

The power-law fit to the M-dwarf data in Figure 16 shows that these data follow the relation

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{cyc}}}{{P}_{\mathrm{rot}}}=(3.6\pm 1.2)\times {{Ro}}^{-1.02\pm 0.06}.\end{eqnarray} \tag{ 4 }$$

Similarly, Figure 16 shows that FGK stars follow the relation

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{cyc}}}{{P}_{\mathrm{rot}}}=(63.6\pm 1.1)\times {{Ro}}^{-1.06\pm 0.09},\end{eqnarray} \tag{ 5 }$$

suggesting that P_cyc ∝ τ for all stars, in contrast to the results of Saar & Brandenburg (1999)—shown by the A and I branches (which we have refit using our updated sample).

Brandenburg et al. (1998) and Saar & Brandenburg (1999) interpreted the ratio of P_cyc/P_rot as a measure of the α effect, based on the simple αΩ dynamo model proposed by Robinson & Durney (1982) and Noyes et al. (1984). In this framework, Equations (4) and (5) suggest that M dwarfs have weaker α effects than FGK type stars at equivalent Rossby numbers (note that the scaling constant in Equation (4) is 17.7 times smaller than in Equation (5)). This result could be explained by considering the physical properties of these stars. Since M-type stars have deeper convective zones and consequently longer turnover times, longer rotation periods are needed to obtain similar Rossby numbers. Longer rotation periods therefore mean that the Coriolis forces that act within these convective zones are weaker, resulting in less helical turbulent convection. Furthermore, the average convective velocities themselves are lower in these cooler stars; both these properties act to reduce the α effect.

Note, however, that the above explanation would suggest a smooth transition from F to M in Figure 16. In contrast, this figure shows a gap between FGK and M stars. Perhaps, because M stars are (almost) fully convective and therefore have (almost) no tachocline, this causes a significant change in the α effect. Alternatively, this gap could simply be the result of our limited sample of stars, which includes two K5s stars and one K7 star, but then no stars until M3. Clearly, Figure 16 needs to be extended to cover a wider range of stars, with a wide range of physical properties and rotation rates, to better understand whether this gap is real - and if so, what may be causing it.

The axes of Figure 16 have the benefit of being dimensionless, and Ro in particular is directly related to the mean-field dynamo number, D, as mentioned in Section 4.1. Clusters and trends within the data can thus provide insights into the underlying dynamo properties. Indeed, we have already discussed the implications of the scaling constants of Equations (4) and (5); below, we discuss the implications of the power-law indices of these equations.

According to most theoretical models (e.g., Ruediger & Kichatinov 1993), the α effect should be quenched with growing magnetic fields, based on the reasoning that the stronger a star's magnetic field, the stronger the suppression of cyclonic convection (i.e., the α effect). This paradigm persisted until Saar & Brandenburg (1999) found evidence suggesting that the α effect may be antiquenched on their A and I branches. Here, with the benefit of a more focused sample of stars (no binary or evolved stars), we reinvestigate these claims.

Using a 2D nonlinear dynamo model, Tobias (1998) found

$$\begin{eqnarray}&&{P}_{\mathrm{cyc}}\propto {D}^{\gamma },\end{eqnarray} \tag{ 6 }$$

where D is the dynamo number, and −0.67 ≤ γ ≤ −0.38. To explore the role of quenching in a simpler linear 1D dynamo model, Saar & Brandenburg (1999) defined a quenching index, q, such that

$$\begin{eqnarray}&&D\propto \alpha {\rm{\Omega }}^{\prime} \propto {{Ro}}^{-q-2},\end{eqnarray} \tag{ 7 }$$

where ${\rm{\Omega }}^{\prime}$ is the radial DR, and

$$\begin{eqnarray}&&q={q}_{\alpha }+{q}_{{\rm{\Omega }}}-2,\end{eqnarray} \tag{ 8 }$$

where q_α is the quenching index of the α effect (i.e., $\alpha \propto {{\rm{\Omega }}}^{{q}_{\alpha }}$, with ${\rm{\Omega }}\equiv \tfrac{2\pi }{{P}_{\mathrm{rot}}}$), and q_Ω is the quenching index of the Ω effect (i.e., ${\rm{\Omega }}^{\prime} \propto {{\rm{\Omega }}}^{{q}_{{\rm{\Omega }}}}$). From the solution of the standard dynamo equations (assuming a fixed wavenumber, k), Brandenburg et al. (1998) noted that the cycle frequency was given by

$$\begin{eqnarray}&&{\omega }_{\mathrm{cyc}}\equiv \displaystyle \frac{2\pi }{{P}_{\mathrm{cyc}}}\propto {\left(\displaystyle \frac{\alpha {\rm{\Omega }}^{\prime} {kL}}{2}\right)}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 9 }$$

where L is the dynamo length scale. Combining Equations (7) and (9), we have

$$\begin{eqnarray}&&{P}_{\mathrm{cyc}}\propto {(\alpha {\rm{\Omega }}^{\prime} )}^{-\tfrac{1}{2}}\propto {{Ro}}^{\tfrac{q}{2}+1},\end{eqnarray} \tag{ 10 }$$

or

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{cyc}}}{{P}_{\mathrm{rot}}}\propto {{Ro}}^{\tfrac{q}{2}}.\end{eqnarray} \tag{ 11 }$$

Since Figure 16 suggests that the relation between P_cyc/P_rot-Ro can be described by a power law, where the power-law index can take various values, we can thus write

$$\begin{eqnarray}&&{{Ro}}^{\tfrac{q}{2}}\propto {{Ro}}^{\delta },\end{eqnarray} \tag{ 12 }$$

where δ is the power-law index (i.e., δ = −1.04 from Equation (4)). Thus we find (by combining Equations (8) and (12))

$$\begin{eqnarray}&&{q}_{\alpha }=2\delta -{q}_{{\rm{\Omega }}}+2.\end{eqnarray} \tag{ 13 }$$

With this equation, we can infer the scaling of the α effect using the power laws plotted in Figure 16. First, however, we must estimate q_Ω (the quenching index for radial DR).

DR measurements are difficult and fraught with misidentifications (Aigrain et al. 2015; Basri 2018), but because good evidence suggests that surface DR scales linearly with rotation period for slower rotators (e.g., Saar 2011), it is reasonable to assume q_Ω ≈ 1; for saturated activity stars (discussed in Section 4.7), the dependence reverses, and q_Ω ≈ −1.7 (Saar 2011). We note that we are implicitly assuming that radial DR scales linearly with the measured surface DR. ¹¹ With this assumption and using these q_Ω, we thus have

$$\begin{eqnarray}&&{q}_{\alpha }=2\delta +1,\end{eqnarray} \tag{ 14 }$$

in the unsaturated regime, and

$$\begin{eqnarray}&&{q}_{\alpha }=2\delta +3.7,\end{eqnarray} \tag{ 15 }$$

in the saturated regime. Using the power-law indices from Equations (4) and (5), we find q_α,M = −1.03 and q_α,FGK = −1.13 in the unsaturated regime, and q_α,M = 1.67 and q_α,FGK = 1.57 in the saturated regime. This suggests that the α effect is quenched with faster rotation in the unsaturated regime (but antiquenched in the saturated regime) similarly for FKG and M stars.

As mentioned previously, quenching of the α effect (in the unsaturated regime) is an expected result. In contrast, the α effect antiquenching that we find in the saturated regime is surprising - although antiquenching is supported by some magnetohydrodynamic (MHD) simulations (e.g., Chatterjee et al. 2011). We note that our findings in the unsaturated regime disagree with the results of Saar & Brandenburg (1999), whose limited sample of stars (and different assumptions about Ω effect scaling) led to the opposite conclusion for their A and I branches (note that the fit in Figure 16 is more consistent with their tentative transitional branch for faster rotators; see Lehtinen et al. 2016). That said, it should be noted that these results are based on a very simple dynamo model (Brandenburg et al. 1998), along with some debated assumptions about DR and α and Ω effect quenching (e.g., Barnes et al. 2005; Saar 2011).

The clustering of stars at the end of the FGK sequence in Figure 16 (particularly around the putative I branch), together with the spread of the sequence in Ro, suggests that perhaps the concept of A and I branches is not completely without merit. If we consider possible A and I branches, however, we must note that we have not been self-consistent in our use of τ. If the I branch is dominated by a tachocline dynamo (as assumed for the Sun), a local τ, τ_L , specifically computed for that location, is more suitable. Further, a fit to all the FGK stars together becomes inappropriate; the faster rotators (using τ_G , as before) and the slower I branch stars (using τ_L ) should be treated separately. Using this distinction, we present the results in Figure 17.

Zoom In Zoom Out Reset image size

Due to reasons that will become clear in Section 4.7, for stars with cycles on both the A and I branches, we deem the I branch cycle to be the primary and the A branch cycle the secondary. In a preliminary version of Figure 17, secondary cycles were plotted also using τ_L , but we noted that this caused larger scatter in the A branch fit (σ = 1.07 in log space). If we use τ_G for these secondary cycles, effectively assuming they arise from a convection zone dynamo still in operation, the A branch scatter is reduced by ∼20% (σ = 0.81; Figure 17). These new values modify Equation (5) to

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{\mathrm{cyc}}}{{P}_{\mathrm{rot}}}=(162.6\pm 0.2)\times {{Ro}}^{-0.81\pm 0.17},\end{eqnarray} \tag{ 16 }$$

but do little to change the spirit of the discussion in Section 4.5. A list of different model parameters, their goodness-of-fit metrics, and their implied q_α values is provided in Table 3.

4.6.1. GJ 358

Indeed, our results here and those of Saar & Brandenburg (1999) may not necessarily be at odds; perhaps we are seeing a switch in dominance between two different types of dynamos with different properties. We largely set aside the discussion of the A branch here because more recent work casts doubt on its robustness/reality (e.g., Boro Saikia et al. 2018). We also skip discussion of the superactive (S) branch, which is occupied almost solely by binaries in which tidal rotational locking may alter the dynamo fundamentally and place it beyond the scope of the presumed single-star dynamos considered here.

We propose the following scenario. The most rapid rotators seem to have nearly solid-body rotation, for example, GJ 1243 (M4, P_rot = 0.59 d, k/k_⊙ = (ΔΩ/Ω)/(ΔΩ_⊙/Ω_⊙) ≈ 0.01; Davenport et al. 2020), LQ Hya (K2, P_rot = 1.597 days, k/k_⊙ = 0.028; Kovári et al. 2004), or V889 Her (G2, P_rot = 1.337 d, k/k_⊙ = 0.045; Kővári et al. 2011). With so little DR, their dynamos are likely a turbulent α² type or at best α²Ω, with the dynamo extant throughout the convection zone. This concept is consistent with models (e.g., Guerrero et al. 2019), and some detailed MHD models obtain roughly the same power-law index for P_cyc/P_rot ∝ Ro^δ : for example, Warnecke (2018) obtain −0.89 ± 0.04, and if one does a least-squares fit of the five long cycle models in Brun et al. (2022; which are clearly not on the I branch) using their stellar Ro, the slope is −0.90 ± 0.47. These results both compare well with our values of −1.06 ± 0.09 and −0.81 ± 0.17.

We therefore identify the fits found here as representing the initial, more turbulent, α-antiquenched, full convection zone dynamos. Stars on the ZAMS start with rapid rotation, low DR and saturated magnetic activity, likely with predominantly poloidal large-scale field structure (Donati & Landstreet 2009). As stars lose angular momentum, they evolve from the upper right to the lower left in Figure 16, with P_cyc/P_rot decreasing, and (if Saar 2011 is correct) increasing DR (sketched in Figure 18; red arrows). DR increases to a maximum point, at which point, saturated magnetic activity stops, the α effect starts becoming quenched, and further spindown leads to activity reduction along the well-known rotation-activity relations.

Zoom In Zoom Out Reset image size

(As a brief aside, we note that an interesting alternative to the full convection zone dynamo, particularly in M stars, is the idea of dynamo cycles driven by thermal fluctuations (e.g., Nigro 2022). This could potentially explain the lack of rotational dependence in M stars, as the distribution of thermal fluctuations would primarily depend on mass/T_eff, echoed in the trend seen here of P_cyc ∝ τ. The similar trend in faster-rotating FGK stars could also be influenced, at least in part, by such a thermally driven dynamo. More work on these models would be useful.)

At some point in partially convective stars, the radiative core and convective envelope decouple (e.g., MacGregor & Brenner 1991), allowing a region of strong shear to form at the interface (the tachocline). The Sun and slower FGK rotators are thought to have αΩ dynamos with the primary generation region in this strong shear layer. M stars, lacking significant (or any) tachoclines, cannot make this transition.

However, once the tachocline dynamo dominates in older FGK stars, which happens at least by the time the stars reach the I branch around Ro⁻¹ = 1–3, further cycle evolution proceeds along the I branch, with P_cyc/P_rot slightly increasing, or may be remaining roughly constant (shown by a blue arrow in Figure 18). Perhaps there is a transitional phase of mixed dynamo dominance before this, explaining the diffuse A branch and the abundance of stars with multiple cycles in this region. The recently noted bend in the chromospheric rotation-activity relation (Lehtinen et al. 2021) may be another reflection of this change. ZDI shows that large-scale toroidal fields already dominate well before the I branch is reached, at around Ro⁻¹ ≈ 15 (on the τ_C scale used here, converted from Donati & Landstreet 2009), about where DR peaks and activity saturation ends (shown by the right dashed magenta line in Figure 18). Perhaps this is where core/envelope decoupling also takes place. Large-scale toroidal field dominance only ends at Ro⁻¹ ≈ 2Ro ${}_{\odot }^{-1}\approx 1$ (Donati & Landstreet 2009)—i.e., at about the time when the star first reaches the I branch. Note that no further stars are seen along the −1.06 slope in Figure 16 among the FGK stars once the I branch is reached. We view this as supporting the final dominance of the tachocline dynamo.

Following the arguments using the same simple dynamo model, Saar & Brandenburg (1999) suggest that the change in slope at the I branch is due to α antiquenching, which implies (surprisingly) that here α increases with the magnetic field (q_α = 1.10 on the I branch in Figure 17; Table 3). The reason might be that α in this regime is driven by magnetic buoyancy and pumping effects (Chatterjee et al. 2011). Here again, some detailed MHD models are in reasonably good agreement: a power-law fit to the results of models M1–4 of Warnecke (2018; their slow-moderate rotators, ignoring the very slowest) yields P_cyc/P_rot = 213Ro^0.514, which lies between the A and I branches and is of similar slope (shown by the dotted black line in Figure 18). As noted above, we further suggest that some secondary cycles are remnant convection zone dynamos still functioning as the tachocline dynamo ramps up to dominance. This dual P_cyc—dual dynamo transitional regime corresponds to the purple transitional epoch arrow in the dynamo evolution in Figure 18.

In the end, however, after further evolution along the I branch, below a critical Ro⁻¹, the dynamo may begin to sputter (van Saders et al. 2016), the star begins to experience magnetic grand minima (Metcalfe et al. 2016), differential rotation weakens (Metcalfe et al. 2022) or may even reverse sign (Gastine et al. 2014; Karak et al. 2020), and the dynamo dies, with P_cyc going to infinity (see, e.g., the lengthening P_cyc in the evolutionary sequence of G2 stars 18 Sco, the Sun, and α Cen A; Judge et al. 2017). With the cycling dynamo shut down, large-scale fields driving spindown vanish, stopping further rotational evolution, and the star evolves vertically off the diagram at constant Ro (e.g., Metcalfe & van Saders 2017, shown by a vertical gold arrow in Figure 18).