HYSTERESIS BETWEEN DISTINCT MODES OF TURBULENT DYNAMOS

The solar magnetic activity cycle is not a strictly periodic phenomenon. Its duration and strength vary from cycle to cycle. An impressive example of this aperiodicity is the famous Maunder minimum when sunspots were extremely scarce over about 70 yr (Hoyt & Schatten 1996). Prolonged events of low magnetic activity like the Maunder minimum are typical characteristics of the Sun. Radiocarbon data reveal solar activity variations for the past ∼11,000 yr, with 27 grand minima covering about 17% of the time (Usoskin et al. 2007; Usoskin 2013).

Based on the extensive literature on nonlinear dynamos displaying long-term variability, we can classify two broad theories of grand minima: amplitude modulation through nonlinearity (Spiegel 1977, pp. 267–283; Tavakol 1978; Ruzmaikin 1981) and externally imposed noise (Choudhuri 1992), as has been extensively reviewed by Charbonneau (2010). Amplitude modulation is found to exist in many nonlinear dynamo models. This can result from the coupling between various dynamo modes with close frequencies (e.g., Brandenburg et al. 1989a, 1989b; Sokoloff & Nesme-Ribes 1994; Beer et al. 1998; Brooke et al. 1998) and/or from the interaction between magnetic field and differential rotation (Kitchatinov et al. 1994; Küker et al. 1999). However, the latter is less likely to apply to the Sun on the grounds that the variation in observed differential rotation is weak. Chaotic behavior of nonlinear dynamo models was also identified to be a cause of amplitude modulation in low-order models (e.g., Weiss et al. 1984). Originally this appeared to be a feature of highly truncated models, but it was later also found in two-dimensional models (Covas et al. 1998). On the other hand, since turbulence is the driver of dynamo action in stars, grand minima through the resulting noise could be possible. Indeed, dynamo coefficients such as the α effect and the Babcock–Leighton type α effect (through variations in the tilt angle of bipolar active regions; Dasi-Espuig et al. 2010) are known to fluctuate (Hoyng 1988; Brandenburg et al. 2008). Therefore, fluctuations in the dynamo parameters (e.g., Choudhuri 1992; Moss et al. 1992; Hoyng 1993; Ossendrijver et al. 1996; Charbonneau et al. 2004; Moss et al. 2008; Choudhuri & Karak 2009; Passos et al. 2014) and even the meridional circulation (Karak 2010) are naturally invoked to explain the origin of grand minima. Turbulence also introduces "magnetic noise" that directly affects the mean electromotive force (Brandenburg & Spiegel 2008). The fluctuations cause irregular changes in dynamo cycle amplitudes with occasional wandering into states of low cycle strength. Dynamo models with fluctuating parameters generally reproduce the grand minima statistics (e.g., Choudhuri & Karak 2012; Karak & Choudhuri 2013; Olemskoy & Kitchatinov 2013). Recent analysis of radiocarbon data by Usoskin et al. (2014), however, showed that grand minima do not constitute a low-activity tail of the distribution common for all activity cycles, but represent a separate activity mode that cannot be interpreted as a fluctuation of the "regular" mode. They concluded that solar dynamo regimes in grand minima and in regular cycles are distinct.

The difference in dynamo operation between grand minima and regular activity modes can be interpreted as a consequence of a hysteresis phenomenon found in nonlinear mean-field dynamo models (Kitchatinov & Olemskoy 2010). In the majority of such models of the solar dynamo, suppression of poloidal field generation by the magnetic field (α-quenching) is invoked. This nonlinearity serves well for stabilizing magnetic field growth. The dynamo amplification of the magnetic field takes place when the dynamo number

$$\begin{eqnarray}&&{\mathcal{D}}=\frac{\alpha {\Delta\Omega}R_{\odot }^{3}}{\eta _{{\rm T}}^{2}}\end{eqnarray} \tag{ 1 }$$

exceeds a critical value ${{{\mathcal{D}}}_{{\rm c}}}$, where ${{\eta }_{{\rm T}}}=\eta +{{\eta }_{{\rm t}}}$ is the total (microphysical plus turbulent) magnetic diffusivity, ${\Delta\Omega}$ is the angular velocity variation in the Sun, and α is the measure of the α effect; see Krause & Rädler (1980). A decrease in α with increasing magnetic field strength reduces the effective dynamo number to saturate the field growth.

However, it is not only α but also the eddy diffusivity ${{\eta }_{{\rm t}}}$ that is magnetically quenched. Using predictions of the quasi-linear theory for α- and ${{\eta }_{{\rm t}}}$-quenching results in a non-monotonous dependence of the effective (magnetically modified) dynamo number on the magnetic field (see Figure 1 in Kitchatinov & Olemskoy 2010). This number initially increases with the magnetic field but the dependence changes to a decrease for stronger fields. If the dynamo number (1) is increased from a subcritical value, the saturated field amplitude jumps to a finite value of the order of the equipartition field just after ${\mathcal{D}}$ exceeds ${{{\mathcal{D}}}_{{\rm c}}}$ and then varies smoothly with increasing ${\mathcal{D}}$ (Rüdiger et al. 1994). If the dynamo number is then decreased, dynamo-generated finite fields survive for ${\mathcal{D}}\lt {{{\mathcal{D}}}_{{\rm c}}}$, but for sufficiently small values of ${\mathcal{D}}$, the field eventually falls to zero. The saturated field amplitude, therefore, depends on the pre-history of the ${\mathcal{D}}$ variation (the dynamo hysteresis). In a finite range of ${\mathcal{D}}$-values, there are two stable solutions with considerably different characteristic strengths of the magnetic field. Fluctuations in dynamo parameters provoke irregular transitions between these two solutions (Kitchatinov & Olemskoy 2010). The dynamo hysteresis can thus explain the distinction between grand minima and regular activity modes found by Usoskin et al. (2014).

Zoom In Zoom Out Reset image size

Kitchatinov & Olemskoy (2010) used a mean-field dynamo model that cannot be free from arbitrary prescriptions. Apart from magnetic quenching of α and ${{\eta }_{{\rm t}}}$, there are other nonlinearities that all are implicitly present in direct numerical simulations. This paper probes the dynamo hysteresis with such simulations. Using a shearing box setup, we perform simulations of helically forced turbulence, which produce oscillating (solar-type) dynamos. In principle, this can also be studied in realistic global rotating magneto-convection simulations in spherical geometry (e.g., Racine et al. 2011; Nelson et al. 2013; Karak et al. 2015), but those simulations are computationally more demanding and would benefit from guidance through simpler turbulence simulations. By varying the amount of relative kinetic helicity, the hysteresis-type dependence of the oscillation amplitude on the pre-history of helicity variations is clearly seen. Similar behavior is also found when magnetic diffusion is varied. The simulations generally confirm the presence of two distinct regimes of large-scale dynamos in the vicinity of the dynamo threshold.

We note that turbulent large-scale dynamos near onset have already been studied by Rempel et al. (2009), who used ABC flow forcing. They found intermittent large-scale fields right after dynamo onset, but in their case no cyclic dynamos were possible nor did they find evidence of two distinct states.

In our model, we assume the fluid to be isothermal and compressible. It obeys the equation of state $p=c_{{\rm s}}^{2}\rho$, with constant sound speed ${{c}_{{\rm s}}}$. Hence we solve the following equations:

$$\begin{eqnarray}\begin{array}{rcl} \frac{D{\boldsymbol{U}} }{Dt} & = & -S{{U}_{x}}\hat{y}-c_{{\rm s}}^{2}\nabla {\rm ln} \rho \\ {} & {} & +{{\rho }^{-1}}[{\boldsymbol{J}} \times {\boldsymbol{B}} +\nabla \cdot (2\rho \nu {\boldsymbol{S}} )]+{\boldsymbol{f}} , \\ \end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\frac{D{\rm ln} \rho }{Dt}=-\nabla \cdot {\boldsymbol{U}} ,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{A}} }{\partial t}+{{{\boldsymbol{\bar{U}}} }^{(S)}}\cdot \nabla {\boldsymbol{A}} =-S{{A}_{y}}\hat{x}+{\boldsymbol{U}} \times {\boldsymbol{B}} +\eta {{\nabla }^{2}}{\boldsymbol{A}} .\end{eqnarray} \tag{ 4 }$$

Here $D/Dt=\partial /\partial t+({\boldsymbol{U}} +{{{\boldsymbol{\bar{U}}} }^{(S)}})\cdot \nabla$ is the advective time derivative, ${{{\boldsymbol{\bar{U}}} }^{(S)}}=(0,Sx,0)$ with $S={\rm constant}$ is the imposed uniform large-scale shear flow, ν is the constant kinematic viscosity, ${\boldsymbol{A}}$ is the magnetic vector potential, ${\boldsymbol{B}} =\nabla \times {\boldsymbol{A}}$ is the magnetic field, ${\boldsymbol{J}} =\mu _{0}^{-1}\nabla \times {\boldsymbol{B}}$ is the current density, η is the constant microscopic diffusivity, and ${\boldsymbol{f}}$ is a forcing function to be specified below. The traceless rate of strain tensor ${\boldsymbol{S}}$ is given by ${{S}_{ij}}=\frac{1}{2}({{U}_{i,j}}+{{U}_{j,i}})-\frac{1}{3}{{\delta }_{ij}}\nabla \cdot {\boldsymbol{U}} ,$ where the commas denote partial differentiation with respect to the coordinate (j or i). The contribution of ${{{\boldsymbol{\bar{U}}} }^{(S)}}$ to ${\boldsymbol{S}}$ is omitted, because it would only introduce a small contribution.

Turbulence is sustained by supplying energy to the system through a forcing function ${\boldsymbol{f}} ={\boldsymbol{f}} ({\boldsymbol{x}} ,t)$, which is helical and random in time (δ-correlated). It is defined as

$$\begin{eqnarray}&&{\boldsymbol{f}} ({\boldsymbol{x}} ,t)=\operatorname{Re}\{N\;{{{\boldsymbol{f}} }_{{\boldsymbol{k}} (t)}}{\rm exp} [i{\boldsymbol{k}} (t)\cdot {\boldsymbol{x}} +i\phi (t)]\},\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{x}}$ is the position vector. At each timestep the wavevector ${\boldsymbol{k}} (t)$ randomly takes any value from many possible wavevectors in a certain range around a given forcing wavenumber ${{k}_{{\rm f}}}$. The phase $-\pi \lt \phi (t)\leqslant \pi$ also changes randomly at every timestep. On dimensional grounds, we choose $N={{f}_{0}}{{c}_{{\rm s}}}{{\left( |{\boldsymbol{k}} |{{c}_{{\rm s}}}/\delta t \right)}^{1/2}}$, where f₀ is a non-dimensional forcing amplitude. The transverse helical waves are produced via Fourier amplitudes (Haugen et al. 2004)

$$\begin{eqnarray}&&{{{\boldsymbol{f}} }_{{\boldsymbol{k}} }}={\boldsymbol{R}} \cdot {\boldsymbol{f}} _{}^{{\boldsymbol{k}} ({\rm nohel})}\quad {\rm with}\quad {{R}_{ij}}=\frac{{{\delta }_{ij}}-i\sigma {{\epsilon }_{ijk}}{{{\hat{k}}}_{k}}}{\sqrt{1+{{\sigma }^{2}}}},\end{eqnarray} \tag{ 6 }$$

where σ is a measure of the helicity of the forcing; for positive maximum helicity, $\sigma =1$. The nonhelical forcing function, ${\boldsymbol{f}} _{}^{{\boldsymbol{k}} ({\rm nohel})}=\left( {\boldsymbol{k}} \times {\rm \hat{e}} \right)/\sqrt{{{{\boldsymbol{k}} }^{2}}-{{\left( {\boldsymbol{k}} \cdot {\rm \hat{e}} \right)}^{2}}},$ where ${\rm \hat{e}}$ is an arbitrary unit vector not aligned with ${\boldsymbol{k}}$. Note that $|{{{\boldsymbol{f}} }_{{\boldsymbol{k}} }}{{|}^{2}}=1$ and ${{{\boldsymbol{f}} }_{{\boldsymbol{k}} }}\cdot {{\left( i{\boldsymbol{k}} \times {{{\boldsymbol{f}} }_{{\boldsymbol{k}} }} \right)}^{*}}=2\sigma k/(1+{{\sigma }^{2}})$.

The fluid and magnetic Reynolds numbers and the magnetic Prandtl number are defined as

$$\begin{eqnarray}&&\operatorname{Re}={{u}_{{\rm rms}}}/\nu {{k}_{{\rm f}}},\;{{R}_{{\rm m}}}={{u}_{{\rm rms}}}/\eta {{k}_{{\rm f}}},\;{{P}_{{\rm m}}}=\nu /\eta ,\;\end{eqnarray} \tag{ 7 }$$

where ${{u}_{{\rm rms}}}={{\langle {{{\boldsymbol{u}} }^{2}}\rangle }^{1/2}}$ is the rms value of the velocity in the statistically stationary state with $\langle \cdot \rangle$ denoting the average over the whole domain and ${{k}_{{\rm f}}}$ is the mean forcing wavenumber.

The boundary conditions are shearing–periodic in the x direction and periodic in the y and z directions, with dimensions ${{L}_{x}}={{L}_{y}}={{L}_{z}}=2\pi$. We always choose $S=-0.2$, ${{f}_{0}}=0.01$, and ${{k}_{{\rm f}}}=5{{k}_{1}}$ or $3{{k}_{1}}$, where ${{k}_{1}}=2\pi /{{L}_{x}}=1$ is the smallest possible wavenumber of the box. We use non-dimensional units by setting ${{c}_{{\rm s}}}={{\rho }_{0}}={{\mu }_{0}}=1$, where ${{\rho }_{0}}=\langle \rho \rangle$ is the volume-averaged density, which is constant in time. As initial conditions we take ${\boldsymbol{u}} ={\rm ln} \rho =0$ and a small-scale low amplitude (${{10}^{-4}}$) Gaussian noise for ${\boldsymbol{A}}$. All computations are performed using the Pencil Code.⁵ The grid resolution of all runs presented in this paper is 96 × 96 × 96.

We begin by focusing on the following two sets of simulations. In one, we have simulations for different values of σ, by taking weak fields as initial conditions (described in Section 2). These simulations are performed to identify the onset of the dynamo. In the other set, we take as initial conditions a snapshot of a previous simulation right after the onset of dynamo action. The corresponding dynamo solution is oscillatory. We perform several simulations by successively reducing the helicity by a small value and using the resulting fields of the previous simulation as the initial condition. We continue this procedure until there is only a decaying solution. In this way we identify the regime of dynamo hysteresis as the location where, depending on the initial conditions, both non-decaying oscillatory dynamos and decaying solutions are possible, i.e., the system becomes bistable. Finally, we repeat the whole procedure in a different parameter regime and explore the robustness of the results.

3.1. Onset of Dynamo Action

We perform a set of simulations by increasing the strength of the helicity parameter σ of the turbulent forcing, starting from 0 to 1 (Set I). For this set we fix $\eta =\nu =0.005$ and ${{k}_{{\rm f}}}=5{{k}_{1}}$. Runs A–M in Table 1 show these simulations. Along with other important parameters, we show a rough measure of the dynamo number defined as $D={{C}_{\alpha }}{{C}_{{\Omega}}}$, where ${{C}_{\alpha }}={{\alpha }_{0}}/{{\eta }_{{\rm T}}}{{k}_{1}}$ and ${{C}_{{\Omega}}}=|S|/{{\eta }_{{\rm T}}}k_{1}^{2}$, with ${{\alpha }_{0}}=-\frac{1}{3}\tau \langle {\boldsymbol{ \omega }} \cdot {\boldsymbol{u}} \rangle$, $\tau ={{\left( {{u}_{{\rm rms}}}{{k}_{{\rm f}}} \right)}^{-1}}$, ${{\eta }_{{\rm T}}}=\eta +{{\eta }_{{\rm t}0}},$ and ${{\eta }_{{\rm t}0}}=\frac{1}{3}\tau \langle {{{\boldsymbol{u}} }^{2}}\rangle$. In Figure 1 we show the temporal mean in the statistically stationary state of the large-scale magnetic field over the whole domain, ${{\bar{B}}_{{\rm rms}}}=\langle \langle {{B}_{x}}\rangle _{y}^{2}+\langle {{B}_{y}}\rangle _{y}^{2}+\langle {{B}_{z}}\rangle _{y}^{2}\rangle _{xzt}^{1/2},$ normalized by ${{B}_{{\rm eq}}}={{u}_{{\rm rms}}}$. We see that the large-scale field is zero when σ is below about 0.32, implying there is no dynamo action. For $\sigma =0.32$ (Run D) in Figure 2 we show the spatio-temporal variation of the y-component of the mean magnetic field ${{\bar{B}}_{y}}={{\langle {{B}_{y}}\rangle }_{xy}}$ (which corresponds to the toroidal field in spherical coordinates) and the time series of $\bar{B}_{y}^{2}$ at an arbitrarily chosen mesh point, normalized by $B_{{\rm eq}}^{2}$ (which may be considered as a measure of sunspot number). Here we do not see clear magnetic oscillations. A few cycles started to appear at around $t=2$, but they did not survive. The overall field is also very weak. On increasing σ slightly we observe a dynamo transition at $\sigma =0.322$ (Run E) and the magnetic field becomes strong (${{\bar{B}}_{{\rm rms}}}\gt {{B}_{{\rm eq}}}$). Hence the critical value of σ for dynamo action is ${{\sigma }_{{\rm c}}}\approx 0.322$. The spatio-temporal variation for this case is shown in Figure 3, where we see clear magnetic cycles with dynamo wave propagation along the positive z direction. Together with positive helicity (which results in a negative α) and negative shear, migration in the positive z direction is indeed expected. However, the cycles are not regular; the amplitude varies from cycle to cycle, similar to the observed solar cycle. We recall that in stochastically forced mean-field dynamo models, the cycle irregularity is related to the amount of imposed fluctuations and the corresponding coherence times (see examples in Charbonneau 2010; Karak & Choudhuri 2011). In the present simulations, however, this cycle irregularity is naturally coming because of the finite number of eddies and it is directly related to the scale-separation ratio (${{k}_{{\rm f}}}/{{k}_{1}}$), which has been demonstrated in Brandenburg & Guerrero (2012). For $\sigma \gt {{\sigma }_{{\rm c}}}$ we always observe clear magnetic cycles and the value of ${{\bar{B}}_{{\rm rms}}}$ remains around ${{B}_{{\rm eq}}}$; see Figure 1.

Zoom In Zoom Out Reset image size

3.2. Dynamo Hysteresis

Now we take a snapshot of a simulation at $\sigma ={{\sigma }_{{\rm c}}}=0.322$ (Run E, which is shown in Figure 3) and perform a set of simulations by reducing σ slowly and taking the output of the previous simulation as the initial conditions. Runs E1–E12 in Table 1 represent such cases and the corresponding ${{\bar{B}}_{{\rm rms}}}$ are shown in Figure 1 as red diamonds. We observe oscillatory solutions over a broad parameter range, $0.22\leqslant \sigma \lt 0.322$, in which there are otherwise decaying solutions when started from weak fields. Therefore, in this range, the results depend on the initial conditions, i.e., system becomes bistable. All the simulations are run for a sufficiently long time to ensure that they remain in the same state. We recall that Brummell et al. (2001) studied the linear and nonlinear dynamo properties using time-dependent ABC flows forcing in triply periodic Cartesian geometry. Their simulations are similar to those of Rempel et al. (2009), but for an incompressible fluid. In the nonlinear regime, Brummell et al. (2001) found two distinct classes of behavior depending on the initial hydromagnetic properties of the forced ABC flow, similar to earlier results by Fuchs et al. (1999) in spherical geometry. One produces the stationary solution followed by an initial exponential growth of the magnetic field, whereas the other initially produces a dynamo solution but later turns into a decaying one because the flow itself evolves to a non-dynamo stage through hydrodynamic instability. However, our study of hysteresis is different from Brummell et al. (2001) because we take different initial conditions for velocity as well as magnetic fields and we believe that the magnetic quenching rather than the hydrodynamic instability is the cause of the bistability.

Zoom In Zoom Out Reset image size

Table 1.  Summary of the Runs

Set I (σ Varied)									Set II (Pm Varied)
Run	σ	${{u}_{{\rm rms}}}/{{c}_{{\rm s}}}$	${{R}_{{\rm m}}}$	D	$\tilde{{\bar{B}}}$	${{\tilde{{\bar{B}}}}_{x}}$	${{\tilde{{\bar{B}}}}_{y}}$	Osc	Run	Pm	${{u}_{{\rm rms}}}/{{c}_{{\rm s}}}$	${{R}_{{\rm m}}}$	D	$\tilde{{\bar{B}}}$	${{\tilde{{\bar{B}}}}_{x}}$	${{\tilde{{\bar{B}}}}_{y}}$	Osc
A	0.25	0.32	12.7	3.28	0.00	0.00	0.00	N	N	0.062	0.23	0.6	0.80	0.00	0.00	0.00	N
B	0.27	0.30	11.9	4.10	0.00	0.00	0.00	N	O	0.100	0.34	1.3	0.86	0.00	0.00	0.00	N
C	0.31	0.35	13.6	3.33	0.00	0.00	0.00	N	P	0.167	0.28	1.8	2.49	0.00	0.00	0.00	N
D	0.32	0.32	12.5	4.26	0.00	0.00	0.00	N	Q	0.250	0.31	3.1	3.07	0.00	0.00	0.00	N
E	0.322	0.12	4.8	22.0	3.14	0.21	3.08	Y	R	0.263	0.27	2.8	4.35	0.00	0.00	0.00	N
F	0.325	0.12	4.9	22.1	3.14	0.21	3.09	Y	S	0.278	0.25	2.7	5.57	0.00	0.00	0.00	N
G	0.327	0.12	4.8	22.7	3.17	0.21	3.16	Y	T	0.294	0.16	1.8	6.70	1.68	0.23	1.59	Y
H	0.33	0.12	4.8	23.0	3.19	0.21	3.14	Y	U	0.312	0.18	2.3	5.21	1.51	0.20	1.48	Y
I	0.35	0.21	8.3	7.31	1.09	0.18	0.93	Y	V	0.357	0.19	2.6	5.74	1.66	0.22	1.75	Y
J	0.40	0.22	8.6	7.47	1.00	0.20	0.76	Y	W	0.500	0.19	3.8	7.18	2.03	0.29	2.93	Y
K	0.45	0.22	8.7	7.79	0.97	0.21	0.73	Y	X	0.556	0.23	5.0	6.14	1.12	0.22	0.93	Y
L	0.50	0.23	8.9	7.89	0.98	0.22	0.71	Y	Y	0.625	0.23	5.8	6.51	1.05	0.23	0.81	Y
M	1.00	0.23	9.1	9.14	0.97	0.27	0.64	Y	Z	0.714	0.24	6.7	6.69	1.03	0.23	0.87	Y
E1	0.32	0.13	5.0	20.0	3.01	0.20	2.96	Y	T1	0.278	0.13	1.4	8.52	1.92	0.29	1.88	Y
E2	0.31	0.24	9.6	4.88	1.47	0.10	1.43	Y	T2	0.250	0.13	1.3	7.45	1.69	0.27	1.65	Y
E3	0.30	0.28	11.2	3.30	1.13	0.09	1.06	Y	T3	0.200	0.13	1.0	6.25	1.34	0.26	1.29	Y
E4	0.29	0.25	10.0	4.16	1.39	0.10	1.37	Y	T4	0.179	0.14	1.0	5.29	1.10	0.23	1.04	Y
E5	0.28	0.30	11.9	2.74	1.00	0.08	0.93	Int	T5	0.167	0.15	1.0	4.54	0.92	0.20	0.86	Y
E6	0.27	0.34	13.4	2.12	0.78	0.07	0.72	Int	T6	0.161	0.15	0.9	4.43	0.88	0.20	0.82	Y
E7	0.26	0.36	13.9	1.91	0.72	0.06	0.66	Int	T7	0.156	0.33	2.0	1.60	0.00	0.00	0.00	N
E8	0.25	0.37	14.7	1.77	0.67	0.06	0.63	Int	T8	0.154	0.35	2.1	1.44	0.00	0.00	0.00	N
E9	0.24	0.34	14.6	1.97	0.35	0.04	0.29	Int	T9	0.152	0.32	1.9	1.71	0.00	0.00	0.00	N
E10	0.23	0.33	12.8	2.36	0.52	0.05	0.45	Int	T10	0.147	0.21	1.2	3.25	0.00	0.00	0.00	N
E11	0.22	0.28	11.0	3.17	0.52	0.05	0.43	Int	T11	0.143	0.39	2.2	1.06	0.00	0.00	0.00	N
E12	0.21	0.32	12.5	2.93	0.00	0.00	0.00	N	T12	0.132	0.31	1.6	1.47	0.00	0.00	0.00	N

Note. Runs A–E12 belong to Set I in which σ is varied. Runs A–M are started from weak seed fields. Run E1 is performed from a snapshot of Run E (bold), but at slightly reduced σ. A similar procedure is continued from Runs E1 $\to$ E2 $\to$ E3 ... $\to$ E12. Runs N–T12 belong to Set II in which ${{P}_{{\rm m}}}$ is varied. Runs N–Z are started from weak seed fields. Run T1 is started from a snapshot of Run T (bold), but at decreased ${{P}_{{\rm m}}}$, and a similar procedure is performed for Runs T1 $\to$ T2 $\to$ T3 ... $\to$ T12. $\tilde{{\bar{B}}}={{\bar{B}}_{{\rm rms}}}/{{B}_{{\rm eq}}}$, ${{\tilde{{\bar{B}}}}_{x}}={{\langle \bar{B}_{x}^{2}\rangle }^{1/2}}/{{B}_{{\rm eq}}}$, and $\tilde{{{{\bar{B}}}_{y}}}={{\langle \bar{B}_{y}^{2}\rangle }^{1/2}}/{{B}_{{\rm eq}}}$. The columns "Osc" indicate whether there are oscillations (Y) or not (N) and "Int" denotes intermittent behavior.

Download table as:  ASCIITypeset image

In Figure 4 we show the magnetic oscillations from the last run (Run E11) at $\sigma =0.22$, below which the oscillations die completely. We see that the magnetic cycles persist most of the time in this simulation. The interesting feature is that occasionally some of the cycles disappear or become weaker. By comparing the first two panels of Figure 4 we note that during weaker cycles (for example, at $t\sim 130$ and 154), ${{\bar{B}}_{x}}$ is not reduced as much as ${{\bar{B}}_{y}}$, implying that the α effect dominates over the Ω effect. This kind of intermittent behavior somewhat resembles the grand minima observed in the Sun. We see a similar behavior for many runs in the bistable region, particularly in Runs E5–E11.

Zoom In Zoom Out Reset image size

3.3. Dynamo Hysteresis: Dependence on ${{P}_{{\rm m}}}$

To explore the robustness of the existence of dynamo hysteresis, we repeat the same procedure in different parameter regimes of the simulations. Here we fix σ at 1 (fully helical flow) but vary the magnetic diffusivity η in each simulation; see Runs N–Z of Set II in Table 1. Hence, ${{P}_{{\rm m}}}$ varies, but $\nu =0.005$ is unchanged. This is similar to the experiments of Rempel et al. (2009), who used ABC flow forcing.

The black points in Figure 5 show ${{\bar{B}}_{{\rm rms}}}$ from different simulations started with weak seed fields as the initial conditions (Runs N–Z). We see that when ${{P}_{{\rm m}}}$ just exceeds about 0.29, the dynamo is excited and the magnetic field becomes oscillatory. The critical ${{P}_{{\rm m}}}$ for dynamo action is $P_{{\rm m}}^{{\rm c}}\approx 0.294$.

Zoom In Zoom Out Reset image size

Next, as before, we take an oscillatory dynamo solution (Run T) as the initial condition for the new simulation and decrease ${{P}_{{\rm m}}}$ by a small value progressively in each simulation by taking as initial conditions the last snapshot from the previous simulation. Runs T1–T12 in Table 1 are such examples and the corresponding ${{\bar{B}}_{{\rm rms}}}$ are shown as red diamonds in Figure 5. We see that, up to about ${{P}_{{\rm m}}}=0.16$, we obtain an oscillatory large-scale magnetic field. Figure 6 shows the typical magnetic cycles from a simulation at ${{P}_{{\rm m}}}=0.1613$ (Run T6) below which the oscillation dies. Note that in the bistable stage, unlike in the previous set of simulations, for example in Figure 4 where some of the magnetic cycles disappear occasionally, here we observe cycles all of the time. However, when we repeat the whole procedure at $\sigma =0.5$ instead of 1, we see this kind of intermittent behavior in the bistable regime. To obtain even more confidence in the results we repeated another set of simulations (Set III) at ${{k}_{{\rm f}}}=3$ and $\nu =8\times {{10}^{-3}}$. These simulations are at slightly higher ${{R}_{{\rm m}}}$. Table 2 gives a summary of these runs, and Figure 7 shows the corresponding dynamo hysteresis. We clearly see a similar behavior.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Same as Set II in Table 1 but Simulations are Performed at Forcing Wavenumbers ${{k}_{{\rm f}}}=3$ and $\nu =8\times {{10}^{-3}}$

Set III
Run	${{P}_{{\rm m}}}$	${{u}_{{\rm rms}}}/{{c}_{{\rm s}}}$	${{R}_{{\rm m}}}$	D	$\tilde{{\bar{B}}}$	${{\tilde{{\bar{B}}}}_{x}}$	${{\tilde{{\bar{B}}}}_{y}}$	Osc
A'	0.50	0.28	5.5	2.63	0.00	0.00	0.00	N
B'	1.00	0.38	15.0	1.71	0.00	0.00	0.00	N
C'	1.14	0.38	17.4	1.76	0.00	0.00	0.00	N
D'	1.33	0.32	17.2	2.77	0.00	0.00	0.00	N
E'	1.60	0.14	9.1	14.81	2.23	0.19	2.08	Y
F'	2.00	0.13	10.8	18.66	2.50	0.19	2.39	Y
G'	2.67	0.13	14.0	23.13	2.59	0.18	2.36	Y
H'	4.00	0.13	21.3	25.01	2.85	0.15	2.63	Y
I'	8.00	0.13	42.6	29.19	2.52	0.12	2.25	Y
E'1	1.14	0.13	6.0	13.99	2.20	0.22	2.09	Y
E'2	1.00	0.14	5.5	11.15	1.95	0.21	1.85	Y
E'3	0.89	0.14	4.8	10.74	1.89	0.22	1.80	Y
E'4	0.80	0.15	4.7	8.36	1.72	0.21	1.61	Y
E'5	0.73	0.14	4.0	8.79	1.72	0.22	1.64	Y
E'6	0.67	0.15	3.9	7.44	1.56	0.21	1.48	Y
E'7	0.62	0.14	3.4	7.59	1.51	0.22	1.43	Y
E'8	0.57	0.14	3.1	7.48	1.46	0.22	1.36	Y
E'9	0.53	0.15	3.2	5.64	1.28	0.21	1.23	Y
E'10	0.50	0.36	7.2	1.36	0.00	0.00	0.00	N

Note. Runs E'1–E'10 have been restarted from Run E' (bold).

Download table as:  ASCIITypeset image

To understand why the hysteresis discussed in this paper has not been seen before, we need to assess more carefully the parameter regimes of our solutions. Given that we use (shearing) periodic boundary conditions, there are no magnetic helicity fluxes in or out of the domain, so we can describe the solutions by comparing with the analytic results of Blackman & Brandenburg (2002, hereafter BB02). The dynamical quenching theory used in BB02 was already compared with numerical solutions by Käpylä & Brandenburg (2009). One of the predictions they tested was the saturation level of the mean magnetic field, which they gave in the form

$$\begin{eqnarray}&&\frac{\bar{B}_{{\rm rms}}^{2}}{B_{{\rm eq}}^{2}}\approx \frac{{{\epsilon }_{{\rm f}}}{{k}_{{\rm f}}}}{{{\epsilon }_{{\rm m}}}{{k}_{{\rm m}}}}-\left( 1+\frac{\eta }{{{\eta }_{{\rm t}0}}} \right),\end{eqnarray} \tag{ 8 }$$

where ${{k}_{{\rm f}}}$ and ${{k}_{{\rm m}}}$ are, respectively, the effective wavenumbers of the fluctuating and mean fields, defined via $k_{{\rm f}}^{2}=\langle {\boldsymbol{j}} \cdot {\boldsymbol{b}} \rangle /\langle {\boldsymbol{a}} \cdot {\boldsymbol{b}} \rangle$ and $k_{{\rm m}}^{2}=\langle {\boldsymbol{\bar{J}}} \cdot {\boldsymbol{\bar{B}}} \rangle /\langle {\boldsymbol{\bar{A}}} \cdot {\boldsymbol{\bar{B}}} \rangle$, and ${{\epsilon }_{{\rm f}}}$ and ${{\epsilon }_{{\rm m}}}$ are their fractional helicities, defined via

$$\begin{eqnarray}&&{{\epsilon }_{{\rm f}}}{{k}_{{\rm f}}}={{\mu }_{0}}\langle {\boldsymbol{j}} \cdot {\boldsymbol{b}} \rangle /\langle {{{\boldsymbol{b}} }^{2}}\rangle ,\quad {{\epsilon }_{{\rm m}}}{{k}_{{\rm m}}}={{\mu }_{0}}\langle {\boldsymbol{\bar{J}}} \cdot {\boldsymbol{\bar{B}}} \rangle /\langle {{{\boldsymbol{\bar{B}}} }^{2}}\rangle .\end{eqnarray} \tag{ 9 }$$

Figure 8 shows our data from two sets of simulations (Sets I–II) that produce significant large-scale magnetic fields, i.e., Runs E–E11 from Set I and Runs T–T6 from Set II. We have added labels to some of the data points to identify the runs. We note that, unlike Figures 1 and 5, where we have plotted ${{\bar{B}}_{{\rm rms}}}/{{B}_{{\rm eq}}}$, here we plot $\bar{B}_{{\rm rms}}^{2}/{{B}_{{\rm eq}}}^{2}$, but in a smaller range, which is why the data in Figure 8 show a nearly linear variation. Even in a limited range, apart from a small offset, there is a reasonable agreement between our data and the theory given by Equation (8); see the dotted line in Figure 8. However, Käpylä & Brandenburg (2009) had data in a wider range and found better agreement for higher field strengths. In Figure 8 we observe that, when increasing the helicity parameter σ (from Runs E to M in Set I), both the wavenumber ratio of fluctuating to mean fields, as defined by Equation (9), and the strength of the mean field decrease. A similar trend is followed while decreasing σ (Runs E1–E8), except for the last few runs (Runs E9–E11), which deviate significantly. Qualitatively similar behavior is observed in Set II, when ${{P}_{{\rm m}}}$ is decreased from Runs T1 to T6, although they consistently deviate from the other runs. However for Runs T–Z, the trend is not monotonous.

Zoom In Zoom Out Reset image size

Another prediction of BB02 was that ${{\epsilon }_{{\rm m}}}$ is directly proportional to the ratio of poloidal to toroidal magnetic field amplitudes via

$$\begin{eqnarray}&&{{\epsilon }_{{\rm m}}}={{\left( 2\langle \bar{B}_{x}^{2}\rangle /\langle \bar{B}_{y}^{2}\rangle \right)}^{1/2}}.\end{eqnarray} \tag{ 10 }$$

From Equation (9) we compute ${{\epsilon }_{{\rm m}}}$ by assuming ${{k}_{{\rm m}}}=-{{k}_{1}}$ and show in Figure 9 a scatter plot of ${{\epsilon }_{{\rm m}}}$ versus ${{(2\langle \bar{B}_{x}^{2}\rangle /\langle \bar{B}_{y}^{2}\rangle )}^{1/2}}$. Here we see better agreement with Equation (10), as indicated by the dotted line. In Set I, with increasing σ from Runs E to M, the ratio of the poloidal to toroidal field increases, but the same happens from Runs E1 to E11 with decreasing σ as well, which was unexpected. The same trend is observed as we decrease ${{P}_{{\rm m}}}$ from Runs T1 to T6 in Set II, but in Runs T–Z the variation is not monotonous. Furthermore, we note that the data points corresponding to two different regimes—subcritical and supercritical—lie on different lines (compare red and black points in Figure 9).

Zoom In Zoom Out Reset image size

We note that in Figure 8, and also to some extent in Figure 9, the simulation data are systematically below the analytically expected values. A smaller value of ${{\bar{B}}_{{\rm rms}}}$ could readily be explained as being a combination of several modes, which results in reduced averages. This is reminiscent of the fratricide $\alpha {\Omega}$ dynamos of Hubbard et al. (2011), who found that they can be destroyed by their growing ${{\alpha }^{2}}$ dynamo siblings. An important difference, however, is that they never found the recovery of the $\alpha {\Omega}$ dynamo. This might be related to different orderings of the onsets of ${{\alpha }^{2}}$ and $\alpha {\Omega}$ dynamo action, but this has not been investigated further.

As discussed in the Introduction, dynamo hysteresis predicted by the nonlinear mean-field model of Kitchatinov & Olemskoy (2010) may be relevant to the understanding of distinct modes of solar activity found by Usoskin et al. (2014). Our simulations of turbulent dynamos in shearing boxes with helically forced turbulent flows, which resemble the equivalent $\alpha {\Omega}$ solar dynamo, demonstrate, for the first time, hysteresis behavior. By performing several simulations, either by varying the helicity parameter σ or the magnetic Prandtl number ${{P}_{{\rm m}}}$, we observe two stable states with largely different characteristic field strengths depending on the initial conditions of the simulation. A decaying solution is obtained when the simulation is started with weak random seed fields, but otherwise an oscillatory solution is obtained when the simulation is started from a snapshot of a previous oscillatory dynamo.

We emphasize that our simulations show hysteresis close to the dynamo onset only. This raises the question of whether the Sun may be close to the marginal state for the dynamo. Stellar observations indicate that this might indeed be the case. Magnetic activity is known to be correlated with rotation rate, but leads to angular momentum loss through a magnetically coupled stellar wind (Kraft 1967; Hartmann & Noyes 1987). Spin-down of solar-type stars does not continue above a certain rotation period that depends on the spectral type (Rengarajan 1984). The maximum period probably corresponds to the rotation rate where the global dynamo ceases. The maximum period for G2 dwarfs is only slightly larger than the rotation period of the Sun (see Figure 1 in Rengarajan 1984). The stars showing low magnetic activity similar to the solar grand minima are typically old and slow rotators (Saar & Baliunas 1992).

We note that the grand minima in models based on stochastic fluctuations (e.g., Hoyng 1993; Hoyng et al. 1994; Brandenburg & Spiegel 2008; Moss et al. 2008) always result in random deviations from a single (regular) state. Thanks to the reconstructed solar activity record (Usoskin et al. 2014), however, a different picture emerges in that the grand minima cannot be described in terms of random fluctuations of a single solar-activity mode, but are distinct from the regular mode and produced as a result of sudden transitions from the regular mode to weak-field mode.

Transitions between the distinct dynamo regimes may be caused by small-scale hydromagnetic fluctuations inherent to 3D simulations. The transitions do indeed happen in the simulations with relative kinetic helicity below unity (Figure 4).

The intermittency of distinct dynamo regimes, however, disappeared for maximally helical forcing (Figure 6). In mean-field language, an increase in fractional helicity changes the dynamo from $\alpha {\Omega}$ toward ${{\alpha }^{2}}$ type. There is no hysteresis for the ${{\alpha }^{2}}$ dynamo. In any case, maximal helical forcing would not be realistic for the Sun. The dependence of intermittency on the (fractional) helicity should be explored further in future work.

Besides the demonstration of a hysteresis phenomenon in a turbulent dynamo, we have analyzed the simulation data to compare with the dynamical theory of BB02. We observe that the data from subcritical and supercritical dynamos behave in a qualitatively similar way and that they are in reasonable agreement with the theoretical predictions. We end by remarking that, although our findings of hysteresis between two distinct modes of dynamos is relevant to the recently discovered bi-modal solar activity of Usoskin et al. (2014), our simulations are far from the real Sun. Therefore, future research is necessary to explore similar behaviors in more realistic setups.

We thank an anonymous referee for careful review and valuable comments. L.L.K. is thankful to the Russian Foundation for Basic Research (project 13-02-00277) for the support and A.B. acknowledges support through the Swedish Research Council grants 621-2011-5076 and 2012-5797 and the Research Council of Norway under the FRINATEK grant 231444. The computations have been carried out at the National Supercomputer Centres in Linköping and Umeå, the Center for Parallel Computers at the Royal Institute of Technology in Sweden, and the Nordic High Performance Computing Center in Iceland.