On the Origin of the Dichotomy of Stellar Activity Cycles

Despite the enormous progress made in recent years in the observations of stellar magnetism, its origin and generation are still not completely understood. It is believed that the large-scale magnetic field observed in cool dwarfs can be explained by the idea of "cyclonic" turbulence (Parker 1955): in a nutshell, a toroidal field produced from a poloidal one by the shear (Ω effect) present in a stellar convection zone further decays in helical turbulence, which (α effect) can then reinforce the dipole field (Steenbeck et al. 1966). The magnetic field produced by the αΩ dynamo mechanism is in general nonstationary and its late-time behavior is characterized by a power law of the type

$$\begin{eqnarray}&&{\omega }_{\mathrm{cyc}}\propto {{\rm{\Omega }}}^{\nu },\end{eqnarray} \tag{ 1 }$$

where ω_cyc ≡ 2π/P_cyc is the rate of the activity cycle, Ω ≡ 2π/P_rot is the rotation rate, and ν is a scaling exponent. During the unstable phase, assuming ${\rm{\Omega }}\propto {\rm{\Omega }}^{\prime} {\ell }$, ${\rm{\Omega }}^{\prime}$ being the shear and ℓ the characteristic scale of the turbulence, it turns out that ν = 1/2 (Stix 1976), while in the opposite limit, in the saturated phase, the cycle period is fixed by the characteristic turbulent diffusive timescale and ν = 0. For this reason it is reasonable to expect that in general ν will lay between these two limits, at least according to kinematic dynamo theory. At the nonlinear level this picture can drastically change and it is not clear if a simple scaling law like Equation (1) holds at all (Warnecke 2018; Guerrero et al. 2019; Pipin 2021).

Pioneering observational studies devoted to finding possible correlation between ω_cyc and Ω mostly agreed on ν ≈ 1 (Wilson 1978; Noyes et al. 1984; Baliunas et al. 1996) albeit the presence of multiple scaling laws (branches) has further complicated the discussion (Saar & Brandenburg 1999; Saar 2002; Böhm-Vitense 2007; Hall et al. 2009; Metcalfe et al. 2016; but see also Brun & Browning 2017 for an extended review).

Stars can be classified in active or nonactive according to the offset of two nearly identical scaling laws in the ${\mathrm{log}}_{10}({\omega }_{\mathrm{cyc}}/{\rm{\Omega }})$–${\mathrm{log}}_{10}{R}_{\mathrm{HK}}^{{\prime} }$ plane: the so-called RCRA (ratio of cycle over rotation versus activity) diagram. This classification plays a fundamental role in the analysis of stellar activity cycles because it determines the existence of possible correlations in the P_cyc–P_rot plane (Boro Saikia et al. 2018).

Note that this has nothing to do with the existence of the Vaughan–Preston gap (Vaughan & Preston 1980); a mild depletion in the number of stars around ${\mathrm{log}}_{10}{R}_{\mathrm{HK}}^{{\prime} }\,=\,-4.75$ caused by the weak bimodality of the distribution of ${\mathrm{log}}_{10}{R}_{\mathrm{HK}}^{{\prime} }$ in main-sequence stars (Gomes da Silva et al. 2021).

The question has been further addressed in Olspert et al. (2018) where a Gaussian classification algorithm has been used to determine the presence of multiple populations in the RCRA diagram. The problem with a direct clustering analysis in the RCRA diagram is the existence of possible correlations between ω_cyc and Ω that would render the interpretation of the results ambiguous.

In this work we follow a different route. In Section 3 we search for scaling laws of the type given by Equation (1) in the ${\mathrm{log}}_{10}{\omega }_{\mathrm{cyc}}$–${\mathrm{log}}_{10}{\rm{\Omega }}$ plane by considering randomly distributed couples of lines in this plane. We thus employ a Bayesian classification algorithm based on a generalized likelihood describing a Gaussian mixture of two populations, so that the existence of branches is determined by the distance from the lines in this plane and not by the cluster center as in the approach followed by Olspert et al. (2018). We think that in the search for possible correlations between P_cyc and P_rot, the ${\mathrm{log}}_{10}{\omega }_{\mathrm{cyc}}$–${\mathrm{log}}_{10}{\rm{\Omega }}$ plane has some advantages compared with the P_cyc–P_rot plane used by some authors (e.g., Böhm-Vitense 2007; Metcalfe et al. 2016; Boro Saikia et al. 2018) because it allows for a direct probe of a large class of power-law dependencies as the one in Equation (1). Nonetheless, the data could be better reproduced by more complicated nonlinear laws: one has to keep in mind that Equation (1) is only a suggestion inspired by kinematic dynamo theory.

We present the data sample in Section 2 and the data analysis in Section 3. We show our results in Section 4 and discuss them in the light of available stellar properties in Section 5. Finally, we draw our conclusions on the origin of the two populations in Section 6.

The analysis proposed in this work is based on measurements of rotation period, P_rot and magnetic activity cycle period, P_cyc that were published by Boro Saikia et al. (2018) (hereafter B18) and Olspert et al. (2018) (hereafter O18). The adopted cycle periods are obtained from chromospheric Ca, H, and K measurements and are therefore excluding fast rotators, which are typically very young objects (e.g., Lehtinen et al. 2016; Oláh et al. 2016; Distefano et al. 2017). The catalogs by B18 and O18 contain a substantial intersection of stars, for which multiple measurements are thus available. We therefore compiled the largest catalog of activity cycle measurements available to date (67 stars, including the Sun), by taking into account all stars in common between B18 and O18 (31 in total, including the Sun), as well as stars that belong to B18 only (14 stars) and to O18 only (22 stars). While for stars that are not common between the two catalogs the merging procedure is straightforward because only one set of measurements is available, for the stars falling in the intersection we applied the procedure described in the following. O18 provides three estimates of P_cyc: one based on a periodic model (P_cyc,P), one on a harmonic model (P_cyc,H), and another one on a quasi-periodic model (P_cyc,QP). When possible we privilege the estimates obtained by means of a quasi-periodic model because (i) it is the most commonly available and (ii) it provides a more accurate fit to the time series thanks to its higher complexity. For obtaining the final estimates of P_cyc we computed a weighted mean between the estimates from O18 and those from B18, i.e., taking into account the uncertainty in P_cyc, except for the stars that we discuss below. HD 100180 has no P_cyc measurements from O18 so we considered the one from B18. HD 101501 was removed because the P_cyc measurements between the two catalogs deviate by more than 30%. HD 149661 has only P_cyc,P from O18, while B18 provides two different estimates. Here we considered as a final P_cyc the weighted mean between O18 and the estimate from B18 that is closer to that of O18. HD 156026 has P_cyc,H and P_cyc,P from O18, which are almost identical, and an estimate from B18 that is also very close. We considered the weighted mean between P_cyc,P and the value from B18 as a final estimate. HD 155886 has P_cyc,H and P_cyc,QP from O18, and an estimate from B18, which are all substantially different from one another. We removed this star to avoid ambiguities. HD 1835 and HD 20630 have no P_cyc measurement from O18 and they have each two different estimates from B18. We removed these targets to avoid ambiguities. HD 190007 has P_cyc,P and P_cyc,QP from O18, which are almost identical, while a substantially different one comes from B18. We considered P_cyc,QP from O18 as the final one. HD 190406 has P_cyc,QP from O18, which is very close to one of the two estimates provided by B18. For the final value, we computed the weighted mean of the two. Despite the fact that HD 201092 has P_cyc,P from O18, which is very close to the value provided by B18, we discarded this star from our sample as its determination of the activity cycle period appears rather uncertain due to the pronounced caothicity of its activity indicator time series (Hempelmann et al. 2006; Robrade et al. 2012). HD 76151 and HD 78366 both have only P_cyc,H from O18 and an estimate from B18 that is very close. The final estimate is taken as a weighted mean of the two values. HD 82443 has P_cyc,P and P_cyc,QP from O18, which are not in agreement, and another value substantially different from these two from B18. We removed this target to avoid ambiguities. The final uncertainties on P_cyc are computed following a standard error propagation from those listed by the two catalogs.

In relation to P_rot, to obtain a final estimate we computed a simple arithmetic mean between the two catalog estimates, except for the stars HD 26923 and HD 32147, which were excluded from our final catalog because of the pronounced disagreement between the two measurements (>30%). Moreover, we computed values for the chromospheric activity index ${\mathrm{log}}_{10}{R}_{\mathrm{HK}}^{{\prime} }$ by taking the arithmetic mean of the measurements available from O18 and B18 for those stars that are in common. For this parameter we find an excellent agreement (well within 10%) between the two catalogs.

By relying on a calibration of the Rossby number (Ro) provided by Corsaro et al. (2021) using the (B − V) color index, we find that our stellar sample has a mean Ro = 0.523 ± 0.481. This shows that, although our catalog comprises stars spanning from K to F spectral types and from main-sequence to early giants, their Rossby number is on average close to that found in the Sun (Ro_⊙ = 0.496), meaning that it is reasonable to study the sample as a whole.

2.1. Stellar Metallicity and Age

The period measurements were complemented with additional measurements of stellar metallicity and age by crossmatching our compiled catalog with the newest available ESA Gaia (Gaia Collaboration et al. 2016) DR3 catalog of additional astrophysical parameters from the Apsis processing chain (paramsup; Fouesneau et al. 2022), and with high-resolution spectroscopy measurements.

The Gaia DR3 catalog provides metallicities $\left[{\rm{M/H}}\right]$ from low-resolution spectra for 47 stars of our sample. For the high-resolution measurements instead, we considered literature values from works exploiting a consistent data analysis procedure, as obtained by Valenti & Fischer (2005) (31 stars), Brewer et al. (2016) (20 stars), Luck (2015) (6 stars), and Luck (2017) (39 stars). Multiple high-resolution $\left[{\rm{M/H}}\right]$ measurements for an individual star, when available, are combined through an arithmetic mean.

Stellar ages are obtained to a large extent from the Gaia DR3 catalog, which provides isochrone-based estimates (related to BaSTI stellar evolution models; Hidalgo et al. 2018) for 33 stars of our sample. For HD 146233 (18 Sco) we replaced the Gaia DR3 estimate with an age based on asteroseismology (Li et al. 2012). Additional measurements are taken from the StarHorse Gaia EDR3 catalog of astrophysical parameters (Anders et al. 2022). In this catalog we found ages for the stars HD 152391, HD 166620, HD 219834A, HD 219834B, HD 26913, HD 21749, HD 154577, HD 75332, and HD 88373, which were instead missing in the Apsis Gaia DR3 catalog. We complete our set of ages by including measurements for the Sun (Bonanno and Fröhlich 2015), and the well-known stars HD 124897 (Arcturus; Ramírez, & Allende Prieto 2011) and HD 201091 (61 Cyg A; Kervella et al. 2008), thus leading to a total of 45 stars having an age estimate. All of the collected characteristics for the ensemble are provided in Table 1 in a machine-readable format.

Table 1. Metadata for Compiled Catalog

Unit	Label	Description
⋯	Name	Star name, HD catalog, or Sun
⋯	DR3	Gaia DR3 identifier
Days	Prot	Rotational period, days
Yr	Pcyc	Activity cycle, years
Yr	e_Pcyc	Standard deviation of Pcyc
Dex	[M/H]–LR	Gaia DR3 low-resolution-based metallicity
Dex	[M/H]–HR	High-resolution (literature) metallicity
Dex	[M/H]	Final metallicity
Gyr	Age	Age, gigayears
Dex	logRHK	Chromospheric activity index, ${\mathrm{log}}_{10}{R}_{\mathrm{HK}}^{{\prime} }$
M⊙	Mstar	Mass of star
Dex	logg	log, surface gravity, $\mathrm{log}g$
L⊙	Lstar	Luminosity of star
Mag	B − V	B − V color
Days	τ	Convective turnover time

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

3.1. Models

As originally presented by Saar & Brandenburg (1999) and then further discussed by, e.g., Brandenburg et al. (2017) (see also O18), the magnetic activity cycles are distributed along two potentially different regimes that have been termed active and inactive branches. In this work we focus on the observable quantities ω_cyc and Ω, and following Equation (1) we assume that these two observables are related by a simple power-law relation of the type

$$\begin{eqnarray}&&{\omega }_{\mathrm{cyc}}=\beta {{\rm{\Omega }}}^{\nu }.\end{eqnarray} \tag{ 2 }$$

We therefore set our analysis in the ${\mathrm{log}}_{10}{\omega }_{\mathrm{cyc}}$–${\mathrm{log}}_{10}{\rm{\Omega }}$ plane, hereafter OCOR (omega cycle versus omega rotation) plane for simplicity, where the power-law relation can be linearized. Following the literature, first we consider two regimes that can be modeled through the presence of two simultaneous linear relations of the type

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\omega }_{\mathrm{cyc}}={\nu }_{1}{\mathrm{log}}_{10}{\rm{\Omega }}+{\mathrm{log}}_{10}{\beta }_{1}\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\omega }_{\mathrm{cyc}}={\nu }_{2}{\mathrm{log}}_{10}{\rm{\Omega }}+{\mathrm{log}}_{10}{\beta }_{2},\end{eqnarray} \tag{ 4 }$$

which we term model ${{ \mathcal M }}_{1}$. In the second case we consider

$$\begin{eqnarray}&&{\mathrm{log}}_{10}{\omega }_{\mathrm{cyc}}={\nu }_{0}{\mathrm{log}}_{10}{\rm{\Omega }}+{\mathrm{log}}_{10}{\beta }_{0},\end{eqnarray} \tag{ 5 }$$

which we term model ${{ \mathcal M }}_{2}$. This second model has the purpose of allowing us to verify whether a single power law, hence a single regime in the evolution of the magnetic activity cycle, ought to be preferred over two separate ones.

In all cases here presented, the coefficients (ν_i , β_i ) for i = 0, 1, 2, will be estimated through a statistical inference (see Section 3.2 for details). Besides, a model comparison process will allow us to draw conclusions on whether one or two regimes ought to be preferred in the light of the existing data (Section 3.3).

3.2. Likelihood

For the simultaneous fitting of two regimes as presented through model ${{ \mathcal M }}_{1}$ we do not perform an a priori division of the data sample in two parts, which would otherwise result in imposing a somewhat arbitrary cut. We instead automatically identify the presence of any internal partition by means of a purely statistical approach. For this purpose we consider a Gaussian mixture likelihood, consisting of two components of the type ${y}^{(1)}={\nu }_{1}x+{\mathrm{log}}_{10}{\beta }_{1}$, ${y}^{(2)}={\nu }_{2}x+{\mathrm{log}}_{10}{\beta }_{2}$, with ν₁, ν₂ the two line slopes, and ${\mathrm{log}}_{10}{\beta }_{1},{\mathrm{log}}_{10}{\beta }_{2}$ the two line offsets. The corresponding likelihood can be expressed as

$$\begin{eqnarray}&&{ \mathcal L }\left({\boldsymbol{\theta }}\right)=\prod _{i=1}^{N}\left[{A}_{i}\left({\boldsymbol{\theta }}\right)+{B}_{i}\left({\boldsymbol{\theta }}\right)\right],\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{A}_{i}\left({\boldsymbol{\theta }}\right)=\displaystyle \frac{{\sin }^{2}\psi }{\sqrt{2\pi }{\sigma }_{i}}\exp \left[-\displaystyle \frac{{\left({y}_{i}-{\nu }_{1}{x}_{i}-{\mathrm{log}}_{10}{\beta }_{1}\right)}^{2}}{2{\sigma }_{i}^{2}}\right]\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{B}_{i}\left({\boldsymbol{\theta }}\right)=\displaystyle \frac{{\cos }^{2}\psi }{\sqrt{2\pi }{\sigma }_{i}}\exp \left[-\displaystyle \frac{{\left({y}_{i}-{\nu }_{2}{x}_{i}-{\mathrm{log}}_{10}{\beta }_{2}\right)}^{2}}{2{\sigma }_{i}^{2}}\right]\end{eqnarray} \tag{ 8 }$$

are the two components of the Gaussian mixture, with θ = (ν₁, β₁, ν₂, β₂, ψ) the parameter vector containing the free parameters that need to be fitted, and σ_i the uncertainties on the dependent variable y_i . We therefore have to estimate five free parameters in the case of two simultaneous regimes; ψ being the parameter that controls the relative weight of each of the two components of the mixture. In the way we presented the models, we further assume that the measurements ω_cyc are log-normally distributed (e.g., see Corsaro et al. 2013; Bonanno et al. 2014). This is because the residuals arising from the difference between the observed and predicted logarithms of the activity cycle rate are assumed to be Gaussian distributed.

For avoiding numerical overflows, which may arise from the exponential term in the likelihood, we switch to the log-likelihood and rely on the identity

$$\begin{eqnarray}\begin{array}{rcl} & & \mathrm{ln}(A+B)=\mathrm{ln}\left[\exp \left(\mathrm{ln}A-\mathrm{ln}B\right)+1\right]+\mathrm{ln}B.\end{array}\end{eqnarray} \tag{ 9 }$$

Therefore, we compute each term i using two different expressions. By defining the quantity

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{ij}}\equiv -\frac{1}{2{\sigma }_{i}^{2}}{\left[{y}_{i}-{y}_{i}^{(j)}\right]}^{2},\end{eqnarray} \tag{ 10 }$$

we obtain for Δ_i1 ≥ Δ_i2 ,

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Lambda }}}_{i}^{(a)}\left({\boldsymbol{\theta }}\right) & = & \mathrm{ln}\left(\displaystyle \frac{{\sin }^{2}\psi }{\sqrt{2\pi }{\sigma }_{i}}\right)+{{\rm{\Delta }}}_{{\rm{i}}1}\\ & & +\,\mathrm{ln}\left[1+{\cot }^{2}\psi \exp \left({{\rm{\Delta }}}_{{\rm{i}}2}-{{\rm{\Delta }}}_{{\rm{i}}1}\right)\right],\end{array}\end{eqnarray} \tag{ 11 }$$

and vice versa for Δ_i1 < Δ_i2

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Lambda }}}_{i}^{(b)}\left({\boldsymbol{\theta }}\right) & = & \mathrm{ln}\left(\displaystyle \frac{{\cos }^{2}\psi }{\sqrt{2\pi }{\sigma }_{i}}\right)+{{\rm{\Delta }}}_{{\rm{i}}2}\\ & & +\,\mathrm{ln}\left[1+{\tan }^{2}\psi \exp \left({{\rm{\Delta }}}_{{\rm{i}}1}-{{\rm{\Delta }}}_{{\rm{i}}2}\right)\right].\end{array}\end{eqnarray} \tag{ 12 }$$

In this way we can compute a final log-likelihood as

$$\begin{eqnarray}&&{\rm{\Lambda }}\left({\boldsymbol{\theta }}\right)=\sum _{i=1}^{N}{{\rm{\Lambda }}}_{i}^{(a,b)}\left({\boldsymbol{\theta }}\right),\end{eqnarray} \tag{ 13 }$$

where Λ^(a,b) can be either Λ^(a) or Λ^(b) depending on the conditions defined above.

For treating the case of model ${{ \mathcal M }}_{2}$, consisting of a single linear relation, we instead rely on a standard Gaussian likelihood. We note that the Gaussian mixture likelihood here presented and often adopted in Bayesian cluster analysis is a more generalized form of the standard Gaussian likelihood. The two likelihood functions become identical by imposing ψ = 0 or π/2, i.e., by assuming that the sample is not internally split in two separate subsets.

3.3. Bayesian Inference and Model Comparison

By means of the log-likelihood function presented in Section 3.2, we perform a Bayesian inference by adopting uniform prior distributions for each free parameter involved in the fitting process. We note that although the free parameter ψ is not directly appearing in the model equations, it is an important fitting parameter because it controls how the data sample is clustered in the light of the two components of the mixture. Given that ψ has to satisfy the identity ${\sin }^{2}\psi +{\cos }^{2}\psi =1$, we have adopted a uniform prior in the range [0, π/2] to allow for all the possible outcomes. We further point out that for considering more reliable uncertainties on ω_cyc than those provided in the literature, which are, in some cases, unrealistically more precise than what can be obtained for our Sun, we have inflated them for our analysis by incorporating a systematic term due to the intrinsically stochastic character of the dynamo mechanism (Mininni et al.2000). We can estimate this term by assuming it is of the same order of the dispersion of the solar activity cycle period. For computing a reliable estimate, we consider the longest available set of solar activity cycle maxima, as obtained from sunspot time series covering a total of 18 different activity cycles (in the years from 1818 to 2014). ¹ The result yields ${\sigma }_{{P}_{\mathrm{cyc},\odot }}^{\mathrm{sys}}=0.24\,$ yr.

To perform the Bayesian parameter estimation we make use of the public tool Diamonds ² (Corsaro & De Ridder 2014), which exploits a nested sampling Monte Carlo algorithm (Skilling 2004). The results of our analysis are presented in Section 4. In addition to the parameter estimation problem we perform a Bayesian model comparison through the odds ratio, which for two competing models ${{ \mathcal M }}_{i}$ and ${{ \mathcal M }}_{j}$ is defined as

$$\begin{eqnarray}&&{{ \mathcal O }}_{{ij}}=\displaystyle \frac{{{ \mathcal E }}_{i}}{{{ \mathcal E }}_{j}}\displaystyle \frac{\pi \left({{ \mathcal M }}_{i}\right)}{\pi \left({{ \mathcal M }}_{j}\right)}.\end{eqnarray} \tag{ 14 }$$

In our application, we can discard the term containing the ratio of the model priors $\pi \left({{ \mathcal M }}_{i}\right)$, which we assume to be the same for each of the models considered, so that the odds ratio is entirely based on the evaluation of the Bayesian evidence ${{ \mathcal E }}_{i}$. We refer to this odds ratio as the Bayes factor, which we can express in logarithmic form as $\mathrm{ln}{{ \mathcal B }}_{{ij}}=\mathrm{ln}{{ \mathcal E }}_{i}-\mathrm{ln}{{ \mathcal E }}_{j}$. According to Jeffreys' scale of strength (Trotta 2008), we are therefore able to understand whether a given model could be favored or not against its competitor. The Bayesian evidence is computed for the pair of models $({{ \mathcal M }}_{1}$, ${{ \mathcal M }}_{2})$ and is provided as a direct output of the computation by Diamonds. We note that the evaluation of Equation (14) is not based on any approximation such as that of the Bayesian Information Criterion, but on the full numerical calculation of the Bayesian evidence as obtained by Diamonds from the sampling of the multidimensional posterior probability distributions.

The output of the model comparison process yields a clear preference for the model that incorporates two different regimes in the OCOR plane, namely favoring model ${{ \mathcal M }}_{1}$ with $\mathrm{ln}{{ \mathcal B }}_{\mathrm{1,2}}\simeq 2.6\cdot {10}^{4}$, largely beyond a strong evidence condition. The result of the fit is shown in Figure 1, while the estimates of the free parameters of model ${{ \mathcal M }}_{1}$ are ${\overline{\nu }}_{1}={0.109}_{-0.005}^{+0.005}$, ${\mathrm{log}}_{10}{\overline{\beta }}_{1}=-{2.417}_{-0.004}^{+0.004}$, ${\overline{\nu }}_{2}=-{0.108}_{-0.008}^{+0.007}$, and ${\mathrm{log}}_{10}{\overline{\beta }}_{2}\,=-{2.903}_{-0.005}^{+0.005}$. We obtain a weight parameter $\overline{\psi }={0.69}_{-0.06}^{+0.07}$, showing that in the mixture both groups are playing equally important roles (${\sin }^{2}\overline{\psi }\simeq 0.49$).

Zoom In Zoom Out Reset image size

We classify the stars in two groups according to their relative differences in ${\mathrm{log}}_{10}{\omega }_{\mathrm{cyc}}$ with respect to each line of the fit, and identify stars having a larger ω_cyc as group 1 (G1) stars, and stars having a lower ω_cyc as group 2 (G2) stars. The Sun clearly belongs to G2 as it is well in line with the regime indicated by Equation (4).

In order to understand the potential origin of the dichotomy shown in Figure 1 we analyze it in relation to the fundamental stellar properties of metallicity $\left[{\rm{M/H}}\right]$, age, and chromospheric activity index $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ for each of the two groups identified.

In relation to metallicity, we find a significant correlation (Pearson's correlation coefficient ρ = 0.78) between the low-resolution values obtained from Gaia DR3 and the high-resolution ones, for the stars that are in common (a total of 33). We therefore combine the two sets by first shifting the zero-point of the Gaia $\left[{\rm{M/H}}\right]$ set to match that of the high-resolution $\left[{\rm{M/H}}\right]$ set (−0.185 dex as obtained from a linear regression over the two sets), and by subsequently computing an arithmetic mean between the two. The final metallicity sample therefore represents stars having both Gaia DR3 and high-resolution measurements available. The presence of an offset in the Gaia metallicities may have an impact on the adopted stellar ages. We therefore caution the reader that this represents a potential limitation of our work in relation to the interpretation of the age distributions.

To highlight any difference between groups 1 and 2 we additionally show their corresponding distributions as a function of stellar properties in Figure 2. We find that the stars in G2 are about 40% more metal-rich than the stars in G1, with median metallicity values ${\left[{\rm{M/H}}\right]}_{{\rm{G}}1}\simeq -0.09\,$ dex and ${\left[{\rm{M/H}}\right]}_{{\rm{G}}2}\simeq +0.05\,$ dex. According to a Kolmogorov–Smirnov test, the two metallicity distributions appear different to a p-value level of 0.015, close to the usual limit adopted for a strong detection (which is set to 0.003). Interestingly, we find that the stars in G2 are about 25% older than those in G1, having median values age _G1 ≃ 2.8 Gyr and age _G2 ≃ 3.5 Gyr. Although the age difference between the two groups cannot be considered statistically significant (the associated p value of the two distributions is 0.17); the trend appears in line with that found for metallicity because stars with lower metallicity have a reduced opacity, which leads to stronger radiative fluxes, hence a faster evolution. For the sake of completeness, we also performed a comparison with ages computed from gyrochronology (Barnes 2007). Here we note that while ages from gyrochronology are in agreement with those from isochrones for G1 (where we obtain a median gyroage of 2.6 Gyr), there is a discrepancy for stars in G2, with gyroages pointing to a median of 1.6 Gyr. According to gyrochronology, this would imply that stars in G2 should rotate slower than observed to match the isochrone ages. While the investigation of this difference is beyond the scope of this work, we cannot exclude that stars in G2 may be subject to a different gyrochronological calibration than the one presented by Barnes (2007).

Zoom In Zoom Out Reset image size

In relation to the activity index ${\mathrm{log}}_{10}{R}_{\mathrm{HK}}^{{\prime} }$ instead, we find two rather similar median values, ${\mathrm{log}}_{10}{R}_{\mathrm{HK},\ {\rm{G}}1}^{{\prime} }\simeq -4.75\,$ dex and ${\mathrm{log}}_{10}{R}_{\mathrm{HK},\ {\rm{G}}2}^{{\prime} }\simeq -4.69\,$ dex. Here one could in principle attempt to rely on an age–activity relationship such as the one proposed by Mamajek & Hillenbrand (2008) to estimate an age difference between the two groups from the chromospheric activity level alone. However, we point out that according to a Kolmogorov–Smirnov test the two distributions in ${\mathrm{log}}_{10}{R}_{\mathrm{HK}}^{{\prime} }$ are essentially indistinguishable (p value > 0.9). This implies that the variance in chromospheric activity between G1 and G2 is just too small to be reliably used to infer any dissimilarity in age. Besides, the general differences identified by our computation of the median estimators appear clearly visible in the distributions shown in Figure 2 in relation to $\left[{\rm{M/H}}\right]$.

We note that we have also inspected available measurements for mass, $\mathrm{log}g$, and luminosity of the stars belonging to G1 and G2, and found that while the median values for mass and $\mathrm{log}g$ of the two groups are very similar, pointing to M ≃ 1 M_⊙ and $\mathrm{log}g\simeq 4.44\,$ dex, there is a ∼14% difference in luminosity (median values L_G1 ≃ 0.87 L_⊙ and L_G2 ≃ 0.76 L_⊙). This not only suggests that it is the metallicity, and not the stellar mass, that could be playing a role in the different dynamo mechanism in action between the two groups, but also that there is a qualitative agreement among the trends observed in $\left[{\rm{M/H}}\right]$, age, and luminosity.

In this study we have shown that stars can be classified into two groups according to whether the period of the activity cycle decreases or increases with rotation. In both cases the simple scaling law represented in Equation (1) seems to hold for a wide range of rotations. Albeit the average rotation rate of the two groups is nearly the same; cycle periods for the group G1 are about a factor 2 shorter than for the group G2. We suggest that this dichotomy is likely originating from a different chemical composition of the stars. Finally, we would like to stress that the Sun does not occupy any special position in this context.

Can these two scaling laws be understood in terms of dynamo theory? As discussed in the introduction the positive slope can be explained already in terms of simple kinematic dynamo theory. The negative slope does not have this interpretation. However if we assume that the dynamo in the Sun is ruled by the meridional circulation at the bottom of the convection zone, then a higher rotation rate would imply weaker meridional circulation (by angular momentum conservation). This would produce a smaller ω_cyc because ω_cyc is determined by the strength of the meridional circulation (Dikpati et al. 2001; Bonanno 2013). What is the role of metallicity in this context? Stars with reduced metallicity have higher luminosity, therefore higher eddy diffusivity, at least according to mixing-length theory. On the contrary, a smaller luminosity implies a reduced diffusivity and therefore a stronger Reynolds number of the meridional circulation for this type of dynamo action.

We plan to further explore these arguments in detail in a future work that will also exploit a larger set of metallicity measurements as obtained by the soon to come ESA PLATO mission (Rauer et al. 2014).

We thank Gustavo Guerrero for comments on the manuscript, Arman Khalatyan, Hans-Erich Froehlich, Antonino Lanza, Sergio Messina, Santi Cassisi, and Isabella Pagano for helpful discussions. A special thanks to the anonymous referee who has helped us in greatly improving the manuscript.