A STOCHASTICALLY FORCED TIME DELAY SOLAR DYNAMO MODEL: SELF-CONSISTENT RECOVERY FROM A MAUNDER-LIKE GRAND MINIMUM NECESSITATES A MEAN-FIELD ALPHA EFFECT

Sunspots, which are strongly magnetized regions, play a key role in governing the activity of the Sun. The number of sunspots observed on the solar surface waxes and wanes with time generating the 11 yr solar cycle. While there is a small variation in this periodicity, fluctuations in the amplitude of the solar cycle are large. Extreme fluctuations are manifest in grand maxima episodes, when the cycle amplitudes are much higher than normal, and grand minima episodes, when the cycle amplitudes fall drastically, even leading to the disappearance of sunspots for an extended period of time. The most striking evidence of such a minimum in the recorded history of sunspot numbers is the so-called Maunder minimum between 1645 and 1715 AD (Eddy 1988). The lack of sunspots during this period is statistically well-proven and is not due to the lack of observations, which covered 68% of the days during this period (Hoyt & Schatten 1996). The occurrence of these solar activity extremes is correlated with temperature records over a millennium scale (Usoskin et al. 2005); the solar Maunder minimum coincided with the most severe part of the Little Ice Age—a period of global cooling on Earth.

Over the last decade, solar activity reconstructions, based on cosmogenic isotopes and geomagnetic data (Usoskin et al. 2000, 2003, 2007; Miyahara et al. 2004; Steinhilber et al. 2010; Lockwood & Owens 2011), which are indirect proxies for probing long-term solar activity, have brought to the fore various properties of these grand minima episodes. These observations show that there have been many such activity minima in the past; however, the solar cycle has recovered every time and regained normal activity levels. There is some evidence for persistent, but very weak, amplitude cycles during the Maunder minimum and a slow strengthening of cycle amplitudes to normal levels during the recovery phase. While the general perception was that the onset of the Maunder minimum was sudden, a recent reconstruction based on historical sunspot records has challenged that notion indicating that the onset phase of the minimum may have been gradual (Vaquero et al. 2011).

A magnetohydrodynamic (MHD) dynamo mechanism involving interactions of plasma flows and magnetic fields drives the solar cycle. Our understanding of the solar dynamo, see, e.g., the reviews by Ossendrijver (Ossendrijver 2003) and Charbonneau (Charbonneau 2010), is based on the generation and recycling of the toroidal and poloidal components of the Sun's magnetic field. The toroidal magnetic field is produced by the stretching of poloidal field lines by differential rotation—a process called the Ω effect (Parker 1955). It is thought that this process occurs all throughout the solar convection zone (SCZ), but is intensified near the base of the SCZ, where the upper part of the tachocline (a region of strong radial gradient in the rotation) and overshoot layer (which is stable to convection) offers an ideal location for toroidal field amplification and storage. Sufficiently strong toroidal flux tubes are magnetically buoyant and erupt radially outwards producing sunspots where they intersect the solar surface.

For the dynamo to function, the poloidal component has to be regenerated back from the toroidal component, a step for which diverse propositions exist. The first such proposition invoked helical turbulent convection as a means of twisting rising toroidal flux tubes to regenerate the poloidal component (a mechanism traditionally known as the mean-field α effect; Parker 1955). Numerous dynamo models based on the mean-field α effect were constructed and such models enjoyed a long run as the leading contender for explaining the origin of the solar cycle (Charbonneau 2010). However, subsequent simulations of the dynamics of buoyant toroidal flux tubes and observational constraints set by the tilt angle distribution of sunspots pointed out that the toroidal magnetic field at the base of the SCZ must be as high as 10⁵ G (D'silva & Choudhuri 1993; Fan et al. 1993; Caligari et al. 1995); such strong toroidal flux tubes being one order of magnitude stronger than the equipartition magnetic field in the SCZ would render the mean-field α effect ineffective. This consideration revived interest in an alternative mechanism of poloidal field production based on the flux transport mediated decay and dispersal of tilted bipolar sunspot pairs in the near-surface layers (Babcock 1961; Leighton 1969), hereby, referred to as the Babcock–Leighton mechanism.

In the last couple of decades, multiple dynamo models have been based on this idea (Durney 1997; Dikpati & Charbonneau 1999; Nandy & Choudhuri 2002; Chatterjee et al. 2004; Muñoz-Jaramillo et al. 2009) and have successfully reproduced many nuances of the solar cycle. Some (Tobias et al. 2006; Bushby & Tobias 2007; Cattaneo & Hughes 2009) have criticized the usage of such mean-field dynamo models to predict the solar cycle, however, it should be noted that recent studies (Simard et al. 2013; Dube & Charbonneau 2013) indicate that if input profiles are extracted from three-dimensional full MHD simulations and fed into two-dimensional mean-field dynamo models, they are capable of producing qualitatively similar solutions to those found in full MHD simulations. The major advantage of the mean-field dynamo framework is that it allows for much faster integration times compared to the full MHD simulations and are therefore computationally efficient as well as physically transparent. Recent observations also lend strong support to the Babcock–Leighton mechanism (Dasi-Espuig et al. 2010; Muñoz-Jaramillo et al. 2013) and this is now believed to be the dominant source for the Sun's poloidal field. Surface transport models (Wang et al. 1989; van Ballegooijen et al. 1998) also provide theoretical evidence that this mechanism is in fact operating in the solar surface. Randomness or stochastic fluctuations in the Babcock–Leighton poloidal field generation mechanism is an established method for exploring variability in solar cycle amplitudes (Charbonneau & Dikpati 2000; Charbonneau et al. 2004, 2005; Passos & Lopes 2011; Passos et al. 2012; Choudhuri & Karak 2012) as are deterministic or non-linear feedback mechanisms (Wilmot-Smith et al. 2005; Jouve et al. 2010). Stochastic fluctuations within the dynamo framework are physically motivated from the random buffeting that a rising magnetic flux tube endures during its ascent through the turbulent convection zone and from the observed scatter around the mean (Joy's law) distribution of tilt angles. It is to be noted that similar fluctuations are to be expected in the mean-field α effect as well (Hoyng 1988) and such phenomena can be explored within the framework of truncated mean-field dynamo models (Yoshimura 1975).

Since the two source layers for toroidal field generation (the Ω effect) and poloidal field regeneration (the α effect) are spatially segregated in the SCZ, there must be effective communication to complete the dynamo loop. Magnetic buoyancy efficiently transports the toroidal field from the bottom of the convection zone to the solar surface. On the other hand, meridional circulation, turbulent diffusion, and turbulent pumping share the role of transporting the poloidal flux from the surface back to the solar interior (Karak & Nandy 2012) where the toroidal field of the next cycle is generated, thus keeping the cycles going. Therefore, there is a time delay built into the system due to the finite time required for transporting magnetic fluxes from one source region to another within the SCZ.

Based on delay differential equations and the introduction of randomness on the poloidal field source, here we construct a novel, stochastically forced, non-linear time delay dynamo model for the solar cycle to explore long-term solar activity variations. We particularly focus our investigations on the recovery from grand minima phases and demonstrate that the Babcock–Leighton mechanism alone,which is believed to be the dominant source for the poloidal field, cannot restart the solar cycle once it settles into a prolonged grand minimum. The presence of an additional poloidal field source capable of working on weak magnetic fields, such as the mean-field α effect, is necessary for recovering the solar cycle.

The model is an extension of the low-order time delay dynamo equations explored in an earlier study involving one of us (Wilmot-Smith et al. 2006). This model was derived considering only the source and dissipative mechanisms in the dynamo process. All space dependent terms were removed and instead the physical effect of flux transport through space was captured through the explicit introduction of time delays in the system of equations.

The time delay dynamo equations are given by

$$\begin{eqnarray} &&\frac{dB_\phi (t)}{dt} = \frac{\omega }{L} A(t-T_0)-\frac{B_\phi (t)}{\tau } \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\frac{dA(t)}{dt} = \alpha _0 {f_1} (B_\phi (t-T_1)) B_\phi (t-T_1) -\frac{A(t)}{\tau },\end{eqnarray} \tag{ 2 }$$

where B_ϕ represents the toroidal field strength and A represents the poloidal field strength. The evolution of each magnetic component is due to the interplay of the source and dissipative terms in the system. In the toroidal field evolution equation, ω is the difference in rotation rate across the SCZ and L is the length of SCZ. The dissipative term is governed by turbulent diffusion, characterized by the diffusion timescale (τ). The parameter T₀ is the time delay for the conversion of the poloidal field into the toroidal field and is justified by the finite time that the meridional circulation or turbulent pumping takes to transport the poloidal magnetic flux from the surface layers to the tachocline. T₁ is the time delay for the conversion of the toroidal field into the poloidal field and accounts for the buoyant rise time of toroidal flux tubes through the SCZ. The meridional circulation timescale is about 10 yr for a peak flow speed of 20 m s⁻¹ (from mid-latitudes at near-surface layers to mid-latitudes above the convection zone base; see Yeates et al. 2008 for a detailed discussion on meridional flow timescales in the SCZ). Another dominant flux transport mechanism for the downward transport of the magnetic field could be turbulent flux pumping with a timescale of about one year (with a relatively high pumping speed of 5 m s⁻¹). The buoyant rise time of flux tubes from the SCZ base to surface is about three months (assuming the rise timescale is of the order of the Alfvénic timescale, which is in general agreement with simulations; see also Fan et al. 1993). As the magnetic buoyancy timescale is much shorter compared to the meridional circulation (or turbulent diffusion or flux pumping) timescale, we assume T₁ ≪ T₀. Since it is not clear which is the most dominant flux transport mechanism, meridional circulation or turbulent pumping, we explore our model in two different regimes of operation to test for robustness. In one setup, we consider T₀ = 4T₁ (if T₀ corresponds to pumping timescale) and T₀ = 40T₁ (if T₀ corresponds to meridional circulation timescale). This model setup mimics spatial separation between two source layers in the Sun's convection zone and the role of magnetic flux transport between them, and is therefore physically motivated. On the other hand, due to its nature, this model is amenable to long time-integration without being computationally expensive.

To account for quenching of the Babcock–Leighton poloidal source, α, we take a general form of α, i.e, α = α₀f₁, where α₀ is the amplitude of the α effect and f₁ is the quenching factor approximated here by a nonlinear function

$$\begin{eqnarray} f_1=& &\frac{\left[1+{\rm erf}\left(B^2_{\phi }(t-T_1)-B^2_{{\rm min}} \right) \right]}{2} \nonumber \\ & &\cdot \frac{\left[1-{\rm erf}\left(B^2_{\phi }(t-T_1)-B^2_{{\rm max}}\right) \right]}{2}. \end{eqnarray} \tag{ 3 }$$

Figure 1 depicts this quenching function, constructed with the motivation that only flux tubes with field strength above B_min (and not below) can buoyantly rise up to the solar surface and contribute to the Babcock–Leighton poloidal field source, i.e., sunspots (Parker 1955) and that flux tubes stronger than B_max erupt without any tilt, therefore quenching the poloidal source (D'silva & Choudhuri 1993; Fan et al. 1993). Accounting for these lower and upper operating thresholds for the Babcock–Leighton poloidal source is fundamentally important for the dynamics.

Zoom In Zoom Out Reset image size

Our aim here is to explore the impact of stochastic fluctuations in this time delay solar dynamo model. For α = α₀, we get a strictly periodic solution. In order to introduce stochastic fluctuations, we redefine α as

$$\begin{equation} \alpha =\alpha _0 \left[1+\frac{\delta }{100} \sigma (t,\tau _{{\rm cor}})\right], \end{equation} \tag{ 4 }$$

where σ(t, τ_cor) is a uniform random function lying in the range [+1, −1], changing values at a coherence time, τ_cor. Statistical fluctuations are characterized by δ and τ_cor, which correspond to percentile level of fluctuation and coherence time correspondingly. Figure 2 shows a typical α fluctuation generated by our random number generation program. Stochastic variations in the Babcock–Leighton α coefficient are natural because they arise from the cumulative effect of a finite number of discrete flux emergences, i.e., active region eruptions, all with various degrees of tilt randomly scattered around a mean Joy's law distribution.

Zoom In Zoom Out Reset image size

In this system the dynamo number (N_D = α₀ωτ²/L) is defined as the ratio between the source and dissipative terms, which is a measure of the efficiency of the dynamo mechanism. The product of the source terms is |α₀ω/L| while that of the dissipative terms is 1/τ². In terms of physical parameters, the expected diffusion timescale (L²/η) in the SCZ is 13.8 yr for a typical diffusivity of 10¹² cm² s⁻¹, implying that the dissipative term (1/τ²) is of the order of 10⁻¹⁸ s⁻². Now, if we take the value of ω as the difference in rotation rate across the SCZ in nHz (as measured; for details, see Howe 2009), L as the length of the SCZ and α₀ as 1 ms⁻¹ then the source term, |α₀ω/L|, is of the same order as the dissipative term and the dynamo number can be made higher than unity by slightly adjusting the α coefficient. In fact, if the tachocline is considered as the interface across which flux transport is occurring, then the dynamo number becomes even greater as the radial differential rotation is about the same while the length scale reduces further. In this model, we always take the value of |α₀ω/L| (source term) to be greater than 1/τ² (decay term), and set the magnitude of |ω/L| and |α₀| in a way such that the strength of the toroidal field is greater than the strength of the poloidal field (as suggested by observations). In summary, keeping all of the other physically motivated parameters fixed, the dynamo number can be varied by adjusting the value of α₀. Since B_min corresponds to the equipartition field strength (on the order of 10⁴ G) above which magnetic flux tubes become buoyant while B_max is on the order of 10⁵ G (above which flux tubes emerge without any tilt, thus shutting off the Babcock–Leighton source; D'silva & Choudhuri 1993), we take the ratio of B_max/B_min as 7 for all of our calculations. Here, we explore our low order time delay model in two parameter space regimes to test for robustness. In the first case, we fix the parameters as τ = 15, B_min = 1, B_max = 7, T₀ = 4T₁, T₁ = 0.5, and ω/L = −0.34, while in the second case, we take τ = 25, B_min = 1, B_max = 7, T₀ = 40T₁, T₁ = 0.5, and ω/L = −0.102. Initial conditions are taken to be (B_min + B_max)/2 for both A and B_ϕ. Our choice of parameters ensures that in both cases the diffusive timescale is much higher than flux transport timescales (τ > T₀  +  T₁). The simulations are robust over a range of negative N_D values; for a detailed parameter space study of the underlying model without stochastic fluctuations, please refer to Wilmot-Smith et al. (2006). Below, we present the results of our stochastically forced dynamo simulations focussing on entry and exit from grand minima episodes.

We first perform simulations without the lower operating threshold in the Babcock–Leighton α effect (setting B_min = 0) in Equation (3). A majority of Babcock–Leighton dynamo models, including many that have explored the dynamics of grand minima do not use this lower operating threshold. As already known (Charbonneau & Dikpati 2000; Choudhuri & Karak 2012), we find that the Babcock–Leighton dynamo with this setup generate cycles of varying amplitudes, including episodes of higher than average activity levels (grand maxima) and occasional episodes of very low amplitude cycles reminiscent of Maunder-like grand minima (Figure 3, upper panel; Figure 4, upper panel). When we do switch on the lower operating threshold, however, we find that the Babcock–Leighton dynamo is unable to recover once it settles into a grand minimum (Figure 3, lower panel; Figure 4, lower panel). This striking result can be explained invoking the underlying physics of the solar cycle. When a series of poloidal field fluctuations lead to a decline in the toroidal field amplitude below the threshold necessary for magnetic buoyancy to operate (with a consequent failure of sunspots to form), the Babcock–Leighton poloidal field source, which relies on bipolar sunspot eruptions, completely switches off resulting in a catastrophic quenching of the solar cycle. Earlier simulations, which did not include the lower quenching, missed out on this physics because even very weak magnetic fields, which in reality could never have produced sunspots, continued to (unphysically) contribute to poloidal field creation. Earlier, it was shown that the lower threshold due to magnetic buoyancy plays a crucial amplitude-limiting role in the Babcock–Leighton solar cycle (Nandy 2002) and this study indicates that this should be accounted for in all Babcock–Leighton solar dynamo models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To circumvent this problem faced by the stochastically forced Babcock–Leighton dynamo, we explicitly test an idea (Nandy 2012) for the recovery of the solar cycle based on an additional poloidal source effective on weak toroidal fields. Since the tachocline is the seat of strong toroidal field, any weak field α that is only effective on a sub-equipartition strength field will get quenched there. Thus, this α effect must reside above the base of the SCZ (Parker 1993) in a layer away from the strongest toroidal fields. Motivated by this, we devise a new system of dynamo equations governed by

$$\begin{eqnarray} &&\frac{dB_\phi (t)}{dt} = \frac{\omega }{L} A(t-T_0)-\frac{B_\phi (t)}{\tau } \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \frac{dA(t)}{dt} =& & \alpha _0 f_1(B_\phi (t-T_1)) B_\phi (t-T_1) \nonumber \\ & &+ \alpha _{{\rm mf}} f_2(B_\phi (t-T_2)) B_\phi (t-T_2) -\frac{A(t)}{\tau } \end{eqnarray} \tag{ 6 }$$

where f₂, the quenching function for the weak field poloidal source, α_mf, is shown in Figure 1 and is parameterized by

$$\begin{equation} f_2 = \frac{{\rm erf}c \left(B^2_\phi (t - T_2)-B^2_{{\rm eq}}\right)}{2}. \end{equation} \tag{ 7 }$$

Taking B_eq = 1 ensures that the weak field source term gets quenched at or below the lower operating threshold for the Babcock–Leighton α and the former, therefore, can be interpreted to be the mean-field α effect. In Equation (6), the time delay, T₂, is the time necessary for the toroidal field to enter the source region where the additional, weak-field α effect is located.

If T₂ = 0, i.e., the generation layer of the additional α effect is coincident with the Ω effect (layer), then we find that the stochastically forced dynamo again fails to recover from a grand minimum. This is reminiscent of the original motivation behind the introduction of the interface dynamo idea with spatially segregated source regions (Parker 1993). However, if T₂ is finite and T₁ > T₂ (i.e., there is some segregation between the Ω effect toroidal source, the additional weak-field α and the Babcock–Leighton α), we find that the solar cycle can recover from grand minima-like episodes in a robust manner. Figures 5 and 6 depict such solutions (for two different sets of parameter), where we explicitly demonstrate self-consistent entry and exit from grand minima-like episodes. We note that this recovery of the solar cycle from grand minima-like episodes is possible with or without fluctuations in the additional, weak field poloidal source term α_mf.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, we have constructed a new model of the solar dynamo for exploring solar cycle fluctuations based on a system of stochastically forced, non-linear, delay differential equations. Utilizing this model for long-term simulations, we have explicitly demonstrated that the currently favored mechanism for solar poloidal field production, the Babcock–Leighton mechanism, alone, cannot recover the solar cycle from a grand minimum. We have also demonstrated that an additional, mean-field-like α effect capable of working on weaker fields is necessary for self-consistent entry and exit of the solar cycle from grand minima episodes. We have demonstrated that our results and conclusions hold over two very diverse regimes of parameter choices. We note that simulations motivated from this current study and based on a spatially extended dynamo model in a solar-like geometry supports the results from this mathematical time delay model (Passos et al. 2014). Taken together, these strengthen the conclusion that a mean-field-like α effect effective on weak toroidal fields must be functional in the Sun's convection zone and that this is vitally important for the solar cycle, even if the dominant contribution to the poloidal field comes from the Babcock–Leighton mechanism during normal activity phases.

We are grateful to the Ministry of Human Resource Development, the Council for Scientific and Industrial Research, the University Grants Commission, and the Ramanujan Fellowship award of the Department of Science and Technology of the Government of India for supporting this research. D.P. acknowledges support from the Fundação para a Ciência e Tecnologia grant SFRH/BPD/68409/2010, CENTRA, and the University of the Algarve and is also grateful to IISER Kolkata for hosting him during the performance of this research.