© ESO, 2014

1. Introduction

For hundreds of years mankind has kept a close eye on the evolution of our star. Since the advent of the telescope, solar activity has been monitored systematically, creating one of the longest continuous observation programmes in the history of science (Owens 2013), the sunspot number measurements. These solar features continue to be used as the most common proxy for solar magnetic activity as their appearance on the Sun’s surface is thought to be connected with the build-up of large-scale, coherent magnetic fields in the solar convection zone (SCZ) mainly in the azimuthal direction, the so called toroidal magnetic field component. This field component is produced by a magnetohydrodynamic (MHD) dynamo that converts a portion of the kinetic energy of the solar differential rotation into magnetic energy (Parker 1955). This process takes place in a thin shear layer at the bottom of the SCZ, known as the tachocline, between the turbulent convection zone and the underlying stable overshoot layers. The most common theoretical framework used to explain the origin and evolution of the large-scale solar magnetic field is the mean-field dynamo theory (Steenbeck et al. 1966). Models based on this theory solve a simplified version of the MHD induction equation, usually in spherical-polar coordinates under some physically inspired simplifications such as axisymmetry and spatio-temporal averages of the magnetic field (see review by Charbonneau 2010). Considered by some as the most successful of these models, Babcock-Leighton (BL) flux transport dynamo models (Wang et al. 1991; Durney 1995; Dikpati & Charbonneau 1999; Nandy & Choudhuri 2002) have become very useful tools for explaining the main features of the solar cycle such as its periodicity, parity, the migration pattern of sunspots towards the equator and the phase difference between the toroidal and dipolar magnetic fields. These models incorporate the solar meridional circulation and emulate the effect of the decay of active regions for the production of the poloidal field based on the original ideas of Babcock and Leighton (Babcock & Babcock 1955; Leighton 1969).

For many decades, the importance of the BL surface mechanism was overshadowed by another type of poloidal field source, the mean-field α effect (or Parker’s α effect) in mean-field αΩ dynamos. In this type of dynamo model, the source for transforming a toroidal field into a poloidal field is the action of helical turbulence in the twisting of flux tubes as they buoyantly rise through the SCZ (Charbonneau 2010). The mean-field α effect is thought to be quenched by super-equipartition magnetic fields and to be effective only on relatively weaker fields. Works such as D’Silva & Choudhuri (1993) have raised questions regarding the effectiveness of this type of alpha effect by asserting that Joys’ Law distribution (tilt angle dependence on latitude of bipolar sunspot pairs) requires very strong, super-equipartition toroidal field strengths; see also Fan (2001). Moreover, recent observational evidence points out that the BL mechanism is the predominant source for the solar cycle (Dasi-Espuig et al. 2010). Following these developments the stakes seem to be loaded against the mean-field α effect mechanism. Nevertheless, recent, realistic 3D MHD simulations of the SCZ (Käpylä et al. 2006; Brown et al. 2010; Racine et al. 2011), show clear evidence that a mean-field-type α effect should exist in the bulk of the convection zone. The question is whether this mean-field α plays a dynamically important role in the solar cycle.

Important constraints on dynamo models have emerged in the last decade based on long-term reconstructions of solar activity using cosmogenic isotopes (Miyahara et al. 2006; Usoskin et al. 2007; Usoskin 2008). These reconstructions show that in addition to the typical variability observed in the 11-year solar cycle, there are extended periods of time in which the Sun behaves in a very unusual way. Solar magnetic activity seems to exhibit intermittency (fluctuations and episodes of stronger than normal and quieter than normal activity), as it is usually defined in the jargon of Dynamic Systems. This means that there are periods where the normal cyclic activity ceases to exist and the Sun enters a grand minimum period. During these grand minima, the Sun exhibits a very low number (or lack) of sunspots in its photosphere and it is believed that the overall magnetic field also becomes very weak. Thanks to recent data, e.g. Steinhilber et al. (2012), we now know that these periods are much more common than previously thought and therefore they should be considered as a constraint on any model of the solar cycle. The most famous of these quiescent periods of activity is called the Maunder minimum that took place between 1645 and 1715 (Maunder 1904; Eddy 1976). This period of low solar activity coincided with a period of very low temperatures in the Earth’s northern hemisphere during the last part of the Little Ice Age. This and other interesting correlations between solar activity and the Earth’s climate have been found in the past (Haigh 2007). A vitally important aspect of the grand minima enigma is to understand how the solar dynamo reawakens from these quiescent phases and the dynamics leading to, during, and exiting the phases.

Over the years some authors have focussed their efforts in modelling and understanding the mechanisms behind the observed short- and long-term variability (e.g. amplitude variations from cycle to cycle and grand minima-maxima) of the solar cycle. Two major approaches are usually used in these studies, one that looks for dynamical effects and a second that studies the impact of stochasticity on the system. The dynamical mechanism approach is based on non-linear interactions between the magnetic field and plasma flows. One of the most important ingredients in any dynamo model is the differential rotation. So, it is only natural to assume that perturbations in this strong, large-scale flow can induce fluctuations in activity. The so called Malkus-Proctor effect, i.e. the feedback of the large-scale magnetic field on the rotation is such an example, wherein, strong magnetic fields quench the very flow which sustains dynamo action. This results in a weakening of magnetic activity, subsequently to which the Reynolds stresses restore the differential rotation to its initial state, thereby resuming activity (Tobias 1997; Brooke et al. 2002; Bushby 2006). While this remains a plausible scenario for activity modulation, recent analysis of global 3D MHD simulations of solar convection that present dynamo action, indicates that modulation effects on the differential rotation could be much smaller than previously thought (Beaudoin et al. 2012). Other examples of dynamical processes are the modulation of the meridional circulation by the magnetic field as seen in dynamo models working in the non-kinematic regime (Rempel 2006; Passos et al. 2012) and the time delay dynamics that are inherent in the dynamo system with spatially segregated source-layers (Charbonneau et al. 2005; Wilmot-Smith et al. 2006; Jouve et al. 2010).

The modelling approach based on stochastic fluctuations is motivated by the fact that the solar dynamo resides in a very turbulent environment, i.e. the SCZ. So it is highly plausible that some of the intervening physical mechanisms are affected by “noise” imparted by random processes that are the basis of turbulent convection. This is usually modelled as stochastic fluctuations around an average typical value of some model’s ingredient such as the α effect (Hoyng 1988; Choudhuri 1992; Hoyng 1993, 1994; Charbonneau et al. 2004; Usoskin 2008) or the meridional circulation (Charbonneau & Dikpati 2000; Lopes & Passos 2009; Karak 2010). Even if noise is introduced in the system at correlation timescales that are much shorter than the solar cycle timescale they tend to induce modulations on timescales ranging from decades to centennia. These models are also very robust because they can handle a wide range of noise levels. Some of these stochastically forced models have been very successful in providing a simple explanation for the observed short-term solar variability, and a few of them also display grand minima-like episodes as well. For the interested reader, a comprehensive review of modulation mechanisms of the solar dynamo can be found in Tobias (2002).

Here, we base our study on a stochastically forced, kinematic, BL dynamo model. Inspired by the study of Hazra et al. (2014, henceforth Paper I), we explore the role of stochastically forced, dual poloidal field sources in the context of grand minima and maxima in solar activity. In Paper I the authors present a novel solution to how the dynamo re-emerges from a grand minima episode. In that paper they used a time-delay low-order dynamo model to explore what happens when two α effects work simultaneously, one effective on weak toroidal fields (mean-field α) and another effective on stronger fields (BL α). They demonstrate that the weak mean-field effect has a crucial role in driving the dynamo out of a grand minimum and postulate the usage of both an upper and lower threshold for the BL α in order to avoid fake, i.e physically incorrect, recovery from grand minima episodes. Here we follow this postulate, test the original idea, and also study additional dynamics associated with extreme solar activity with a spatially extended solar dynamo model with solar-like differential rotation and other physically inspired parameterizations. We introduce a secondary, mean-field α effect in the bulk of the convection zone in our BL flux transport model and explore the consequent dynamics in the stochastically forced regime. Our results support the main conclusions presented in Paper I and we explore the subject further. In the following section we discuss the base model and show how fluctuations in the surface BL α effect induce the dynamo to enter a grand minima-like state. In the subsequent section we add the secondary mean-field α effect to the model, discuss the modifications introduced to the parametrization of both α effects, and demonstrate scenarios wherein stochastic fluctuations trigger grand maxima and minima behaviour and self-consistent recovery from them. We also analyse the resultant butterfly diagrams and demonstrate that hemispheric coupling during quiescent phases can range from strong asymmetry to near simultaneous occurrence of quiescent phases in both the northern and southern solar hemispheres. We end with a concluding section discussing the implications of our results and putting them in a broader context.

2. The model

The axisymmetric dynamo Eqs. (1) and (2) are solved in the kinematic regime using a modified version of the Surya - code (Nandy & Choudhuri 2002; Chatterjee et al. 2004). This code solves the equations

where s = rsin(θ), v is the meridional flow, Ω the differential rotation, and α is a surface source term that emulates the BL mechanism. This numerical model assumes different magnetic diffusivities for the poloidal and toroidal field components, η_p and η_t, respectively. The use of two distinct diffusivity profiles in the Surya code is inspired by the reasoning that strong magnetic fields will suppress turbulence. Hence it is expected that in the strong field regions where the toroidal field resides turbulent diffusivity will be much less effective, while in the weak field regions occupied by the poloidal component the turbulent diffusivity will be higher (closer to values predicted by mixing length theory). This is a relatively simple, although effective technique for numerically capturing this physics. The model also has a built-in buoyancy algorithm that searches for a toroidal field exceeding a certain threshold (10⁵G by default) at the base of the SCZ (at r = 0.71  R_⊙) at certain time intervals, removes half of it, and deposits this field at near-surface layers where the BL poloidal source is located. This is done to emulate the eruption of magnetic flux tubes. For most of our simulations, Eqs. (1) and (2) are solved in a 256 × 256 grid between R_b = 0.55  R_⊙ < r < R_⊙ and 0 < θ < π under appropriate boundary conditions.

This code also uses a differential rotation profile, Ω(r,θ), that is an analytical fit to the solar internal rotation data provided by Helioseismology. The meridional circulation, v, has a one cell per hemisphere configuration, with latitudinal flows drifting towards the poles in the surface layers and towards the equator just below the tachocline. The amplitude of the surface component of this flow at mid latitudes is defined by v₀ = −29 m s^-1 (in the standard version of the code). Detailed description of mathematical parametrization of these quantities including their justifications, for this well tested dynamo code, are available elsewhere (Nandy & Choudhuri 2002; Chatterjee et al. 2004; Yeates et al. 2008). The standard version of this code implements a surface αeffect that mimics the BL mechanism and is parameterized through (3)where r₁ = 0.95  R_⊙, r₂ = R_⊙, and d₁ = d₂ = 0.025  R_⊙ are scaling factors and α₀ = 25 m s^-1 is the amplitude of the α effect.

As a proxy for the solar cycle, in this study we use the amplitude of the toroidal field component given by (or | B_φ |) just above tachocline depth and at usual active latitudes (~14°). Figure 2 shows a reference solution depicting the evolution of the toroidal field in both hemispheres (red and black) just above the tachocline depth at r = 0.706  R_⊙ and the corresponding synoptic representation (similar to a butterfly diagram) where the positive (negative) erupted field is identified by solid (dashed) contour lines. This reference solution was obtained using α₀ = 27 m s^-1, and the threshold for the buoyancy algorithm B_c = 8 × 10⁴ G. All other parameters are set to the default values as in the public version of the Surya dynamo code.

3. Modelling grand minima episodes with stochastic fluctuations in the α effect

A physical mechanism that has been proposed as responsible for grand minima is the stochastic nature of the poloidal source, traditionally the α effect (Choudhuri 1992; Charbonneau & Dikpati 2000; Proctor 2007; Brandenburg & Spiegel 2008; Moss et al. 2008). As noted before, as flux tubes rise through the CZ they are acted upon by the Coriolis force that imparts to them the tilt angle observed in bipolar sunspot pairs. The values for the average tilt angle presents a uniform distribution centered around 6° (Howard 1991; Dasi-Espuig et al. 2010). The idea behind this tilt dispersion is that during the buoyant rise, the turbulent buffeting of the tubes adds a random component to the tilt angle, contributing to the scatter in the distribution. In turn, this scattering in the tilt angles of active regions will contribute to variations in the efficiency of the surface BL mechanism. This phenomena (variations in the BL poloidal source amplitude) can, in principle, be modelled by introducing stochastic fluctuations in the amplitude of the BL α effect.

With this in mind, we now redefine this coefficient by splitting it into a constant part and a fluctuating part so that α = α₀ + α′σ(t,τ), where σ is a function (with random values between − 1 and 1) that depends on the time t, on the correlation time for fluctuations τ, and the amplitude of the fluctuations α′. Flux tube simulations suggest that their rise time through the CZ is of the order of a few months (e.g. Caligari et al. 1995) and we know that surface flows take of the order of months to redistribute the sunspot flux. Therefore, we choose τ = 6 months. As a first test, we implement fluctuations at 100% level in α independently in both hemispheres. Since the fluctuation levels can be difficult to estimate, the amplitudes we use are motivated by previous works (Charbonneau & Dikpati 2000; Brandenburg et al. 2008), and on the eddy velocity distributions present in some highly turbulent global 3D MHD simulations of solar convection (e.g. Racine et al. 2011; Passos et al. 2012). In addition, significant fluctuation in the BL source term is expected because of the large variation in the tilt angle distribution of bipolar solar active regions (Dasi-Espuig et al. 2010).

An important result of this experiment is that fluctuations on a timescale much smaller (few months) than the solar cycle itself (approximately 11 years) can induce variability on the timescale of decades. Although the main cyclic activity persists, this shows how susceptible the system is to forcing factors.

As an example of the obtained behaviour we present the result of one of the simulations in Fig. 3. In this figure we can identify periods (e.g. between 110 and 140 or between 340 and 390) where two or three cycles are suppressed in the synoptic representation. These periods would correspond to grand minima in which the dynamo is still operating, but in a regime where the produced fields are not strong enough to produce surface eruptions, i.e. the field amplitude stays below the buoyancy threshold (blue dashed line in the graphic).

3.1. Two poloidal sources at play

Although the previous simulation set-up with fluctuations in a dynamo system with a BL poloidal source alone represents a scenario for entry and exit from grand minima-like conditions, as has been demonstrated in earlier simulations (Charbonneau & Dikpati 2000; Karak 2010; Choudhuri & Karak 2012), the physics of the recovery from a grand minimum in these model set-ups is questionable (Hazra et al. 2014). The problem lies in the parametrization of the BL mechanism itself. Since the main idea behind this mechanism¹ is the decay of active regions, one should ensure that this α_BL is only acting on strong fields that have buoyantly erupted and not on weaker fields that might reach the surface layers by advection or diffusion. This implies that a lower threshold in the α_BL must be introduced to prevent this (Nandy 2002b; Charbonneau et al. 2005; Hazra et al. 2014). Therefore we introduce a new parametrization in the α_BL that depends on the field intensity, (4)where r₁, r₂, d₁, and d₂ have the same values as defined for Eq. (3); α_0BL controls the amplitude of the effect; a₁ = 0.393 is a normalization constant; B_1lo = 10³ G corresponds to the lower threshold; d₃ = 10² G; B_1up = 10⁵ G is an upper threshold; and d₄ = 10⁶ G. This formulation ensures that only the fields that reach the surface layers with magnitudes between B_1lo and B_1up contribute to the BL mechanism.

To determine what would be an appropriate value for the lower threshold of α_BL in our model, B_1lo, we conducted the following experiment. We performed simulations where we let the model evolve until it reached a maximum concentration of toroidal field in the base of the convection zone (akin to a cycle maximum). In the following step we use this final state as the initial condition of a new simulation and, to ensure that there is no toroidal field in the main body of the convection zone, we added a condition to wipe out fields above r = 0.725  R_⊙. In other words, we initiate a new simulation with a high concentration of toroidal field just below the tachocline, but none above. We switched off the buoyancy algorithm and let the dynamo solution evolve just in the presence of advective flows and diffusion. The time that the peak associated with this strong field (~10⁵ G at r = 0.725  R_⊙) takes to reach near surface layers at active latitudes, more specifically at a latitude θ = 15° N and r = 0.97  R_⊙ around the radius where α_BL peaks, is approximately 5 years. The field reaching this radius had an amplitude around 750 G. This implies that this field is being transported primarily by the meridional flow towards the surface. The diffusion timescale for the toroidal field is given by τ_t = ℓ²/η_t, and direct substitution (ℓ ≃ 0.245  R_⊙ and η_t = 4 × 10¹⁰ cm² s^-1) yields that τ_t is around 230 years. Even if we consider an average value of η one order of magnitude higher, we still obtain a diffusion time several times higher than the advection time found.

This numerical experiment and associated physical insights point out that even if buoyant eruption is absent, for example during solar grand minima episodes, in such model set-ups, the toroidal field is still dredged up from beneath the base of the convection zone to the surface by meridional circulation. Even though this toroidal field takes 5 years to move up and is weak, it ends up at the near-surface region where the spatially distributed, α-coefficient representing the BL poloidal source is located. Not being able to distinguish whether this toroidal field represents a contribution from sunspot eruptions or from dredging up of the field due to meridional circulation, the BL α source ends up producing poloidal fields in a manner that is un-physical and not in keeping with established ideas of the dynamics of magnetic flux tubes. It is well known, that if it takes many years for a toroidal flux tube to rise up, it will be completely shredded by turbulence and will not have the required tilt angle distribution to contribute to the BL mechanism. Thus, to circumvent this un-physical behaviour in the model set-up, a lower threshold of B_1lo = 10³ G is necessary in such kinematic dynamo models.

The effect of this lower threshold when compared to the reference solution is to lower the amplitude of the solution’s magnetic field and consequently decrease the number of eruptions and the latitudes at which they appear. The value of α_BL = 27 now corresponds to a value just above critical with eruptions first appearing at 25 degrees in latitude. This situation makes the dynamo solution especially susceptible to amplitude fluctuations in this coefficient. In Fig. 4 we show a dynamo solution which decays after half a century with just 10% fluctuations in this redefined BL α-coefficient with the lower threshold. The decay time of the solution scales inversely to the fluctuation levels. For the current set-up, fluctuation levels between 50% and 150% induce the decay of the solution after a few decades, but for lower fluctuation level the solution takes more time to decay. The observed decoupling between northern and southern hemispheres that appear in these solutions are also plausibly connected to the fluctuations levels.

Therefore, while a BL mechanism with both upper and a lower operating thresholds allows for the possibility of entering a grand minima period, it does not seem to permit a self-consistent recovery from grand minima episodes; this conclusion is also supported by simulations with the independent, time-delay dynamo model presented in Paper I. Since the Sun has always emerged from such quiescent periods, one needs to come up with a solution that models this behaviour. The solution we present here is inspired by classical αΩ models and its implementation is based on ideas discussed in Paper I in the context of a reduced dynamo model; i.e. we add a weak mean-field α effect working in conjunction with the surface α_BL. This α_MF operates in the bulk of the convection zone and operates on weak flux tubes that are below a certain threshold and which do not contribute to the formation of sunspots (thus this α effect could be in operation even during grand minima phases).

We parameterize this mean-field α effect as (5)where α_0MF is the amplitude of the effect, r₃ = 0.713  R_⊙, r₄ = R_⊙, and B_2up = 10⁴ G is an upper threshold. The radial profile and quenching profiles for the two poloidal sources (parameterized through two α effects as described above) are presented in Fig. 5. The combined final α effect acting in our system is then defined as α = α_BL + α_MF, the first acting on strong fields that erupt buoyantly to the surface and the second acting on weaker fields that are either advected by meridional flow or diffuse out into the SCZ from the tachocline.

In the absence of any fluctuations, the solution obtained using this combined α with α_BL = 27, B_1lo = 10³ and α_MF = 0.4 is similar to the reference solution (Fig. 2) but with a slight decrease (~1 year) in the period of the cycle and with eruptions starting from somewhat higher latitudes. They are qualitatively similar in other aspects.

3.2. Stochastic fluctuations in the poloidal sources

To test if the combination of these two poloidal sources can produce grand minima-like episodes (under the influence of stochastic fluctuations) and recover self-consistently from these episodes, we apply different levels of fluctuations to both α effects, separately and in conjunction to perform further simulations.

We start by repeating the simulation with 100% fluctuations in the α_BL = 27, but now with the additional presence of α_MF = 0.4 (see Fig. 6). We observe that there are periods when the field falls below the buoyancy threshold, but after some time it regains strength and begins producing sunspot eruptions again. The underlying dynamics probably involves fluctuations in α_BL triggering grand minima periods, and the presence of α_MF facilitating recovery. Hints of these dynamics between two interacting dynamo α effects were already present in Ossendrijver (2000), but where instead of the near-surface BL α effect, a buoyancy instability induced α effect in the overshoot layer was considered.

In these simulations we continue to observe the decoupling between the northern and southern hemispheres, especially in simulations with a high level of fluctuations in the BL mechanism. Occasionally, small changes in the parity occur as well (deviations from dipolar parity). In some extreme cases, both hemispheres sustain dynamo action but in an almost independent manner. During the deepest phase of a grand minimum the hemispheres loose synchronization, but interestingly it is quickly regained after the system emerges from this quiescent period. There is a small, variable phase shift between the two hemispheres that was already present in Fig. 3. This type of weak hemispheric decoupling during normal activity is plausibly due to the presence of fluctuations and is also observed in the solar butterfly diagram. Different levels of fluctuations in the α_BL between 25% and 200% were tested. We observe that the decoupling between hemispheres and the number of grand minima increases with the fluctuation level. Independent simulations in a different context also point out that hemispheric decoupling may be a phenomena that is naturally associated with, and symptomatic of grand minima episodes (Olemskoy & Kitchatinov 2013). In Fig. 6 it can also be seen that, depending on the circumstances, it is possible to have a quiescent minimum-like phase hemisphere, but a failed minimum in the other hemisphere.

In the next batch of simulations we apply fluctuations just to the mean-field poloidal source α_MF using a correlation time of 1 year. The simulations with α_BL = 27 and α_MF = 0.4 + 100% fluctuations returned very few (and short, two cycles at the most) minima episodes. Increasing fluctuations to 200% increases the number of minima, but not their length. To test the importance of α_MF in the dynamics, we decrease the strength of the BL source to α_BL = 21 (barely supercritical). With this lower value, values of 100% to 200% fluctuations in the mean-field produce longer minima episodes. Both hemispheres seem to evolve in a similar way although some variability between north and south is observed. Even for longer simulated periods, this coupled behaviour is maintained as is the parity of the solution. This implies that the frequency of occurrence, and duration of grand minima episodes depend on the relative amplitudes of the two poloidal sources and the level of fluctuations in them. Moreover, it appears that the mean-field poloidal source – distributed across the convection zone – plays a more important role in the maintenance of hemispheric coupling and dipolar parity than the BL poloidal source.

To complete the set of possible numerical experimentation, different levels of fluctuations can be applied to both the poloidal sources. The overall set of obtained solutions, not surprisingly, shows a wide range of behaviours in terms of frequency and duration of grand minima, hemispheric decoupling, and parity change. A representative solution from numerical simulations with intermittency in both poloidal sources is presented in Fig. 8. Here we note a strong hemispheric asymmetry in the activity just prior to entering the grand minimum phase, continued hemispheric decoupling throughout the minimum phase, and regaining of hemispheric coupling following recovery from the minimum. The period during the grand minimum is characterized by an extended period of low magnetic activity; however, there are occasional sunspot eruptions during this period and the solar polar field (dipolar component) continues its regular cycle of reversal albeit with weak field amplitudes. These features are in broad qualitative agreement with observed features of the solar Maunder minimum (Beer et al. 1998; Miyahara et al. 2006).

4. Concluding summary

Solar dynamo models can, in general, explain the main features of the solar cycle. Explanations of solar cycle fluctuations, such as grand minima and maxima episodes, should also be encompassed in these models. Radio-isotope records have shown that our Sun regularly goes through quiescent periods and higher than normal activity phases. To study this activity modulation two different approaches are usually used. In the first approach, one utilizes the role of dynamic, non-linear feedback mechanisms in the solar dynamo (Tobias 1997; Brooke et al. 2002; Bushby 2006) including time delay dynamics (Charbonneau et al. 2005; Wilmot-Smith et al. 2006; Jouve et al. 2010). In the second approach, one utilizes stochastic forcing of the dynamo system motivated from the turbulent nature of the SCZ. In this work we follow the second approach, where we use a stochastically forced kinematic solar dynamo model based on the flux transport paradigm, to focus on solar cycle dynamics related to grand minima episodes.

We also point out that in this model, turbulent diffusion, and not meridional circulation, is the dominant process for transporting the poloidal field in the solar interior where the toroidal field is produced and stored. However, meridional circulation does play an important role in the near-surface poloidal field dynamics and latitudinal transport of toroidal fields deep in the interior. In this context it should be noted that another mechanism, downward turbulent pumping of magnetic flux (Tobias 2001), is also expected to play an important role in solar cycle dynamics (Guerrero & de Gouveia Dal Pino 2008; Nandy & Karak 2012). Moreover, turbulent pumping in a rotating medium may also contribute to the reprocessing of magnetic flux and thus be an effective mechanism for the continuation of weak cycles; we expect that the current model set-up with a mean-field α effective on weak fields would encompass the consequences of such a mechanism implicitly, even in the absence of an explicit treatment of turbulent pumping.

Since the convection zone is highly turbulent, stochastic fluctuations arise naturally in the dynamo source terms. Motivated by the results obtained in a low-order dynamo system in Paper I, showing the importance of a dual source formalism in the context of self-consistent entry and exit from grand minima, we started by adding a lower quenching threshold to the parametrization of the BL poloidal source. Adding this low-threshold to the Babcock-Leighton mechanism ensures that this source only acts on strong fields that erupt to the surface and not on weak fields that leak into the convection zone by advection and diffusion. These weak fields do not produce sunspots, hence they do not participate in the BL mechanism. We believe that this presents a more realistic approach to modelling the BL poloidal source. Although stochastic fluctuations (in α_BL) in models without the lower threshold present grand minima like episodes, we argue that these are not real and are un-physical. When the threshold is included, the model returns solutions with a weaker amplitude of the magnetic field that are very susceptible to stochastic fluctuations and models with this physically justified BL source parametrization fail to recover from a grand minima as pointed out in Paper I.

An additional weak, classical mean-field effect, α_MF, operating in the bulk of the convection zone is introduced in the model. This source acts on toroidal fields that are not strong enough to produce sunspots and are advected by the meridional flow (and diffuses) into the convection zone. Our kinematic dynamo model with the dual poloidal sources presents typical dynamo solutions with emerging patterns similar to the classical sunspot butterfly diagram. This model presents more robust solutions when stochastic fluctuations are introduced to both α effects (separately or simultaneously). For several combinations in the parameter space, we observe behaviours analogous to grand minima, where the eruption of toroidal fields at the surface is suppressed. However, with this dual poloidal source, physically consistent recovery from a grand minimum is possible as demonstrated in this spatially extended dynamo model. Thus, the results obtained herein corroborate and extend those already presented in Paper I (using a completely different modelling technique), which justifies the usage of appropriate low-order models such as the time delay dynamo model used in Paper I to explore long-term phenomena which are computationally expensive for spatially extended dynamo models.

Two possible scenarios for entering and exiting a grand minimum were identified. The first is the case where small fluctuations in the surface α_BL are combined with fluctuations in α_MF in a model with a near critical solution. Since the field is near the limit necessary for the production of sunspots, fluctuations in α_MF can take it below the buoyancy threshold. These grand minima tend to be small with a duration of just a few cycles. The other case is in models with solutions above critical when both source terms are subjected to considerable levels of fluctuation (above 75%). In this case, fluctuations in α_BL tend to make the toroidal field fall abruptly below the threshold necessary for the production of sunspot eruptions; this renders the BL source completely ineffective with consequent disruption of dynamo activity. The combination of different fluctuation levels also controls the onset of the grand minima (abrupt or gradual). In several cases we observe that one of the hemispheres emerges from a minimum a few cycles before the other. This was also apparently the case during the exit from the Maunder minimum (Nagovitsyn et al. 2010), so this model behaviour is entirely plausible.

We find that fluctuations in α_BL tend to decouple the magnetic activity in the two hemispheres as observed in these spatially extended simulations because in our model we allow for different randomly generated fluctuations in the northern and southern hemispheres, which is expected in reality. Nevertheless, for fluctuations of α_BL of the order of 25% the asymmetry between hemispheric magnetic activity is small, but with increasing fluctuation levels this asymmetry increases. We find that the mean-field poloidal source plays a more important role in keeping the two hemispheres coupled, while fluctuations in the BL source tends to introduce parity shifts and decouple the two hemispheres. It is possible, as demonstrated, that under certain conditions, one hemisphere may undergo minimum like conditions with no eruptions while there is a failed minimum in the other hemisphere where sunspot eruptions continue. Independent simulations point out that hemispheric asymmetry may be symptomatic of grand minima episodes; Olemskoy & Kitchatinov (2013) and Charbonneau et al. (2005) provide a possible physical explanation of the role of hemispheric coupling vis-a-vis solar grand minima episodes.

We note that for a large number of combinations of poloidal source amplitudes and fluctuations, we see regular polar field reversals (in radial field evolution) even when the surface BL mechanism remains ineffective over a long period during a grand minimum (because of the lack of bipolar sunspot eruptions). This behaviour in our model is due to the presence of the mean-field poloidal source which maintains a low activity cycles. This result is supported by solar activity reconstructions based on cosmogenic isotopes which indicate regular polar field reversals occurred during the Maunder minimum (Miyahara et al. 2004).

One of the important aspects of extreme solar activity fluctuations is the observed statistics (i.e. the frequency distribution) of solar grand minima and maxima episodes (Usoskin et al. 2007). These constraints are on the duration and waiting periods for minima and maxima episodes. While we have not addressed this aspect in our first comprehensive study with a dynamo driven by dual poloidal field sources, we point out that an independent study does focus on a comparative analysis of the observed statistics of grand minima with a kinematic dynamo model set-up in a similar way (Karak & Choudhuri 2013). Another important aspect that we do not directly focus on here is the phase that immediately precedes a solar grand minimum and the defining characteristics of this phase. Since such studies require detailed investigations, and significant computational investment in simulations, we defer this to the future.

---

Acknowledgments

We are grateful to the Ministry of Human Resource Development, Council for Scientific and Industrial Research, 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. Passos acknowledges support from the Fundação para a Ciência e Tecnologia grant SFRH/BPD/68409/2010, CENTRA-IST, the University of the Algarve for providing office space and he is also grateful to IISER Kolkata for hosting him during the final phases of the research work presented here. We also acknowledge the referee for his critical comments which helped us improve an earlier version of our manuscript.

References

All Figures