Stellar mid-life crisis: subcritical magnetic dynamos of solar-like stars and the breakdown of gyrochronology

ABSTRACT

Recent observations have revealed the surprising breakdown of stellar gyrochronology relations at about the age of the Sun hinting the middle-aged, solar-like stars transition to a magnetically inactive future. We provide a theoretical basis for these intriguing observations inspired by simulations with a mathematical-dynamo model that can explore long-term solar cycle fluctuations. We reproduce the observed bimodal distribution of sunspot numbers, but only for subcritical dynamos. Based on a bifurcation analysis, we argue that the ageing of solar-like stars makes the magnetically weak dynamo regime readily accessible. Weak magnetic field production in this regime compromises wind-driven angular momentum losses, thus disrupting the hegemony of magnetic braking on stellar rotational spin-down. This hypothesis of subcritical magnetic dynamos of solar-like stars provides a self-consistent, unifying physical basis for a diversity of solar-stellar phenomena such as why stars beyond their mid-life do not spin-down as fast as in their youth, the break-down of stellar gyrochronology relations, the observed bimodal distribution of long-term sunspot observations, and recent findings suggesting that the Sun may be transitioning to a magnetically inactive future.

dynamo, MHD, sunspots, stars: evolution

1 INTRODUCTION

The widely used diagnostic of stellar rotation rate versus age relationship – that is the basis of gyrochronology – has recently been found to be unreliable after stellar mid-life (van Saders et al. 2016). The standard understanding of angular momentum loss mediated braking of rotation due to the magnetically coupled stellar wind cannot explain this. It has been reported that the Sun is currently going through a weakened braking phase at about its mid-life of 4.5 billion years (van Saders et al. 2016). What causes this weakened magnetic braking is an intriguing question. The latest measurements confirm that the Sun is less magnetically active than similar stars (Reinhold et al. 2020b). This has spurred debates and controversy on whether the Sun is already transitioning to a magnetically inactive future (Metcalfe & van Saders 2020; Reinhold et al. 2020a); studies exploring rotation period, magnetic cycle period relationships support this hypothesis (Metcalfe, Egeland & van Saders 2016; Metcalfe & van Saders 2017). Since stellar magnetic activity is ultimately a product of the magnetohydrodynamic (MHD) dynamo mechanism, a causal basis could exist between the nature and evolution of stellar dynamo action and these observations. Inferring the nature of magnetic dynamos in solar-like stars at about their mid-life, therefore, could be the Rosetta stone to explain these apparent disparate observations.

With an aim to characterize solar-like dynamos at around their mid-life – i.e. the current age of the Sun we look at long-term solar observations. These reveal intriguing phases of reduced sunspot activity that are termed as grand minima episodes, the most recent being the Maunder minimum between 1645 and 1715 AD (Eddy 1976). Several long-term solar activity reconstructions based on the cosmogenic isotopes like Be¹⁰ in ice cores and C¹⁴ in tree rings suggest that such grand minima have also occurred in the past (Usoskin, Solanki & Kovaltsov 2007); reconstructed data reveal 20 grand minima in the last 9000 yr (Usoskin et al. 2016). Such grand minima have recently been confirmed as a special mode of solar dynamo operation, distinct from the regular activity cycle mode (Usoskin et al. 2014).

Similar low activity phases have been observed in other Sun-like stars, which have been dubbed as ‘Maunder-minima stars’ (Baliunas et al. 1995); the nature of the magnetic dynamos in these stars is hotly debated. A crucial observation is that the main sequence age of these low activity ‘Maunder-minima stars’ are typically similar to, or higher, compared to the Sun. This is thought to be due to the rotational spin-down through steady angular momentum losses mediated via stellar winds and consequent reduction in the efficiency of the dynamo mechanism (Nandy 2004; Brun et al. 2015). It has been speculated that this might be related to a drastic change in the nature of the solar–stellar dynamos at the current solar age (Metcalfe et al. 2016; Metcalfe & van Saders 2017; Metcalfe & Egeland 2019).

Are the wandering of the magnetic activity of the Sun between the two distinct modes of low and regular activity cycles, the occurrence of ‘Maunder minima stars’, and the observed break-down of the rotation rate–age relationship connected? Based on a reduced dynamo model that imbibes the spiritual essence of Babcock–Leighton (BL) and turbulent dynamos by utilizing stochastically forced, non-linear delay differential equations, we show these diverse phenomena can be explained on the basis of the bimodality of the underlying dynamo mechanism. Based on this analysis, we argue that a stellar evolution mediated shift in the dynamo efficiency enables a transition across supercritical to subcritical states that makes the full dynamical regime – including very weak activity states – accessible to an ageing star.

The reconstructed sunspot number suggests that the occurrence of grand minima is not a cyclic phenomenon but an outcome of an underlying generation mechanism that is inherently stochastic and/or chaotic (Usoskin et al. 2007). In a chaotic dynamo, the long time-scale ‘supermodulation’ could be related to symmetry changes in the magnetic fields (Weiss & Tobias 2016). The current understanding of the large-scale magnetic field generation in Sun-like stars is based on the recycling of toroidal and poloidal magnetic field components (Parker 1955; Charbonneau 2020). The toroidal field is generated from the poloidal field due to the stretching by differential rotation. For the current Sun, observations and data-driven model reconstructions of solar activity (Dasi-Espuig et al. 2010; Cameron & Schüssler 2015; Bhowmik & Nandy 2018) demonstrate that the dominant mechanism for poloidal field creation is the so-called BL mechanism, i.e. the decay and dispersal of bipolar sunspot pairs which contribute to a large-scale dipolar field (Babcock 1961; Leighton 1969). This mechanism relies on the existence of strong toroidal flux tubes in the convection zone of the Sun, which is not the case during grand minima episodes. It has been demonstrated recently that a mean-field dynamo α-effect capable of working on weaker fields can ‘rescue’ the solar cycle from a grand minimum (Hazra, Passos & Nandy 2014; Passos et al. 2014) and therefore, must be operational in stellar convection zones. Direct numerical simulations utilizing the full MHD system of equations or spatially extended mean-field dynamo models are too expensive computationally for the necessary long-term simulations to explore the range of possible dynamical solutions that a solar-like dynamo mechanism may exhibit over long-term. For our numerical bifurcation analysis, we adapt a reduced BL dynamo model based on delay differential equations that have been found very useful in providing insights on long-term solar cycle behaviour (Wilmot-Smith et al. 2006; Hazra et al. 2014).

In BL dynamos, the sources of the toroidal and the poloidal fields are spatially segregated as they act across a distributed domain in the solar convection zone (SCZ). Such spatial segregation necessitates a finite time for flux transport across the SCZ, introducing time delays in the communication between the two source layers (Jouve, Proctor & Lesur 2010b). These time delays have been shown to create a dynamical memory even in stochastically forced dynamo systems which allow a predictive window for the activity (Yeates, Nandy & Mackay 2008; Karak & Nandy 2012; Bhowmik & Nandy 2018).

The evolution of the toroidal field (B_ϕ) and the vector potential (A) for the poloidal field in the reduced BL dynamo model are governed by the following delay differential equations

$$\begin{eqnarray*} &&\frac{\mathrm{ d}B_\phi (t)}{\mathrm{ d}t} = \frac{\omega }{L} A(t-T_0)-\frac{B_\phi (t)}{\tau }, \end{eqnarray*}$$

(1)

 

$$\begin{eqnarray*} &&\frac{\mathrm{ d}A(t)}{\mathrm{ d}t} = \alpha _{\mathrm{ BL}} \cdot {f} \left(B_\phi \left(t-T_1\right)\right) \cdot B_\phi (t-T_1) -\frac{A(t)}{\tau }, \end{eqnarray*}$$

(2)

where T₀ and T₁ are the time delays in the conversion of the poloidal flux to the toroidal one and vice versa. The interplay amongst meridional circulation, turbulent diffusion, and turbulent pumping mediated flux transport (Hazra & Nandy 2016) determine T₀. While there have been suggestions that meridional flow variations can explain the observed stellar activity amplitude–period relationship (Nandy 2004), numerical simulations exploring meridional flow dependence on the stellar rotation rates indicate the necessity of more complex flow profiles or additional mechanisms to explain this (Jouve, Brown & Brun 2010a). Thus we reiterate that T₀ must be interpreted as a generic flux transport time-scale, not necessarily determined by stellar meridional flows alone. The turbulent dissipation time-scale of the field is τ, ω is the difference in rotation rates across the SCZ of length L, and α_BL represents the BL source term for the poloidal field. The quenching function f constrains the poloidal source to work within a finite range of toroidal field strength with a non-zero lower bound. This is motivated from the perspective that the weak toroidal flux tubes are shredded by turbulence and cannot form the sunspots and very strong flux tubes are immune to the Coriolis force, which causes them to emerge at the solar surface without any tilt. This effectively quenches the BL source (for further details and motivation, see Wilmot-Smith et al. 2006; Hazra et al. 2014. The efficiency of the dynamo is measured by a non-dimensional number called the dynamo number, which is defined as N_D = ωα_BLτ²/L.

2 RESULTS AND DISCUSSION

Long-term simulations with this dynamo model show hysteresis (Fig.1 left-hand panel). There is a hint of similar behaviour in earlier studies with non-linear mean-field dynamos and three-dimensional MHD simulations of turbulent dynamos (Kitchatinov & Olemskoy 2010; Karak, Kitchatinov & Brandenburg 2015b; Oliveira et al. 2021). However, these computationally expensive numerical simulations precluded an adequate analysis of the distribution–function of a large number of magnetic cycles. The results from our model (Fig. 1 left-hand panel) substantiate the proposition and demonstrate that hysteresis engenders two distinct modes of dynamo activity. The results are qualitatively robust with respect to the choice of parameters within the regimes established in Wilmot-Smith et al. (2006) and Hazra et al. (2014).

Figure 1.

(Left-hand panel) Bifurcation diagram of toroidal field amplitude. The middle red curve (unstable limit cycles) between the upward and the downward arrow separates the basins of attraction for the oscillating solutions (stable limit cycles) and the decaying solutions (stable fixed points). The dotted rectangular region at the bottom left is zoomed in and is shown in the middle panel. (Middle panel) The region within the dotted vertical lines at |N_D| = 9.1 and |N_D| = 16.9 is the dynamo working region with 30 per cent fluctuation in the dynamo number |N_D^mean| = 13. The dashed–dotted and dashed horizontal lines are the sunspot formation threshold and the dynamo operating threshold (below which the BL process alone cannot recover the decaying cycles). The yellow shaded region represents the basin of attraction for the oscillatory solutions. The numbers (66, 68, 69, 80, 83) are cycle numbers counted from the beginning of the simulation. The fluctuating dynamo numbers (depicted by black filled circles) are noted at each cycle maximum of the time series (shown in the inset). (Right-hand panel) Typical onset of (and exit from) a grand minimum phase when an additive magnetic noise is employed with a typical root-mean-squared (rms) value of about 17 per cent of the rms value of the BL source in equation (3). The parameters used are ω/L = −0.34, τ = 15, T₀ = 2, T₁ = 0.5 yr.

The dynamo cannot self-consistently recover from a grand minimum once it is nudged into such a phase by stochastic fluctuations in the BL source term (Fig. 1 middle panel). Our extensive parameter space study shows that the BL dynamo shut-down threshold (dashed horizontal line in the inset of the middle panel of Fig. 1) shifts upward to engulf more cycles into grand minima episodes when the time-scale T₀ is increased in the diffusion-dominated regime (or decreased in the advection-dominated regime) as defined in Wilmot-Smith et al. (2006).

We find that the addition of a ‘magnetic noise’ allows the dynamo to recover from grand minima episodes as well as sustain hysteresis. Incorporation of additive noise, ϵ(t), modifies equation (2) to the form

$$\begin{eqnarray*} \frac{\mathrm{ d}A(t)}{\mathrm{ d}t} = \alpha _{\mathrm{ BL}} \cdot {f} \left(B_\phi \left(t-T_1\right)\right) \cdot B_\phi (t-T_1) -\frac{A(t)}{\tau } + \epsilon (t), \end{eqnarray*}$$

(3)

where ϵ(t) is uniform and white in time (Charbonneau 2001; Charbonneau, Blais-Laurier & St-Jean 2004; Charbonneau, Beaubien & St-Jean 2007) and has a zero mean. Since the noise has zero mean, it sometimes builds up the poloidal field while at other times, reduces the existing dipolar field of the Sun. Magnetic noise definitely abounds in the solar convection zone, especially near the sub-surface layers of the Sun where the convection is much more turbulent. Small-scale dynamo action in this region (Meneguzzi & Pouquet 1986; Nordlund et al. 1992) produces magnetic field that is highly uncorrelated in space and time. Such noise continuously regenerates and replenishes the magnetic field in the upper layers of the Sun (Cattaneo, Emonet & Weiss 2003; Hagenaar, Schrijver & Title 2003; Charbonneau et al. 2004) and our results (Fig. 1 right-hand panel) indicate this could be enough to recover solar–stellar dynamos from a grand minima over a growth time that is dependent on the level of imposed noise.

A time series of the smoothed cycle-averaged sunspot number from an additive-noise incorporated BL dynamo model is shown in Fig. 2, upper panel. The smoothing has been done using the same procedure used by Usoskin et al. (2016), i.e. we first average the monthly sunspot number (SSN) over the usual cycle period and apply the Gleissberg low pass 1-2-2-2-1 filter. The shaded regions below the horizontal line signify grand minima, which are defined as epochs during which the smoothed SSN falls below 50 per cent of its average value for at least three consecutive decades. We observe 13–18 grand minima (18 in Fig. 2) in different realizations of stochastic fluctuations in a typical 9000 yr simulation, which is in close agreement with the reconstructed data identifying 20 grand minima in the last nine millennia; this is of course sensitive to model parameter choices (Choudhuri & Karak 2012).

Figure 2.

(Upper panel) Temporal variation of the smoothed SSN from a 9000 yr simulation. The shaded regions below the horizontal line represent grand minima episodes. (Lower panel) Variation of monthly SSN for a selected duration of 2100 yr. The monthly SSN is assumed to be a measure of the toroidal field energy (i.e. |$\mathrm{ SSN} = B_\phi ^{2}$|⁠). The (broken) horizontal lines depict the extent of different grand minima epochs. The parameter used is |N_D^mean| = 13 and the level of fluctuations in |N_D| is |$30{{\ \rm per\ cent}}$|⁠. The rms value of the additive noise imposed here is about 24 per cent of the rms value of the Babcock-Leighton source. Other parameters of the model are set at ω/L = −0.34, τ = 15, T₀ = 2, T₁ = 0.5 yr.

What is fascinating is that the assimilation of magnetic noise in our model generates magnetic cycles that qualitatively reproduce the observed bimodal distribution of toroidal field energy; see Fig. 3 and for a comparison with observations, Fig. 3 in Usoskin et al. (2014). This implies that grand minima episodes cannot be considered as events that correspond to a tail of a single regular mode of solar activity, but it is in itself a separate state of dynamo operation. Most previous studies of grand minima relied on fluctuations in a single regular mode of activity cycles (Brandenburg & Spiegel 2008; Moss et al. 2008; Olemskoy, Choudhuri & Kitchatinov 2013) and thus could not explain the observed bimodal distribution. This bimodal nature in the probability density of cycle-averaged SSN is, in fact, due to the hysteresis present in the model (Kitchatinov & Nepomnyashchikh 2015). Hysteresis in a certain parameter regime emerges with the existence of bistable solutions, which emanate in our model due to the non-zero lower operating threshold of the BL process and the shorter poloidal-to-toroidal flux conversion time-scale. Numerical experiments with our dynamo model show that if we remove the lower operating threshold and increase the time delays simultaneously, hysteresis and the bimodal nature of the probability density of cycle-averaged sunspot number is no longer manifested.

Figure 3.

Bimodal distribution of probability density of toroidal magnetic energy obtained from the reduced dynamo model, which is gleaned from a 2 million years simulation run. The relevant dynamo parameters are |N_D^mean| = 13, level of fluctuations in |N_D| is |$20{{\ \rm per\ cent}}$| and the rms value of the additive noise is about |$17{{\ \rm per\ cent}}$| of the rms value of the BL source. Other parameters are fixed at ω/L = −0.34, τ = 15, T₀ = 2, T₁ = 0.5 yr.

Our model also qualitatively reproduces the phase space collapse of toroidal field; compare Fig. 4 upper panel here with fig. 10 in Lopes et al. (2014), which Passos & Lopes (2011) unsuccessfully attempted to achieve using their low dimensional model with stochastic fluctuations in α-effect alone. The collapse in the phase space implies that a grand minimum episode results in strong intermittency in magnetic cycles and that the weak intermittency (or moderate amplitude cycles) cannot be attributed to a grand minimum epoch. It is interesting to note in Fig. 4 (upper panel) that the polarity reversal of magnetic cycles, albeit far less regular, still occurs during the grand minimum when the BL source is dysfunctional. This is because the stochastically injected poloidal field due to the magnetic noise still provides a seed for differential rotation to generate toroidal fields. A combination of statistical fluctuations in the right direction can then nudge the weak toroidal field such that it exceeds the lower operating threshold and revive the large-scale BL dynamo that produces sunspots. The lower panel of Fig. 4 portrays the corresponding temporal variation of toroidal field energy depicting this dynamics.

Figure 4.

(Upper panel) Phase space collapse seen in a typical 400 year simulation when a grand minimum phase was observed. This is the same grand minimum shown in Fig. 2 (lower panel) within the dotted vertical lines at around time  = 4400 yr. The filled circle (square) represents the beginning (end) of the considered 400 yr. (Lower panel) Corresponding temporal variation of sunspot number in the simulation.

We observe the bimodal distribution of sunspot numbers only in the (significantly) subcritical regime of the dynamo. Earlier works in this context were inconclusive and conflicting (Karak et al. 2015b; Kitchatinov & Nepomnyashchikh 2017). The dynamo hysteresis in our model sustains coexisting oscillating and decaying solutions for similar dynamo parameters, which is not possible for supercritical dynamos. Thus, our theoretical analysis establishes that the observed bimodal distribution of decadal sunspot numbers is a direct manifestation of the subcritical nature of the current solar dynamo (at mid-life of its evolution as a star).

Over its lifetime, a Sun-like star spins-down due to angular momentum loss through the magnetically and thermally driven stellar winds. Ageing and lower rotation rates (i.e. higher rotation period P_Ω) increase the stellar Rossby number (R_o  = P_Ω/τ_c, where τ_c is the convective turn-over time) – a parameter used widely in observational studies of stellar activity. The efficiency of the dynamo generation of magnetic fields is represented by the dynamo number, which scales as N_D ∼ 1/R_o². Based on our earlier analysis, we postulate that the ageing of a star – with a concomitant increase in R_o and decrease in N_D – makes the subcritical dynamo regime with a bimodal activity distribution accessible around the mid-life of its evolution on the main sequence.

3 CONCLUSION

Our analysis is relevant only for stars which continue to exhibit solar-like rotation and thus not applicable for stars with anti-solar type differential rotation, see e.g. Karak et al. (2015a, 2020), Brandenburg & Giampapa (2018). A systemic change in global coronal magnetic field configuration (e.g. from closed to open loops) may also impact wind properties and the nature of magnetic braking which we do not address here.

Our findings have several important consequences. We establish that Maunder-like grand minima are a distinct mode of dynamo operation possible only in subcritical dynamo regimes. Stellar observations indicate the existence of low activity ‘Maunder minima stars’, which are similar or older than the Sun (Baliunas et al. 1995) and have high R_o. Our framework suggests these stars are in the subcritical regime, where accessing the low activity mode of the bimodal stellar dynamo becomes possible. Based on our analysis, we postulate that this is not due to a drastic change in the underlying dynamo mechanism itself, but merely the opening up of the possibility of existing in the low activity mode – wherein, very low dynamo output drastically reduces magnetically driven stellar winds and angular momentum loss rates. At about stellar mid-life, the magnetically weak subcritical states are readily attainable. These weak-magnetic epochs disrupt the hegemony of magnetic braking on stellar rotation. This argument is consistent with the suggestion by Metcalfe et al. (2020) – namely, magnetic cycles become weaker while still intermittently exhibiting large-scale dynamo action. This regime allows for the Rossby number (or equivalently, the dynamo number) to meander across the critical Rossby number for transition between different activity states or remain nearly constant, as hypothesized by Metcalfe et al. (2020) for 94 Aqr Aa.

To summarize, the subcritical solar magnetic dynamo hypothesis self-consistently explains the breakdown of stellar gyrochronology, the bimodal nature of the dynamo mechanism in the current Sun, solar cycle transition between grand minima and regular activity modes and the possible transition of the Sun to a less magnetically active future.

ACKNOWLEDGEMENTS

The authors thank Travis S. Metcalfe and Jennifer van Saders for useful discussions. BT wishes to acknowledge the support extended to him during a visiting fellowship at the Center of Excellence in Space Sciences India (CESSI) at the Indian Institute of Science Education and Research Kolkata. CESSI is funded by the Ministry of Education, Government of India under the Frontier Areas of Science and Technology Scheme.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.