https://orcid.org/0000-0002-2296-9141
ABSTRACT
Like the solar cycle, stellar activity cycles are also irregular. Observations reveal that rapidly rotating (young) Sun-like stars exhibit a high level of activity with no Maunder-like grand minima and rarely display smooth regular activity cycles. On the other hand, slowly rotating old stars like the Sun have low activity levels and smooth cycles with occasional grand minima. We, for the first time, try to model these observational trends using flux transport dynamo models. Following previous works, we build kinematic dynamo models of one solar mass star with different rotation rates. Differential rotation and meridional circulation are specified with a mean-field hydrodynamic model. We include stochastic fluctuations in the Babcock–Leighton source of the poloidal field to capture the inherent fluctuations in the stellar convection. Based on extensive simulations, we find that rapidly rotating stars produce highly irregular cycles with strong magnetic fields and rarely produce Maunder-like grand minima, whereas the slowly rotating stars (with a rotation period of 10 d and longer) produce smooth cycles of weaker strength, long-term modulation in the amplitude, and occasional extended grand minima. The average duration and the frequency of grand minima increase with decreasing rotation rate. These results can be understood as the tendency of less supercritical dynamo in slower rotating stars to be more prone to produce extended grand minima.
1 INTRODUCTION
Sun is not the only star that has a complex and variable magnetic field. Chromospheric emission of Ca ii H & K of many stars with spectral types from early F to M observed since 1966 through the HK Project of Mount Wilson Observatory (MWO) revealed the magnetic cycles (Baliunas et al. 1995). Other recent observations from the coronal X-ray emission (Wright et al. 2011; Wright & Drake 2016) and the magnetic field through the Zeeman-Doppler Imaging (ZDI) (Donati, Semel & Rees 1992; Vidotto et al. 2014) provide further evidence for magnetic activity in many stars studied earlier using the Ca ii H & K emission for the longer duration.
Stellar magnetic activity is observed to be largely controlled by rotation. The more rapidly a star rotates, the more active it is (Skumanich 1972; Rengarajan 1984). Noyes et al. (1984a) and Wright & Drake (2016) gave the activity-rotation relation using Ca ii H & K and X-ray emissions, respectively. They showed that the activity increases with the increase in rotation rate (or decrease in rotation period) for slow and moderate rotators, and then the activity tends to saturate for the fast-rotating stars.
While most of the earlier observations focused on the activity-rotation relation, some observations find a trend between the cycle duration and the rotation period of different Sun-like stars (Noyes, Weiss & Vaughan 1984b; Soon, Baliunas & Zhang 1994). Using the chromospheric activity from HARPS (High Accuracy Radial velocity Planet Searcher) and MWO data of 4454 cool stars, Boro Saikia et al. (2018) showed that the trend of the cycle period with the rotation period for the fast rotators is different from the slowly rotating stars. As the rotation period increases, the cycle period somewhat decreases for the rapidly rotating stars and increases for the slowly rotating stars. The trend is, however, quite complicated for fast rotators.
Our Sun shows the magnetic cycle of 22 yr period (11 yr in strength) with amplitude varying somewhat smoothly in the long-term (beyond the 11-yr period) and show occasional extended periods of weaker activity, the grand minima, e.g. Maunder minimum (Usoskin 2017). Different stars, in contrast, show a wide range of variability in the magnetic cycles. Baliunas et al. (1995) observed a smoother variability and occasional grand minima in the magnetic cycles of the slowly rotating (old) stars. On the other hand, they observed much irregular activity and no grand minima in the rapidly rotating young stars. Oláh et al. (2016), Boro Saikia et al. (2018), and Garg et al. (2019) also produced similar pieces of evidence using additional data. Recently, Baum et al. (2022) claimed that a K2V star, HD 166620, has entered into a grand minimum phase, and interestingly, it is a slowly rotating star. Also, Shah et al. (2018) suggests that HD 4915 is a possible Maunder minimum candidate, although its rotation period still needs to be confirmed.
A large-scale dynamo, powered by the helical convection and differential rotation, is responsible for the generation of the magnetic cycle in the Sun (Parker 1955). As the other Sun-like stars have convection zones (CZs) in their outer layers, it is natural to expect that these stars also support dynamo action through which the stellar magnetic cycles are maintained. Some of the stellar cycles (e.g. HD 10476, HD 16160, and so on) are so similar to the solar cycle (in terms of regular cyclic variation and obeying the Waldmeier effect, which says that strong cycles rise faster than the weaker ones; Garg et al. 2019) that it suggests a similar dynamo operating in Sun and other Sun-like stars (also see Jeffers et al. 2022, for a theoretical argument behind this expectation). The motivation of our work is to extract the dependence on the rotation rate of the Sun-like stars for their cycle variability and the occurrence of the grand minima using dynamo models.
Recently, the Babcock–Leighton mechanism (Babcock 1961; Leighton 1969), in which the tilted bipolar magnetic regions produce poloidal field in the Sun, has received strong observational supports (Dasi-Espuig et al. 2010; Kitchatinov & Olemskoy 2011; Priyal et al. 2014; Cameron & Schüssler 2015). Including this process for the generation of the poloidal field, the Babcock–Leighton type dynamo models have produced great successes in providing many observational features of the solar magnetic cycle, including the grand minima (e.g. Choudhuri & Karak 2012; Olemskoy & Kitchatinov 2013; Passos et al. 2014; Karak & Miesch 2017; Cameron & Schüssler 2017; Lemerle & Charbonneau 2017; Inceoglu, Arlt & Rempel 2017; Biswas, Karak & Cameron 2022).
In the past, Babcock–Leighton dynamo model has also been used to study the stellar magnetic cycles. For instance, Nandy & Martens (2007) employed a time-delay dynamo model to investigate the relationship between the magnetic field and cycle period with the dynamo number. Jouve, Brown & Brun (2010) utilized a kinematic Babcock–Leighton dynamo model to observe that the cycle period increases as the rotation rate increases, unless the meridional flow speed is assumed to increase with the rotation rate, which contradicts theoretical results (Miesch 2005; Brown et al. 2008). Later, by employing the non-local and distributed α effects in non-linear α2 dynamo models for moderate to slowly rotating stars, Pipin (2015) find some agreement of the cycle period versus rotation period with observation. Karak, Kitchatinov & Choudhuri (2014a) constructed the Babcock–Leighton type flux transport dynamo model for Sun-like stars with different rotation periods by including differential rotation and meridional circulation from corresponding hydrodynamical models of Kitchatinov & Olemskoy (2011). They managed to reproduce the activity–rotation relation correctly but again failed to reproduce the cycle period versus rotation rate relation. Recently, Hazra et al. (2019) performed simulations using the same model but by including a radial pumping near the surface of the stars. They found an increasing trend of the cycle period with an increase in the rotation period for the slowly rotating stars, and a decreasing trend in the cycle period for the rapidly rotating stars; also see Do Cao & Brun (2011) who included latitudinal pumping and found some agreement with observations. Kitchatinov (2022) studied the stellar activity cycles using a Babcock–Leighton type dynamo model and showed a strong temperature dependence on the cycle period. Karak, Tomar & Vashishth (2020) and Noraz et al. (2022) applied mean-field models in different stars and addressed the effect of antisolar differential rotation (which naturally arises in the high-Rossby number convection, e.g. Gastine et al. 2014; Brun et al. 2017; Karak, Miesch & Bekki 2018) on the polarity reversal and the strength of magnetic field in slowly rotating stars. However, to the best of our knowledge, no previous study was performed to explore the cycle variability and grand minima in stars.
In recent years, global MHD convection simulations have produced some exciting results of the stellar magnetic cycles (Karak et al. 2015; Augustson et al. 2015; Käpylä et al. 2016; Strugarek et al. 2018; Viviani et al. 2019; Brun et al. 2022). Warnecke et al. (2018) analysed how the magnetic cycle period changes as a function of the Rossby number. Viviani et al. (2018) studied the simulations of different stars and showed the transition of the magnetic field from axi- to non-axisymmetric field configuration at around 1.8 times the solar rotation rate, where the differential rotation changes from solar to antisolar. However, these MHD simulations are not capable of reproducing some basic features of the solar cycle (e.g. the 11-yr periodicity with regular reversal, equatorward migration of toroidal field at low latitude, poleward migration of surface radial field, largely dipolar field) and the correct flow (particularly the observed amplitude of the convective flow) robustly, there lies the uncertainty of whether these results hold in stars. Being computationally expensive, these simulations were not applied to study the long-term variability of stellar cycles; however, see Passos & Charbonneau (2014), Augustson et al. (2015), and Käpylä et al. (2016), who have performed simulations for several cycles, which may be used for long-term studies.
In this paper, we apply the models of Karak et al. (2014a) and Hazra et al. (2019) to study the irregularities of the stellar cycle, in particular, how the variability and the frequency of grand minima change with the stellar rotation. For this, we shall include the stochastic fluctuations to capture the inherent randomness in the stellar convection (Choudhuri 1992) as seen in the form of noise in the flux emergence and the tilts of BMRs around Joy’s law (Dasi-Espuig et al. 2010; Stenflo & Kosovichev 2012; McClintock, Norton & Li 2014; Wang et al. 2015; Arlt et al. 2016; Jha et al. 2020) (Section 2). We shall see that our models produce a strong magnetic activity and highly irregular cycles in rapidly rotating stars and, on the contrary, a weak magnetic activity and more regular cycles in slowly rotating stars (Section 3.1). Maunder-like extended grand minima are only produced in slowly rotating stars (with rotation period of 10 d and longer), and the frequency of occurrence of these events increases with the increase in the rotation period of the star (Section 3.3).
2 MODEL
We build our model based on Karak et al. (2014a) where the following equations for the axisymmetric magnetic field are evolved,
where Bp = ∇ × [A(r, θ)eϕ] is the poloidal component of the magnetic field and B(r, θ) is the toroidal component, s = rsin θ, |${\bf v_p} = v_r { \hat{r}} + v_{\theta } {\bf \hat{\theta }}$| is the meridional circulation, and Ω is the angular velocity. While Ω is well-measured in the whole CZ of the Sun, the meridional flow is only constrained in the near-surface layer of the Sun. Recent helioseismic studies for the deep meridional circulation indicate a single-cell flow in the solar CZ (Rajaguru & Antia 2015; Gizon et al. 2020). Observations for other stars, on the other hand, are limited to the surface differential rotation only. Global MHD simulations for the Sun-like stars provide differential rotation, which often shows transition from solar to antisolar profile with increasing Rossby number near the solar value, and the meridional flow, which is multicellular and time varying (e.g. Featherstone & Miesch 2015; Karak et al. 2015; Viviani et al. 2019). Therefore, for vp and Ω, we use the data from a mean-field hydrodynamic model of Kitchatinov & Olemskoy (2011). This numerical model jointly solves the mean-field equations for the angular velocity, meridional flow, and heat transport in a spherical layer of a stellar CZ. The model produces the solar-type differential rotation as a consequence of angular momentum fluxes. The one-cell per hemisphere meridional flow predicted by the model for the Sun agrees with the recent seismological detection (Rajaguru & Antia 2015; Gizon et al. 2020). Also, the computed dependence of differential rotation on stellar rotation rate and spectral type (Kitchatinov & Olemskoy 2012b) is in at least qualitative agreement with observations by Barnes et al. (2005) and Balona & Abedigamba (2016). The model, however, does not produce the tachocline self-consistently; rather, it has the lower boundary at 0.72Rs (Rs being the stellar radius). Therefore, the tachocline in our dynamo model is formed by smoothly varying the Ωmodel(r, θ) from the differential rotation model at r = 0.72Rs to the value of the rotation rate at the core Ωcore in the following way:
The η is the turbulent magnetic diffusivity and is taken as a function of r alone, having the following form:
with rBCZ = 0.7Rs, dt = 0.015Rs, d2 = 0.05Rs, and rsurf = 0.95Rs; ηRZ, ηSCZ, and ηsurf represent the diffusivities, at the inner boundary, within CZ, and at the surface, respectively, having the values as ηRZ = 5 × 108 cm2 s−1, ηSCZ = 5 × 1010 cm2 s−1, and ηsurf = 2 × 1012 cm2 s−1. The diffusivity profile of equation (4) approaches ηRZ at the inner boundary of 0.6Rs, remains at ηSCZ in the bulk of CZ and increases to ηsurf at the surface (see Fig. 5 in Karak et al. 2014a). We have ignored the change of diffusivity with the rotation rate and the magnetic field in our study due to its limited knowledge in estimating its value in different stars.
The term S(r, θ; B) is the source for the generation of poloidal field, which captures the Babcock–Leighton mechanism in our axisymmetric model, and it is given by,
where α0 is the strength of Babcock–Leighton process, |$\overline{B}(r_t,\theta)$| is the toroidal field at latitude θ averaged over the whole tachocline r = 0.685Rs to r = 0.715Rs. From equation (5), we observe that when the magnetic field becomes comparable to B0 (the saturation field strength), the nonlinearity becomes important and the field eventually tends to hover around B0. Therefore, everywhere in our study, we measure the magnetic field in the unit of B0. While in the traditional α effect based on the helical convection, above α quenching is obvious, the Babcock–Leighton α also experiences a quenching due to the fact that BMRs with strong fields rise quickly and Coriolis force gets less time to induce tilt and the strong magnetic field also gives more tension, which causes less tilt (D’Silva & Choudhuri 1993; Fan, Fisher & McClymont 1994). Observations indeed find some evidence of tilt quenching (Jha et al. 2020). Furthermore, the latitudinal variation of BMRs (stronger cycles produce BMRs at higher latitudes) gives rise to quenching in the Babcock–Leighton process (Jiang 2020; Karak 2020).
We have limited knowledge about the Babcock–Leighton process in other stars, and thus we are not sure of how the strength of this process changes with the rotation of the star. On theoretical grounds, we expect the tilt of bipolar magnetic regions to increase with the increase of the rotation rate of stars (D’Silva & Choudhuri 1993; Kitchatinov & Olemskoy 2015). However, there is an opposing effect that arises due to the fact that with the increase of rotation rate, the latitudes of BMR emergences are expected to shift to higher latitudes (Schuessler & Solanki 1992) and higher latitudes BMRs are less efficient in producing polar field (Jiang, Cameron & Schüssler 2014; Karak 2020). In our study, the strength of Babcock–Leighton process α0, is chosen to depend on the rotation in the following way,
where α0, s is the value of α0 for the solar case and P⊙ and Prot are the rotation period of Sun and the star, respectively.
The Babcock–Leighton process includes considerable randomness, primarily due to irregular variations in the tilts and emergence rates of the bipolar active regions (Jiang et al. 2014). Therefore, we include fluctuations in the α appearing in equation (6) in the following way.
where Xn is the Gaussian random number with a mean of unity and standard deviation (σ) of 2.67. This value of σ is inspired by the study of Olemskoy, Choudhuri & Kitchatinov (2013), who computed the fluctuations in the Babcock–Leighton process by estimating the contribution of the sunspot group to the polar field using the data of Royal Greenwich, Kodaikanal, and Mount Wilson Observatories. We keep the value of σ same throughout all the simulations presented in this paper. In our models, the value of α0 is updated randomly after a certain time step which we take to be 1 month.
It may be noted that other model parameters like eddy diffusivity or turbulent pumping can also include fluctuations due to randomness inherent to turbulent stellar convection. However, the fluctuations are relatively small and less consequential compared to the fluctuations in the Babcock–Leighton mechanism. In particular, fluctuations in the angular momentum fluxes produce only small variations in the differential rotation and moderate fluctuations in the meridional flow (Rempel 2005). Thus, these fluctuations are neglected. The differential rotation and meridional flow are steady and equator-symmetric in our model.
2.1 Model I
In this case, we use the same model as given in Karak et al. (2014a) except the value of the strength of the Babcock–Leighton α, and we add fluctuations in the Babcock–Leighton process. In this model, the value of α0 is taken as 0.9 cm s−1 [instead of 1.6 cm s−1, which was used in Karak et al. (2014a)].
As in Babcock–Leighton models, α captures the average effect of the decay of tilted BMRs, it must be non-zero only near the surface, and it must have cos θ dependence due to the angular dependence of the Coriolis force, which is the possible cause of BMR tilt. However, to suppress the poloidal field generation in high latitudes (as BMRs do not appear in high latitudes), a sin θ factor is also introduced (see, e.g. Dikpati & Charbonneau 1999). Therefore in this model, the profile of α is given by,
with r4 = 0.95Rs, r5 = Rs, d4 = 0.05Rs, and d5 = 0.01Rs.
2.2 Model II
In the previous model (Model I), one hemisphere of the Sun was used to study the dynamo, for which a dipolar boundary condition was imposed at the equator. This did not allow us to observe the magnetic field configuration across the equator and the hemispheric asymmetry of the magnetic field, which are very relevant for the solar and stellar observations (DeRosa, Brun & Hoeksema 2012). Therefore, in Model II, we extend the same Model I to the full sphere of the Sun, and we include the fluctuations separately in the two hemispheres. Other than extending Model I to the full sphere, and thus eliminating the equatorial boundary condition, no other changes are made.
2.3 Model III
Finally, we take Model III, which is same as Model II but at increased diffusivity and added radial magnetic pumping. This inclusion of pumping is inspired by Hazra et al. (2019), who found some agreement of the cycle period versus rotation trend with observations. It was realized that a downward magnetic pumping helps to make the magnetic field radial near the surface and reduce the toroidal flux loss through the surface, making the dynamo model in accordance with the surface flux transport models and observations (Cameron et al. 2012). The near-surface pumping also helps the dynamo to operate at a high diffusivity range consistent with the mixing-length theory (Kitchatinov & Olemskoy 2012a; Karak & Cameron 2016; Karak & Miesch 2017), and facilitates the model to recover from the Maunder-like extended grand minima (Karak & Miesch 2018). The pumping has the following form:
where the amplitude of the radial magnetic pumping is given by γ0, which is 24 m s−1 in all the stars.
We note that we do not use the exact same model of Hazra et al. (2019) because, in that model, when we include fluctuations, even the solar case does not produce the dipolar field as seen in the observations. Therefore, to obtain the dipolar field, we reduce the diffusivity for the bulk of the CZ by taking the following parameters (equation 4): ηSCZ = 3 × 1011 cm2 s−1, and ηsurf = 3 × 1012 cm2 s−1. We note that with these parameters, the diffusivity in the whole CZ is about six times stronger than that used in Models I and II.
The α profile used in this model is given by,
where d = 0.01Rs. The α0 has the same form (equation 6) as in Model I, except in this case, α0, s = 4 cm s−1 and fluctuations in this model are included separately in the two hemispheres. We note that above α in equation (10) has a sin 2θcos θ dependence instead of sin θcos θ as used in Models I and II to make the α effect strong (weaker) in low (high) latitudes. Also, the radial extent of this α is a bit wider than that used in Models I and II.
After specifying all the parameters, we solve the above equations (1) and (2) numerically in the SCZ with the radial extent of 0.55Rs to Rs and the following boundary conditions. We take, at the lower boundary: A = B = 0, at the top (surface) layer: B = Bθ = 0 (i.e. radial boundary condition), at poles: B = A = 0 (i.e. no singularity) and at the equator for Model I: B = 0 = ∂A/∂θ = 0 (dipolar condition). Simulations are performed in 129 × 129 grid points in radial and latitudinal directions.
3 RESULTS AND DISCUSSIONS
We initialized our dynamo simulations for stars of solar mass with rotation periods of 1 d, 3 d, 7 d, 10 d, 15 d, 20 d, 25.38 d (Sun), and 30 d, respectively, by computing differential rotation and meridional circulation from the mean-field hydrodynamic model of Kitchatinov & Olemskoy (2011). The following sections discuss the various aspects of magnetic cycles obtained from all three models.
3.1 Magnetic field morphology
Fig. 1 depicts the butterfly diagrams of the toroidal field at the base of CZ for the rotation periods 1, 7, 25.38, and 30 d from Model I. After analysing these panels, we find the regular polarity reversal. However, in Model II (see Fig. 2) the magnetic cycles are a bit irregular. We observe a strong hemispheric asymmetry in the magnetic field. Sometimes the magnetic field in one hemisphere is largely suppressed or enhanced. Hence, in Model II, the polarity reversal is not regular. Finally, for Model III (Fig. 3), we observe that for the rapidly rotating case and the Sun, regular reversal is seen, but for intermediate and slow rotators, the polarity reversal is not regular. The magnetic field distribution for the rotation period of 7 and 30 d show extended cycles. However, for the 7 d case, the magnetic field is largely quadrupolar, while for the 30 d case, it is largely bipolar. We observe that the magnetic field distribution at higher radial layers is largely different. We also note that for the slowly rotating stars with a rotation period of ≥ 10 d, the magnetic field is largely dipolar. In contrast, for the rotation period of 7 d and less, the parity is changed to largely quadrupolar, although there are sometimes when the parity remains dipolar.
Time–latitude plots of toroidal field at r = 0.71Rs (in the unit of B0) for different stars with rotation period of (a) 1 d, (b) 7 d, (c) the solar value (i.e. 25.38 d), and (d) 30 d for Model I.
Same as Fig. 1, but for Model II.
Same as Fig. 1, but for Model III.
In all the models for the slowly rotating stars, we observe an equatorward migration of the toroidal field at the low latitudes as a consequence of the transport by the equatorward meridional circulation (Fig. 3 of Karak et al. 2014a). However, for the rapidly rotating stars, we find a weak poleward migration of the field in high latitudes, particularly see Fig. 3(a). This poleward migration is due to the diffusion of the field from the mid-latitude where the toroidal field generation is strongest (also see Figs 6 and 10 in Hazra et al. 2019). Moreover, the meridional flow is much weaker in rapidly rotating stars.
One obvious feature in these simulations is that the magnetic field becomes strong in fast-rotating stars. This happens because the strength of α increases with the rotation rate of the star. This increase in the magnetic field is congruous with the observations (Noyes et al. 1984a; Wright et al. 2011). While in observations, the magnetic field is saturated in rapidly rotating stars, our model always produces an increasing trend with the rotation rate. This is because, in our model, the latitude of operation of the Babcock–Leighton process (the band of BMR emergences) is fixed in all stars, while in observations, it increases. Thus, the generation of poloidal field is less efficient in rapidly rotating stars (Kitchatinov & Olemskoy 2015).
These time–latitude plots also give a hint about the variability observed in different stars. Slowly rotating stars seem to produce more long-term modulation in their cycles, including extended episodes of weaker magnetic field. In contrast, fast rotators generate less modulation. This result is in agreement with observations. The root cause for such behavior is that the slowly rotating stars have a small dynamo number. Due to this, if the magnetic cycle gets weaker sometimes, then it would take a long time to grow the field. Therefore, we see a long-term modulation in the slowly rotating stars. But for the fast rotators, the cycle recovers its strength quickly after getting into the weak phase. In fast rotators, the dynamo number is high so the growth rate is very high. This trend is also explained in Kumar, Karak & Vashishth (2021) and Vashishth, Karak & Kitchatinov (2021). The results are in accordance with the observations as well (Baliunas et al. 1995). A detailed discussion of the cycle variability is made in Section 3.3.
3.2 Cycle duration versus rotation period
We now compute the cycle periods for all three models. This is done by determining the peak of the Fourier power spectrum of the time series of the toroidal field over the tachocline for both the Northern and Southern hemispheres separately. However, due to the irregular nature of the cycles in Model III, we fail to identify a prominent peak in the power spectrum, and hence cannot identify the dominant cycle period for the stars having a rotation period of less than 10 d. Fortunately, we are able to find a range in which the cycle periods of these stars could lie. The computed cycle periods for all the models are listed in Table 1, and the variations with the rotation rate are shown in Fig. 4. This figure infers an increasing trend of the cycle period with the stellar rotation rate for both Models I and II. This increasing trend is quicker in the fast-rotating stars and milder in the slowly rotating stars. This happens because, as the rotation period decreases (or rotation rate increases), the meridional flow becomes weaker (although the flow speed increases in the thin layers near the top and bottom boundaries).
Variations of the activity cycle period (Pcyc in years) with rotation period (Prot in days) for (a) Models I (filled circles) and II (asterisks) and (b) Model III (solid and dashed lines are for Northern and Southern hemispheres, respectively).
Summary of simulations. Here, Prot is the rotation period of the star in days, and Pcyc is the mean magnetic cycle period in each hemisphere. For each model, the number of grand minima and parity are computed from the surface radial and the toroidal fields at the base of CZ, which are separated by a comma.
| Prot (d) | Pcyc (yr) | Number of grand minima | Parity | |||||
|---|---|---|---|---|---|---|---|---|
| Model I | Model II (N, S) | Model III (N, S) | Model I | Model II | Model III | Model II | Model III | |
| 1 | 13.68 | 12.89, 12.76 | 7.81–9.18, 7.26–9.45 | 0, 0 | 0, 0 | 0, 0 | 0.013, −0.053 | 0.434, 0.410 |
| 3 | 11.27 | 10.64, 12.25 | 7.87–8.78, 7.63–8.91 | 0, 0 | 0, 0 | 0, 0 | 0.003, −0.043 | −0.726, −0.742 |
| 7 | 9.15 | 8.46, 9.15 | 7.69–8.56, 8.26–9.25 | 0, 0 | 0, 0 | 0, 0 | −0.015, −0.079 | −0.685, −0.676 |
| 10 | 7.62 | 7.96, 7.42 | 10.09, 11.16 | 2, 4 | 0, 0 | 2, 2 | 0.009, −0.032 | −0.786, −0.788 |
| 15 | 6.92 | 7.13, 6.82 | 10.47, 10.71 | 6, 15 | 2, 4 | 3, 4 | 0.001, −0.038 | −0.844, −0.846 |
| 20 | 6.44 | 6.15, 6.45 | 12.69, 12.59 | 17, 22 | 6, 12 | 6, 8 | 0.009, −0.011 | −0.831, −0.830 |
| 25.38 (Sun) | 5.65 | 6.15, 5.79 | 11.78, 12.55 | 32, 40 | 16, 21 | 8, 12 | −0.117, −0.194 | −0.699, −0.830 |
| 30 | 5.65 | 5.74, 5.65 | 12.24, 13.80 | 36, 41 | 19, 24 | 10,16 | −0.138, −0.201 | −0.889, −0.888 |
| Prot (d) | Pcyc (yr) | Number of grand minima | Parity | |||||
|---|---|---|---|---|---|---|---|---|
| Model I | Model II (N, S) | Model III (N, S) | Model I | Model II | Model III | Model II | Model III | |
| 1 | 13.68 | 12.89, 12.76 | 7.81–9.18, 7.26–9.45 | 0, 0 | 0, 0 | 0, 0 | 0.013, −0.053 | 0.434, 0.410 |
| 3 | 11.27 | 10.64, 12.25 | 7.87–8.78, 7.63–8.91 | 0, 0 | 0, 0 | 0, 0 | 0.003, −0.043 | −0.726, −0.742 |
| 7 | 9.15 | 8.46, 9.15 | 7.69–8.56, 8.26–9.25 | 0, 0 | 0, 0 | 0, 0 | −0.015, −0.079 | −0.685, −0.676 |
| 10 | 7.62 | 7.96, 7.42 | 10.09, 11.16 | 2, 4 | 0, 0 | 2, 2 | 0.009, −0.032 | −0.786, −0.788 |
| 15 | 6.92 | 7.13, 6.82 | 10.47, 10.71 | 6, 15 | 2, 4 | 3, 4 | 0.001, −0.038 | −0.844, −0.846 |
| 20 | 6.44 | 6.15, 6.45 | 12.69, 12.59 | 17, 22 | 6, 12 | 6, 8 | 0.009, −0.011 | −0.831, −0.830 |
| 25.38 (Sun) | 5.65 | 6.15, 5.79 | 11.78, 12.55 | 32, 40 | 16, 21 | 8, 12 | −0.117, −0.194 | −0.699, −0.830 |
| 30 | 5.65 | 5.74, 5.65 | 12.24, 13.80 | 36, 41 | 19, 24 | 10,16 | −0.138, −0.201 | −0.889, −0.888 |
Summary of simulations. Here, Prot is the rotation period of the star in days, and Pcyc is the mean magnetic cycle period in each hemisphere. For each model, the number of grand minima and parity are computed from the surface radial and the toroidal fields at the base of CZ, which are separated by a comma.
| Prot (d) | Pcyc (yr) | Number of grand minima | Parity | |||||
|---|---|---|---|---|---|---|---|---|
| Model I | Model II (N, S) | Model III (N, S) | Model I | Model II | Model III | Model II | Model III | |
| 1 | 13.68 | 12.89, 12.76 | 7.81–9.18, 7.26–9.45 | 0, 0 | 0, 0 | 0, 0 | 0.013, −0.053 | 0.434, 0.410 |
| 3 | 11.27 | 10.64, 12.25 | 7.87–8.78, 7.63–8.91 | 0, 0 | 0, 0 | 0, 0 | 0.003, −0.043 | −0.726, −0.742 |
| 7 | 9.15 | 8.46, 9.15 | 7.69–8.56, 8.26–9.25 | 0, 0 | 0, 0 | 0, 0 | −0.015, −0.079 | −0.685, −0.676 |
| 10 | 7.62 | 7.96, 7.42 | 10.09, 11.16 | 2, 4 | 0, 0 | 2, 2 | 0.009, −0.032 | −0.786, −0.788 |
| 15 | 6.92 | 7.13, 6.82 | 10.47, 10.71 | 6, 15 | 2, 4 | 3, 4 | 0.001, −0.038 | −0.844, −0.846 |
| 20 | 6.44 | 6.15, 6.45 | 12.69, 12.59 | 17, 22 | 6, 12 | 6, 8 | 0.009, −0.011 | −0.831, −0.830 |
| 25.38 (Sun) | 5.65 | 6.15, 5.79 | 11.78, 12.55 | 32, 40 | 16, 21 | 8, 12 | −0.117, −0.194 | −0.699, −0.830 |
| 30 | 5.65 | 5.74, 5.65 | 12.24, 13.80 | 36, 41 | 19, 24 | 10,16 | −0.138, −0.201 | −0.889, −0.888 |
| Prot (d) | Pcyc (yr) | Number of grand minima | Parity | |||||
|---|---|---|---|---|---|---|---|---|
| Model I | Model II (N, S) | Model III (N, S) | Model I | Model II | Model III | Model II | Model III | |
| 1 | 13.68 | 12.89, 12.76 | 7.81–9.18, 7.26–9.45 | 0, 0 | 0, 0 | 0, 0 | 0.013, −0.053 | 0.434, 0.410 |
| 3 | 11.27 | 10.64, 12.25 | 7.87–8.78, 7.63–8.91 | 0, 0 | 0, 0 | 0, 0 | 0.003, −0.043 | −0.726, −0.742 |
| 7 | 9.15 | 8.46, 9.15 | 7.69–8.56, 8.26–9.25 | 0, 0 | 0, 0 | 0, 0 | −0.015, −0.079 | −0.685, −0.676 |
| 10 | 7.62 | 7.96, 7.42 | 10.09, 11.16 | 2, 4 | 0, 0 | 2, 2 | 0.009, −0.032 | −0.786, −0.788 |
| 15 | 6.92 | 7.13, 6.82 | 10.47, 10.71 | 6, 15 | 2, 4 | 3, 4 | 0.001, −0.038 | −0.844, −0.846 |
| 20 | 6.44 | 6.15, 6.45 | 12.69, 12.59 | 17, 22 | 6, 12 | 6, 8 | 0.009, −0.011 | −0.831, −0.830 |
| 25.38 (Sun) | 5.65 | 6.15, 5.79 | 11.78, 12.55 | 32, 40 | 16, 21 | 8, 12 | −0.117, −0.194 | −0.699, −0.830 |
| 30 | 5.65 | 5.74, 5.65 | 12.24, 13.80 | 36, 41 | 19, 24 | 10,16 | −0.138, −0.201 | −0.889, −0.888 |
Although these two models reproduce various stellar observations, they fail to reproduce the magnetic cycle period versus rotation trend correctly for the slowly rotating stars. Limited observations (Boro Saikia et al. 2018) seem to show a rapid increase in the cycle period with the increase of the rotation rate for fast-rotating stars. This is consistent with the trend found in our Models I and II. However, the observed data for slow rotators show an increasing trend of the activity cycle period with the increase in rotation period, which is opposite to the findings in Models I and II. One way to resolve this discrepancy is to include radial magnetic pumping in the stellar CZs. Hazra et al. (2019), after including the pumping, got a trend somewhat similar to the observations. In our Model III, after including radial pumping, we also got the cycle-rotation period trend closer to the observations. It is also possible that the observed trend is a consequence of strong decrease of the cycle period with stellar effective temperature and on average faster rotation of hotter stars (Kitchatinov 2022). A different trend of the cycle period with the rotation period found in Model III is due to the operation of the dynamo in a pumping-dominated regime. When strong downward magnetic pumping is included in this model, the diffusion of the magnetic field across the surface becomes negligible and then the dynamo allows it to operate at a low α (Karak & Cameron 2016). Lower the α longer is the cycle period. We can see from Fig. 4 that at 30 d rotation period, while Models I and II were producing a cycle period of 6 yr, Model III produced a much longer period of 13 yr. Then with the decrease of the rotation period, the α becomes stronger, and thus the poloidal field generation process becomes more efficient. This makes the reversal of the field faster. This effect in the pumping-dominated regime overpowers the increase of the cycle period due to a decrease in meridional flow speed.
3.3 Variability and grand minima occurrence
We now come across the central question of our study, i.e. how the long-term variability of the stellar cycles changes with the rotation rates of the stars. In our model, the cycle variability is produced due to the randomness in the Babcock–Leighton mechanism. To analyse the long-term variability qualitatively, we carefully observe the time series data of the toroidal magnetic flux of the Northern hemisphere and the absolute radial magnetic field averaged over the whole surface from a simulation of 11 000 yr. Figs 5 and 6 show the discussed time series plots of the toroidal flux for Model I and the radial field from Model III.
Time series plot along with its smoothed variation of toroidal magnetic field for Model I of stars having rotation period of (a) 1 d, (b) 7 d, (c) 25.38 d (the solar value), and (d) 30 d. The dark-red bars highlight the extended weaker activity episodes i.e. grand minima in each case.
Same as Fig. 5 but computed from the absolute radial magnetic field, averaged over the whole surface for Model III.
From these figures, we can easily see that the fast rotators produce irregular cycles with smaller long-term variability, and on the other hand, slowly rotating stars produce more long-term modulation in their cycles with episodes of extended weak magnetic fields.
To make a quantitative estimate of the irregularities and the long-term memory in the time series for each case, we compute the well-known Hurst exponent (H), which gives a measure of the temporal memory or persistence in the time series. A value of H = 0.5 implies that the time series is obtained from a memoryless random process. On the other hand, if a time series gives H > 0.5, then it suggests to have persistence. When a system has a memory that depends on the previous step, it is said to be persistent, and thus in the time series if there was an increase in the value, it is more likely that the following step(s) will increase as well. In this case, the time series will cover more ‘distance’ than a random walk can. When H < 0.5, the opposite applies and it is said to be antipersistent. When there is a long-term modulation in the stellar cycle data, we expect a memory and the value of H should be larger than 0.5. A larger value of H implies large long-term memory in the stellar cycle.
To obtain the Hurst exponent, we use the famous R/S method as given in Mandelbrot & Wallis (1969) and applied to many solar data in the past (Ruzmaikin, Feynman & Robinson 1994; Suyal, Prasad & Singh 2009; Das, Ghosh & Karak 2022). To do so, we first bin the data by using a bin-size of half a year. Then, to evaluate H, the binned time series is divided into several shorter time series of length τ = 50. The average re-scaled range (R/S) is then calculated for each temporal window τ. At last, the slope of the log(R/S) versus log(τ) values gives the value of H.
Fig. 7 shows the values of H evaluated in all the models for each star. A linear fit to all the data shows the overall increase in H with the rotation period. After analysing this figure, we come to the conclusion that for the rapidly rotating stars, there is little long-term modulation and cycles are more irregular. Whereas a long-term memory is seen in slowly rotating stars. Hence, the persistence increases as the rotation period increases. Interestingly, these results are in-tune with the observations (Baliunas et al. 1995; Oláh et al. 2016).
Variation of Hurst exponent with respect to the rotation rate along with the linear-fit curve of all the models.
Finally, we identify the grand minima from the time series of the toroidal field at the base of CZ and the surface radial field. For this, we employ the same method as used in Usoskin, Solanki & Kovaltsov (2007) for the Sun; i.e. we first bin the data by using a bin-size of the duration of one cycle, then we filter the data by using Gleissberg’s low-pass filter 1-2-2-2-1. This gives us the smoothed data. Finally, the portion of this data that falls below 50 per cent of its mean for at least two consecutive cycle periods, is considered as a grand minimum. Later, we count these numbers of grand minima to evaluate the frequency of occurrence of grand minima in each case. The computed number of grand minima for all the models are listed in Table 1.
From Fig. 8(a), we infer that in all the models, the number of grand minima increases with the increase in rotation period. Rapidly rotating stars hardly produce any grand minima, in fact, stars with a rotation period of 7 d or less, do not produce any Maunder-like grand minima. On the other hand, slowly rotating stars produce some grand minima with an increasing trend with the rotation period. This is because, with the increase of rotation period, the supercriticality of the dynamo decreases, and the dynamo is more prone to produce extended grand minima in this regime. This result is as per Vashishth et al. (2021), where we observe a decrease in the frequency of occurrence of grand minima as the supercriticality increases. To check the robustness of this result, we run our simulations with another set of Gaussian random numbers having the same mean and σ as the previous case. We again find the same conclusion that stars with Pcyc ≤ 7 d, do not produce grand minima, and the number of grand minima increases with the increase in rotation period.
Change of (a) the number and (b) the average duration of grand minima with the rotation period of stars. Yellow circles, blue asterisks, and red diamonds depict the trends for Models I, II, and III, respectively. In (b), the error bars are computed from the standard deviation of the durations of the grand minima in each case.
Additionally, we estimated the change in the average duration of grand minima with the stellar rotation. This variation is depicted in Fig. 8(b). Similar to the frequency, the average duration of grand minima increases with the increase of rotation period for all models except for the rotation period of 30 d in Model III. The association of the frequency of occurrence of the grand minima with the duration of these events is shown via histogram in Fig. 9. With the help of Fig. 8(b) and Fig. 9, we can easily infer that the duration of the stellar grand minima falls mostly below 150 yr, while the average value lies below 70 yr only. In Model I, we see a few grand minima occurring for longer duration due to the absence of hemispheric coupling. In full sphere models (Models II and III), hemispheric coupling helps to recover the model from extended grand minima easily (Karak & Miesch 2018), and thus this does not produce very long grand minima. Further, in the Sun, we get about 10–40 grand minima (depending on which model we are considering), while in observations, this number is 27. And the result, that most of the solar grand minima hover below 150 yr is in agreement with observations (Usoskin 2017).
Relation between frequency of occurrence of grand minima with the corresponding duration in Models I–III for the stars in which grand minima are observed.
4 CONCLUSIONS AND DISCUSSION
From our extensive simulations of the kinematic flux transport dynamo model with stochastically forced Babcock–Leighton source for the stars of 1M⊙ mass with rotation periods of 1, 3, 7, 10, 15, 20, 25.38 (solar value), and 30 d, we make the following inferences.
Rapidly rotating stars produce a strong magnetic field and the strength of the field increases with the increase of rotation rate, which is in accordance with the observations (Noyes et al. 1984a). The increase of field with the increase of rotation rate in our model is due to the enhancement of the strength of the Babcock–Leighton source.
The cycle period increases with the increase of the rotation rate of the stars (in Models I and II) due to the weakening of the meridional circulation. However, when the downward turbulent magnetic pumping is included (in Model III), cycles become longer in slowly rotating stars and shorter (although become very irregular) in rapid rotators. Thus, pumping helps to bring the results closer to observations (Boro Saikia et al. 2018) as also suggested by Hazra et al. (2019).
Strong hemispheric asymmetry is produced in the magnetic field for all the stars. In general, the quadrupolar field dominates in the rapidly rotating stars and the dipolar field dominates in the Sun and slowly rotating stars.
In rapidly rotating stars, the stellar magnetic cycles are highly irregular, while in slowly rotating stars cycles are more regular and the cycle amplitude displays a smooth long-term modulation. These results are consistent with the stellar observations (Baliunas et al. 1995; Boro Saikia et al. 2018; Oláh et al. 2016; Garg et al. 2019).
Only slowly rotating stars with rotation period ≥10 d produce grand minima. The number and the average duration of grand minima increase with the increase of the rotation period of the stars. This is again supported by the available observations because the confirmed Maunder minimum candidates are only slow rotators [Sun (Prot = 25.38 d), HD 166620 (Prot = 45 d); Baum et al.2022].
The length of the stellar grand minima lies mostly below 150 yr. However, in the one hemisphere (Model I) model, several grand minima occur with longer duration. The average duration of grand minima in this model is longer than the other two models because the hemispheric coupling is absent in this model. In full sphere models (Models II and III), hemispheric coupling helps to recover the model from extended grand minima easily (Karak & Miesch 2018; Hazra & Nandy 2019), and thus this does not produce very long grand minima. The result that most of the solar grand minima hover below 150 yr is in agreement with the reconstructed solar activity data (Usoskin2017).
Although many results of stellar cycles are robust and congruous with observations, there are limitations to our study. First, we have considered the only nonlinearity through the standard α quenching and we have ignored the nonlinear feedback of the magnetic field on large-scale flows. While, in the Sun, this is not a concern, in the rapidly rotating stars having a strong magnetic field, this nonlinearity can have a serious impact in producing cycle irregularity. Second, the turbulent transport coefficients are expected to change with the rotation and the magnetic field (Kitchatinov, Pipin & Ruediger 1994; Karak et al. 2014b) and they can change the dynamo properties. Due to limited knowledge of their variations in different stars, we have not changed their values in our models. Third, the level of stochastic noise is kept constant in all the stars, again due to its limited knowledge. Fourth, we have not considered the turbulent α effect, which in the Sun is negligible in comparison to the Babcock–Leighton process (Cameron & Schüssler 2015), but may becomes increasingly important in the rapidly rotating stars. However, our results of the trend of the cycle variability and the grand minima are not expected to change with many details of the model (e.g. type of nonlinearity, stochastic fluctuations, and turbulent transport) because they depend on the amount of supercriticality of the model. It is obvious to accept that the dynamo supercriticality decreases with the decrease of the rotation rate of the star (mainly due to the decrease of α). Extended grand minima are easy to produce when the dynamo is near critical (Kitchatinov & Olemskoy 2010; Vashishth et al. 2021). Also, observations hint that the solar dynamo (which produce grand minima) is operating near the critical transition (Rengarajan 1984; Metcalfe, Egeland & van Saders 2016). Furthermore, observing some robust results (mainly the increase of the number and duration of grand minima with the increase of rotation period) in all the models having different parameters, we can have some confidence that our results will be validated in the more realistic stellar dynamo models and observations.
ACKNOWLEDGEMENTS
We express our gratitude to Aparup Ghosh for providing us the code for computing Hurst exponent and Ricky Egeland for encouraging us to work on this problem. Additionally, we would like to acknowledge the contributions of Allan Sacha Brun and the anonymous referee, who provided valuable feedback and raised insightful questions, which enhances the quality of the paper. Financial Support from the Department of Science and Technology (SERB/DST), India, through the Ramanujan fellowship (project number SB/S2/RJN-017/2018) awarded to BBK is acknowledged. VV acknowledges the financial support from the DST through INSPIRE fellowship. LK acknowledges financial support from the Ministry of Science and High Education of the Russian Federation.
DATA AVAILABILITY
In our work, the dynamo calculations are done using a freely available code: Surya (Nandy & Choudhuri 2002; Chatterjee, Nandy & Choudhuri 2004). Data from our dynamo models and the analysis codes can be shared upon a reasonable request.
