© ESO, 2010

1. Introduction

The solar dynamo models developed so far that agree with solar magnetic field observations usually solve the αΩ mean field dynamo equations. The turbulent α effect first proposed by Parker (1955) is believed to be generated by helical turbulence in the convection zone of the Sun. Because α is generated by quadratic correlations of the small-scale turbulence, we need a closure in order to complete the set of mean field equations, e.g., the first order smoothing approximation (FOSA), and express the mean electromotive force in terms of the mean magnetic fields. This turbulent α encounters a critical problem when the energy of the mean field becomes comparable to the equipartition energy of the turbulence in the convection zone, and therefore it becomes increasingly difficult for the helical turbulence to twist rising blobs of magnetic field. The solar dynamo modellers have traditionally used what is referred to as algebraic alpha quenching to mimic this phenomena. This involves replacing α by ${\mathit{\alpha }}_{\mathrm{0}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}{\overline{\mathit{B}}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{eq}}^{\mathrm{2}}\mathrm{)}$, an expression introduced by Jepps (1975), or by ${\mathit{\alpha }}_{\mathrm{0}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{R}}_{\mathrm{m}}{\overline{\mathit{B}}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{eq}}^{\mathrm{2}}\mathrm{)}$, where α₀ is the unquenched value and R_m is the magnetic Reynolds number, $\overline{\mathit{B}}$ is the mean magnetic field and B_eq is the equipartition magnetic field. The latter expression has been discussed since the early work of Vainshtein & Cattaneo (1992). The R_m in the denominator is included because the small-scale fluctuating magnetic field reaches equipartition long before the mean magnetic field does. This has been supported by several numerical experiments to determine the saturation behaviour of α (e.g. Cattaneo & Hughes 1996; Ossendrijver et al. 2002). Given the large magnetic Reynolds numbers of astronomical objects, this phenomenon is referred to as catastrophic quenching.

After the discovery of the layer of strong radial shear (called the tachocline by Spiegel & Zahn 1992) at the bottom of the solar convection zone, Parker (1993) proposed a new class of solar dynamo models called the interface dynamo. In these models the shear is confined to a narrow overshoot layer immediately beneath the convection zone, which is also the region of α effect. The dynamo wave propagates in a direction given by the Parker-Yoshimura rule at the interface between the two layers defined by a steep gradient in the turbulent diffusivity. The toroidal field produced by the stretching effect of the shear is much stronger than the poloidal field and remains confined in the overshoot layer, away from the region where the α effect operates. Note that the interface dynamo model may have serious problems when solar-like rotation with positive latitudinal shear is included (Markiel & Thomas 1999). Similarly, in the Babcock-Leighton class of flux transport models (Choudhuri et al. 1995; Durney 1995) the toroidal and the poloidal fields are produced in two different layers. Unlike in the interface dynamo models, the coupling between the two layers is mediated both by diffusion and the conveyor belt mechanism of the meridional circulation.

It has been proposed that in interface and Babcock-Leighton type dynamos, the α effect is not catastrophically quenched at high R_m because the strength of the toroidal field is very weak in the region of finite turbulent α (e.g. Tobias 1996; Charbonneau 2005). However, according to our knowledge, not much has been done to study the variation of the strength of the saturation magnetic field with the magnetic Reynolds number for these classes of αΩ dynamos. Zhang et al. (2006) made an attempt to reproduce the surface observations of current helicity in the Sun with a 2D mean field dynamo model in spherical coordinates coupled with the dynamical quenching equation. In a separate paper (Chatterjee et al. 2010) we demonstrated that interface dynamo models are also subject to catastrophic quenching.

Magnetic helicity has been identified as a key player in any dynamo process. Because the net magnetic helicity is conserved, any process can only create large- and small-scale magnetic helicities of opposite signs. This was shown in direct numerical simulations of a twisted magnetic flux tube by Blackman & Brandenburg (2003), in good agreement with observations of sigmoid loops in the solar corona. However, the small-scale magnetic helicity backreacts on the helical turbulence and quenches the dynamo (Blackman & Field 2000; Kleeorin et al. 2000).

It has now been shown that this mechanism reduces the saturation strength of the magnetic field (B_sat) with increasing magnetic Reynolds number (R_m). Nevertheless this constraint may be lifted if the system is able to rid itself of small-scale helicity through at least one of several ways such as open boundaries, advective, diffusive and shear-driven fluxes (Shukurov et al. 2006; Zhang et al. 2006; Sur et al. 2007; Käpylä et al. 2008; Brandenburg et al. 2009; Guerrero et al. 2010). Even though the helicity constraint in direct numerical simulations (DNS) of dynamos with strong shear have been clearly identified, the results can be matched with mean field models with a weaker algebraic quenching than α² dynamos (Brandenburg et al. 2001). It is possible to include this process in mean-field dynamo models through an equation describing the evolution of the small-scale current helicity. We shall refer to this equation as the dynamical quenching mechanism.

In this paper we perform a series of calculations with mean field αΩ models in spherical geometry along with a dynamical equation for the evolution of α for magnetic Reynolds numbers in the range 1 ≤ R_m ≤ 2 × 10⁵. An important feature of the calculation is that the region of strong narrow shear is separated from the region of helical turbulence. In addition to providing detailed results not mentioned in Chatterjee et al. (2010), this paper is also aimed at studying somewhat more complicated models like flux transport (FT) models, which include meridional circulation (MC). The role of diffusive helicity fluxes modelled into the dynamical quenching equation by using a Fickian diffusion term is also discussed for various models. Helicity fluxes across an equator can indeed be modelled by this diffusion term, as was shown by Mitra et al. (2010). In Sects. 2.1 and 2.2 we discuss the features of the two classes of αΩ models used. The formulation of dynamical α quenching is given in Sect. 2.3. The results for 1 < R_m < 2 × 10⁵, with or without small-scale helicity fluxes, are presented in Sects. 3.1–3.3, while Sect. 3.4 presents the results for the flux transport dynamo model. Finally, we draw the conclusions of this study in Sect. 4.

2. αΩ dynamo models with dynamical α quenching

2.1. Interface dynamo

We solve the induction equation in a spherical shell assuming axisymmetry. Our dynamo model consists of the induction equations for the mean poloidal potential A_φ(r,θ) and the mean toroidal field B_φ(r,θ) (see Eqs. (2) and (3) of Guerrero et al. 2010, GCB from now on). Axisymmetry demands that for all variables ∂/∂φ = 0. The turbulent coefficients α and η_t in these equations are estimated as follows: from the mixing length theory we know that the turbulent diffusivity in the convection zone is given by (cf. Sur et al. 2008) ${\mathit{\eta }}_{\mathrm{t}}\mathrm{=}\frac{{\mathit{u}}_{\mathrm{rms}}}{\mathrm{3}{\mathit{k}}_{\mathrm{f}}}\mathit{,}$(1)where u_rms is the rms velocity of the turbulent eddies, k_f is the wavenumber of the energy-carrying eddies, corresponding to the inverse pressure scale height near the base of the convection zone. In our models, we have taken k_f = 10k₁, where k₁ is the longest wave that can fit the domain latitudinally. Below the convection zone the magnetic diffusivity has the molecular value, η_r. Both regions are smoothly matched by the step function ${\mathrm{\Theta }}^{\mathrm{\pm }}\mathrm{\left(}\mathit{r,}{\mathit{r}}_{\mathrm{e}}\mathit{,}{\mathit{d}}_{\mathrm{e}}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{\pm }\mathrm{erf}\left(\frac{\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{e}}}{{\mathit{d}}_{\mathrm{e}}}\right)\mathit{,}$(2)where r_e = 0.73   R_⊙, and d_e = 0.025   R_⊙. Thus the magnetic diffusivity profile is given by $\mathit{\eta }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\eta }}_{\mathit{r}}\mathrm{+}{\mathit{\eta }}_{\mathrm{t}}{\mathrm{\Theta }}^{\mathrm{+}}\mathrm{\left(}\mathit{r,}{\mathit{r}}_{\mathrm{e}}\mathit{,}{\mathit{d}}_{\mathrm{e}}\mathrm{\right)}\mathit{.}$(3)The turbulent kinetic α effect considered in this model has the form ${\mathit{\alpha }}_{\mathrm{K}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{0}}{\mathrm{\Theta }}^{\mathrm{+}}\mathrm{\left(}\mathit{r,}{\mathit{r}}_{\mathrm{a}}\mathit{,}{\mathit{d}}_{\mathrm{a}}\mathrm{\right)}\cos \mathit{\theta }}{\mathrm{1}\mathrm{+}{\mathit{g}}_{\mathit{\alpha }}{\overline{\mathit{B}}}^{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{eq}}^{\mathrm{2}}}\mathit{,}$(4)where g_α is a non-dimensional coefficient equal to 1 or R_m depending on the assumed form of algebraic quenching in the models and ${\mathit{B}}_{\mathrm{eq}}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{\pi \rho }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{u}}_{\mathrm{rms}}\mathrm{=}\mathrm{(}\mathrm{4}\mathit{\pi \rho }{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{3}{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}_{\mathrm{f}}$is the value of the magnetic field at equipartition between the magnetic and the turbulent kinetic energies. The amplitude of the kinetic α effect is estimated by considering the first order smoothing approximation (FOSA) as being α₀ = τϵ_fω_rmsu_rms/3, where ω_rms is the rms vorticity of the turbulence and τ ~ (k_fu_rms)^-1 is the eddy correlation time scale. The prefactor ϵ_f, usually of the order of 0.1 or less, is used because (u·ω)_rms < u_rmsω_rms. The case ϵ_f = 1 means that the flow is maximally helical. These approximations give us an estimate of α₀ in terms of eddy diffusivity η_t and forcing scale k_f as ${\mathit{\alpha }}_{\mathrm{0}}\mathrm{=}{\mathit{ϵ}}_{\mathrm{f}}\frac{\mathit{\tau }{\mathit{k}}_{\mathrm{f}}{\mathit{u}}_{\mathrm{rms}}^{\mathrm{2}}}{\mathrm{3}}\mathrm{=}{\mathit{ϵ}}_{\mathrm{f}}{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}_{\mathrm{f}}\mathit{.}$We would consider ϵ_f instead of α₀ as a free parameter in the model apart from η_t. Even though the helical turbulence pervades almost the entire convection zone, here we take r_a = 0.77   R_⊙ and d_a = 0.015   R_⊙, so that we can have a large spatial separation between the shear and turbulent layers. Consequently we consider a differential rotation profile like that in the high latitude tachocline of the Sun given by $\mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{0}}{\mathrm{\Theta }}^{\mathrm{+}}\mathrm{\left(}\mathit{r,}{\mathit{r}}_{\mathrm{w}}\mathit{,}{\mathit{d}}_{\mathrm{w}}\mathrm{\right)}\mathit{,}$(5)where Ω₀ = 14 nHz, r_w = 0.68   R_⊙ and d_w = 0.015   R_⊙. The radial profiles of η_t, α and ∂Ω/∂r are plotted as a function of fractional radius r/R_⊙ in Fig. 1. In that figure it is possible to see that the region of strong radial shear is separated from the region of helical turbulence and the diffusivity has a strong gradient at a radius lying between these two source layers. The reason of this is to decrease the cycle period T_cyl of the oscillatory dynamos to a reasonably small fraction of the diffusion time t_diff.

Fig. 1 Profiles of radial shear ∂Ω/∂r, α and η as a function of fractional solar radius.

2.2. Flux transport dynamo

Apart from the interface dynamo, we also run numerical experiments for models that include meridional circulation and consider the so-called Babcock-Leighton (BL) α effect. In these models the dynamo is thought to be operating in two separated regions, the toroidal component of the magnetic field is produced at the base of the convection zone and the poloidal field is produced by the decay of tilted bipolar active regions on the solar surface. The MC plays an important role in the dynamics of the system by advecting the magnetic flux and connecting both source regions. For this reason, these models are often referred to as flux-transport (FT) dynamo models. In the literature FT models have been studied extensively by several authors (Dikpati & Charbonneau 1999; Chatterjee et al. 2004; Guerrero & Dal Pino 2008, and references therein). These models have now reached a stage where they are able to reproduce the butterfly diagram as well as the correct phase between the polar fields and the toroidal fields.

We recall that the α_K is now not due to the helical turbulence in the bulk of the convection zone, but due to a phenomenological BL α where the poloidal field is produced from the toroidal field by decay of tilted bipolar active regions. For the models in this section, the BL α is assumed to be concentrated only in the upper 0.05% of the convection zone and becomes maximum at  ± 40° latitude and goes to zero at poles. The expression for α_K is $\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{K}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}{\mathit{\alpha }}_{\mathrm{BL}}{\mathrm{\Theta }}^{\mathrm{+}}\mathrm{\left(}\mathit{r,}\mathrm{0.95} {\mathit{R}}_{\mathrm{\odot }}\mathit{,d}\mathrm{\right)}{\mathrm{\Theta }}^{\mathrm{-}}\mathrm{\left(}\mathit{r,}{\mathit{R}}_{\mathrm{\odot }}\mathit{,d}\mathrm{\right)}\cos \mathit{\theta }{\sin }^{\mathrm{2}}\mathit{\theta ,}& & \end{array}$(6)where d = 0.015   R_⊙.

For the meridional circulation we use the analytical expressions given by van Ballegooijen & Choudhuri (1988), where the radial and latitudinal components of the velocity field are given by ${\mathit{u}}_{\mathit{r}}\mathrm{=}{\mathit{u}}_{\mathrm{0}}{\left(\frac{{\mathit{R}}_{\mathrm{\odot }}}{\mathit{r}}\right)}^{\mathrm{2}}\mathit{\zeta }\left(\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{s}\mathrm{1}}}{\mathrm{2}}{\mathit{\zeta }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}\frac{\mathrm{4}{\mathit{c}}_{\mathrm{s}\mathrm{2}}}{\mathrm{9}}{\mathit{\zeta }}^{\mathrm{3}\mathit{/}\mathrm{4}}\right)\mathrm{(}\mathrm{2}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{-}{\sin }^{\mathrm{2}}\mathit{\theta }\mathrm{)}\mathit{,}$(7)${\mathit{u}}_{\mathit{\theta }}\mathrm{=}{\mathit{u}}_{\mathrm{0}}{\left(\frac{{\mathit{R}}_{\mathrm{\odot }}}{\mathit{r}}\right)}^{\mathrm{3}}\mathrm{(}\mathrm{-}\mathrm{1}\mathrm{+}{\mathit{c}}_{\mathrm{s}\mathrm{1}}{\mathit{\zeta }}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathit{c}}_{\mathrm{s}\mathrm{2}}{\mathit{\zeta }}^{\mathrm{3}\mathit{/}\mathrm{4}}\mathrm{)}\sin \mathit{\theta }\cos \mathit{\theta ,}$(8)where ζ = R_⊙/r−1, r_b = 0.71R_⊙, ζ_b = R_⊙/r_b−1, ${\mathit{c}}_{\mathrm{s}\mathrm{1}}\mathrm{=}\mathrm{4}{\mathit{\zeta }}_{\mathrm{b}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$ and ${\mathit{c}}_{\mathrm{s}\mathrm{2}}\mathrm{=}\mathrm{3}{\mathit{\zeta }}_{\mathrm{b}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{4}}$. In the equations above, the flow is poleward at the surface with a maximum amplitude of u₀ = 20 m s^-1, and the equatorward return flow occurs at r_b with an amplitude of around 3 m s^-1. Unlike in flux transport dynamo models, the meridional circulation does not reverse the direction of propagation of the dynamo wave in interface dynamo models as long as the meridional circulation is confined within the convection zone (Petrovay & Kerekes 2004). The shear is still radial and given by Eq. (5) with r_w = 0.7   R_⊙. Finally, the turbulent diffusivity has the same profile as in Eq. (3), but with η_t = 2 × 10¹¹ cm s^-1 and r_e = 0.7   R_⊙.

2.3. Dynamical α quenching

Pouquet et al. (1976) first showed that the turbulent α effect is modified due to the generation of small-scale helicity in the way given by Eq. (9) below. The second term is sometimes referred to as the magnetic α-effect, $\mathit{\alpha }\mathrm{=}{\mathit{\alpha }}_{\mathrm{K}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{M}}\mathrm{=}\mathrm{-}\frac{\mathit{\tau }}{\mathrm{3}}\mathrm{(}\overline{{\omega }\mathrm{·}{u}}\mathrm{-}{\mathit{\rho }}^{-1}{\overline{{j}\mathrm{·}{b}}}^{\mathrm{)}}\mathit{,}$(9)where ω, u, j, b denote the fluctuating component of the vorticity, velocity, current, and magnetic field in the plasma. Also, b = ∇ × a. It is possible to write an equation for the evolution of the magnetic part of α or α_M from the equation for the evolution of the small-scale magnetic helicity density ${\mathit{h}}_{\mathrm{f}}\mathrm{=}\overline{{a}\mathrm{·}{b}}$ using the relation ${\mathit{\alpha }}_{\mathrm{M}}\mathrm{=}\frac{{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}_{\mathrm{f}}^{\mathrm{2}}}{{\mathit{B}}_{\mathrm{eq}}^{\mathrm{2}}}{\mathit{h}}_{\mathrm{f}}\mathit{.}$(10)As the equation for $\overline{{a}\mathrm{·}{b}}$ is gauge-dependent, it makes sense only to write an equation for the volume-averaged quantity in order to avoid dependence on a specific gauge (Blackman & Brandenburg 2002). Our dynamo equations are independent of any gauge because we solve for the magnetic potential component A_φ with an axisymmetric constraint. It is important for us that the equation for α_M is also gauge-independent. Subramanian & Brandenburg (2006) used the Gauss linking formula for the expression for h_f and wrote an equation independent of the gauge for the magnetic helicity density under the assumption that the correlation length for all fluctuating variables remains small compared to the system size at all times. Using Eq. (10) we write the same equation in terms of α_M, $\frac{\mathit{\partial }{\mathit{\alpha }}_{\mathrm{M}}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\mathrm{2}{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}_{\mathrm{f}}^{\mathrm{2}}\left(\frac{\overline{{ℰ}}\mathrm{·}\overline{{B}}}{{\mathit{B}}_{\mathrm{eq}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\alpha }}_{\mathrm{M}}}{{\mathit{R}}_{\mathrm{m}}}\right)\mathrm{-}{\nabla }\mathrm{·}{\overline{{F}}}_{\mathit{\alpha }}\mathit{,}$(11)where $\overline{{ℰ}}$ and $\overline{{B}}$ are the mean field EMF and the mean magnetic field. We recall that in the mean field approach, we do not solve evolution equations for the fluctuating components a and b. Instead we model the effect of the small scales by writing an evolution equation for h_f similar to Eq. (11). Because the evolution of the magnetic field components (A_φ, B_φ) in the large scales is explicitly solved, we do not need a corresponding equation for the large-scale helicity ($\overline{{A}}\mathrm{·}\overline{{B}}$). Moreover, $\overline{{a}\mathrm{·}{b}}$ and $\overline{{A}}\mathrm{·}\overline{{B}}$ are in principle gauge-dependent. However, if there is scale separation, $\overline{{a}\mathrm{·}{b}}$ can be written as a density of linkages and is therefore gauge-independent (Subramanian & Brandenburg 2006), while $\overline{{A}}\mathrm{·}\overline{{B}}$ is not. This was demonstrated in Hubbard & Brandenburg (2010).

The flux term, ${\overline{{F}}}_{\mathit{\alpha }}$, consists of individual components, e.g., advection owing to the mean flow, Vishniac-Cho fluxes (Vishniac & Cho 2001), effects of mean shear, diffusive fluxes, etc. When we make an ansatz for the flux of normalized small-scale helicity, ${\overline{{F}}}_{\mathit{\alpha }}$, the effect on the large-scale helicity flux is also taken into account by the nonlinear coupling between mean-field induction equations and Eq. (11). If small-scale magnetic helicity is removed by stellar winds in Eq. (11), a corresponding advection term in the mean field induction equation also changes the large-scale helicity accordingly.

For the interface models we put F_α = 0 unless mentioned otherwise. For the flux transport dynamo (Sect. 2.2), the normalized helicity flux F_α in Eq. (11) is given by ${\overline{{F}}}_{\mathit{\alpha }}\mathrm{=}{\mathit{\alpha }}_{\mathrm{M}}{{u}}_{\mathrm{p}}\mathrm{-}\mathit{\kappa }{\nabla }{\mathit{\alpha }}_{\mathrm{M}}\mathit{,}$(12)where κ = κ₀η(r) is the diffusion coefficient for α_M. Let us define the diffusion time in the model as ${\mathit{t}}_{\mathrm{diff}}\mathrm{=}\mathrm{1}\mathit{/}{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}$. The decay time in Eq. (11) is then ${\mathit{t}}_{\mathit{\alpha }}\mathrm{=}{\mathit{R}}_{\mathrm{m}}\mathit{/}{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}_{\mathrm{f}}^{\mathrm{2}}\mathrm{=}\mathrm{4.55}\mathrm{\times }{\mathrm{10}}^{-3}{\mathit{R}}_{\mathrm{m}}{\mathit{t}}_{\mathrm{diff}}$. Note that we use g_α = 0 in Eq. (4) whenever we employ the dynamical quenching equation, because dynamical quenching is usually more important.

Our computational domain is defined to be the region confined by 0 ≤ θ ≤ π and 0.55   R_⊙ ≤ r ≤ R_⊙. Unless otherwise stated, the boundary conditions for A_φ are given by a potential field condition at the surface (Dikpati & Choudhuri 1994) and A_φ = 0 at the poles. We also performed some calculations with the vertical field condition at the top boundary, which means that B_θ = B_φ = 0. At the bottom we use the perfect conductor boundary condition of Jouve et al. (2008) with A_φ = ∂(rB_φ)/∂r = 0. But a more realistic perfect conductor boundary condition in our opinion would be ∂(rB_θ)/∂r = ∂(rB_φ)/∂r = 0. Also B_φ = 0 on all other boundaries. The equation for α_M is an initial value problem for F_α = 0. For finite fluxes we also set α_M = 0 at all boundaries. We checked that the results are not very sensitive to the different bottom boundary conditions given above, mainly because the bottom boundary is far removed from the dynamo region.

3. Results

3.1. Nonlinear interface dynamo without helicity fluxes

In order to study the R_m dependence of the saturation of the magnetic field in the interface dynamo of Sect. 2.1, we keep all the dynamo parameters the same for all runs and change η_r from 2 × 10⁵ cm² s^-1 to 2 × 10¹⁰ cm² s^-1 while keeping η_t fixed at 4 × 10¹⁰ cm² s^-1.

Fig. 2 Critical α in terms of a fraction of ηtkf as a function of magnetic Reynolds number Rm for the interface dynamo model of Fig. 1.

To be able to correctly compare the dynamo models for different R_m, it is important to first calculate the critical value of α₀, denoted by α_c for each model. This is done by solving the mean field dynamo equations with an algebraic quenching with g_α = 1. The values of α_c as a function of R_m are shown in Fig. 2. From this figure we observe that the model is most efficient for R_m ~ 20, because for this value the critical dynamo number is the lowest. A similar variation of α_c with the ratio η_t/η_r was obtained analytically for interface dynamos by MacGregor & Charbonneau (1997, see their Fig. 5A). All the simulations from now on are performed by setting α₀ = 2α_c, corresponding to the R_m for each model. The cycle period of these dynamo models (T_cyl) slightly increases with R_m and is found to be of the order of 2t_diff, for the investigated range of R_m. Figure 3 shows the evolution of the volume-averaged magnetic energy as a function of time for an algebraic form of quenching with g_α = 1. Note that in Fig. 3, the slopes in the kinematic phase are almost similar for all R_m within the errors in the numerical determination of α_c.

Fig. 3 Time evolution of the volume-averaged magnetic energy in the domain scaled with the equipartition energy for the case with algebraic quenching and gα = 1 for models with Rm = 1 (diamond), Rm = 20 (solid), Rm = 200 (dashed), Rm = 2 × 103 and Rm = 2 × 105 (triangles).

Fig. 4 Same as Fig. 3 but with dynamical quenching and gα = 0 for the Rm indicated in the figure.

Fig. 5 Saturated value of the volume-averaged magnetic energy scaled with the equipartition energy as a function of Rm for models with dynamical α quenching (triangles+solid) and algebraic quenching with gα = 1 (squares + dashed) and with gα = Rm (cross + dashed-dotted). Note that all runs included in this figure used α0 = 2αc.

Now, we set g_α = 0 and consider a dynamically quenched α by solving Eq. (11) along with the linear induction equations for the mean fields A_φ and B_φ. The time evolution of the volume-averaged magnetic energy for these systems for a range of magnetic Reynolds numbers is shown in Fig. 4. The strong R_m dependence, which is reminiscent of catastrophic quenching in all astrophysical dynamos, can be easily discerned from this figure. For high values of R_m, the saturation phase is clearly different from that obtained for algebraic quenching (see Fig. 3).

In this context it is important to recall an important difference between dynamos in periodic and open domains such as those considered here. In a periodic domain, the dynamo is expected to reach and sustain mean fields of the order of the inverse square root of the scale separation ratio (Blackman & Brandenburg 2002). However, in an open domain the final saturation field strength is always of the order of ${\mathit{R}}_{\mathrm{m}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$ (Brandenburg & Dobler 2001; Brandenburg & Subramanian 2005), although the mean field may reach an early equipartition strength peak just at the time when the small-scale field reaches saturation (Fig. 3). The physical reality of this peak remains somewhat mysterious, although direct simulations also sometimes show this peak (Hubbard & Brandenburg 2010). In Fig. 5 we summarize the results of models with a simple algebraic quenching (Eq. (4)) for both g_α = 1 and g_α = R_m, as well as models with dynamical quenching (Eq. (11)). Both algebraically quenched models with g_α = R_m and those with dynamical α quenching exhibit a monotonic decrease of the saturation magnetic energy with R_m.

In our models we see that the region of strong toroidal field, B_φ, is different from that of the poloidal field produced by the α effect as shown in Fig. 6. We verified that the same dependence of the B_sat on R_m is reproduced in models where the two source regions are further spatially separated by setting r_a = 0.87   R_⊙ instead of 0.77   R_⊙.

Fig. 6 Radial profiles of Aφ and Bφ at two different latitudes (λ) in the saturated phase for Rm = 2 × 103.

The nature of the saturation curves of the magnetic energy (Fig. 4) is strongly governed by the ratio of t_α and T_cyl. For low R_m, t_α ≤ T_cyl, so that the time evolution of $\mathrm{⟨}{\overline{\mathit{B}}}^{\mathrm{2}}\mathrm{⟩}\mathit{/}{\mathit{B}}_{\mathrm{eq}}^{\mathrm{2}}$ is flat in the saturation regime – similar to the algebraic quenching case (see Fig. 3). On the other hand, for R_m = 2 × 10³, the decay time t_α ≫ T_cyl and so the system is underdamped, i.e., there are amplitude modulations of the magnetic energy before it settles to a final saturation value (see dash-dotted line in Fig. 4). For R_m = 2 × 10⁵, the underdamped oscillation has a long period (t_α) so that the model has to be run for more than 500 t_diff before the dynamo field starts becoming “strong” again. Because of the long computational time involved in this exercise we have not continued the calculation beyond 60 t_diff. Hence, the determination of saturation magnetic energy may be inaccurate.

Fig. 7 Time-latitude plot of αM(0.72   R⊙,θ) for a) Rm = 20 and b) Rm = 200.

A similar pattern is observed in the time-latitude plot of α_M as shown in Fig.  7. In this figure we note that for R_m = 20 (upper panel), α_M shows strong oscillations (see the modulation in the colour of the filled contours) since t_α ≪ T_cyl. For R_m = 200 (lower panel), t_α ~ T_cyl, so that the amplitude of the oscillation is weak because the α_M decays at the same rate at which it is produced. From the same figure it may be concluded that the small-scale current helicity, α_M, is predominantly negative (positive) in the Northern (Southern) hemisphere. For high values of R_m (≳2000) the models start showing changes in parity for t > 40t_diff. Nevertheless, the magnetic energy and the dynamo period, T_cyl, remain fairly constant even while the system fluctuates between symmetric and anti-symmetric parity at irregular time intervals (see Fig. 8). This parity oscillation is absent in the corresponding models with algebraic quenching.

Fig. 8 a) Evolution of parity (purely dipolar  = −1 and purely quadrupolar  =  +1) for an interface dynamo model with Rm = 2 × 103 and no fluxes. b) A small part in the time-latitude diagram of Bφ indicated by dotted lines in a) where the parity is changing from quadrupolar to dipolar.

Fig. 9 Time-latitude plot of the toroidal field a) and c) and αMb) and d) for interface dynamo models with α = 4αc and Rm = 20.

3.2. Secondary dynamo waves with F_α = 0

An interesting result emerges when we increase the value of the kinetic α effect (α₀ = 4α_c instead of 2α_c as in Sect. 3.1) in the interface dynamo model with dynamical α quenching and R_m = 20. In this case we observe that in addition to the primary dynamo wave travelling equatorward, a poleward-propagating secondary dynamo wave appears in the butterfly diagram for B_φ at r = 0.72   R_⊙ (see finger-like projections at high wave number in Fig. 9a). A weak signature of this secondary wave can also be seen in the butterfly diagram at 0.8   R_⊙ in Fig. 9b.

The existence of these secondary dynamo waves can be explained by studying the source term in Eq. (11), i.e., $\overline{{ℰ}}\mathrm{·}\overline{{B}}\mathrm{=}\mathit{\alpha }{\overline{{B}}}^{\mathrm{2}}\mathrm{-}{\mathit{\eta }}_{\mathrm{t}}\overline{{J}}\mathrm{·}\overline{{B}}$. Because of a low value of R_m, the toroidal field generated in the shear layer diffuses into the overshoot layer, and by virtue of the term $\mathrm{-}{\mathit{\eta }}_{\mathrm{t}}\overline{{J}}\mathrm{·}\overline{{B}}$, generates a region of α_M that is positive (negative) in the northern (southern) hemisphere. Because α_K → 0 for r ≤ 0.73   R_⊙, the total α effect there has the sign of α_M. The new dynamo wave therefore travels poleward according to the Parker-Yoshimura rule.

Vishniac & Cho (2001) also presented a case of an αΩ dynamo driven by a supercritical helicity flux. Their mechanism requires a finite initial magnetic field unlike here, where the initial field is  ~10^-6B_eq. The difference compared to the case above is that the mean field dynamo is not driven by supercritical Vishniac & Cho fluxes, but by a local generation of small-scale magnetic helicity. Note that in this case the secondary dynamo wave is energetically powered by the kinematic part of the helical convection. Blackman & Field (2004) also found a magnetic-helicity driven dynamo (MHDD) in addition to a kinetic-helicity driven dynamo (KHDD) by solving a coupled set of equations for large-scale helicity and the mean small-scale helicity. They highlight an important difference between KHDD and MHDD, namely that in the former, the α_K produces small-scale and large-scale magnetic helicity of opposite signs, whereas in the latter magnetic helicity has same sign in both scales. Our calculations in this section agree with their observation in the sense that α_M has the same sign as $\overline{{J}}\mathrm{·}\overline{{B}}$. Evidence for dominance of magnetically generated α in stratified magnetorotational turbulence in accretion discs has also been found by Gressel (2010).

We have not observed any evidence of chaotic behaviour in the range of magnetic Reynolds number 20 ≤ R_m ≤ 2 × 10⁵ for supercritical α ≤ 4α_c in agreement with Covas et al. (1997). However, if the α effect is highly supercritical, the dynamical quenching formula for α_M is insufficient for dynamo saturation, and additional algebraic quenching terms are needed (Kleeorin & Rogachevskii 1999). We continue the discussion of secondary dynamo waves driven by diffusive magnetic helicity fluxes in the next section.

3.3. Supercritical diffusive magnetic helicity fluxes

Recently, Brandenburg et al. (2009) showed that catastrophic quenching in one-dimensional α² dynamos can be alleviated by introducing a Fickian diffusive flux in Eq. (11) given by ${\overline{{F}}}_{\mathit{\alpha }}\mathrm{=}\mathrm{-}\mathit{\kappa }{\nabla }{\mathit{\alpha }}_{\mathrm{M}}\mathit{.}$(13)There was an attempt to calculate the diffusion coefficient from direct numerical simulations and it was found that κ ~ 0.3η_t for R_m ~ 20 (Mitra et al. 2010). In the context of of αΩ dynamos GCB showed that even a very small diffusive flux for the mean small-scale helicity alleviates catastrophic quenching for a scale separation k_f/k₁ ~ 10. Their Fig. 7 moreover shows that ${\overline{\mathit{B}}}_{\mathrm{sat}}\mathit{/}{\mathit{B}}_{\mathrm{eq}}$ levels off at high R_m for finite diffusive fluxes even though the value of ${\overline{\mathit{B}}}_{\mathrm{sat}}$ at which the curve levels off may vary with k_f.

Fig. 10 a) Time evolution of the volume averaged magnetic energy in the domain scaled with the equipartition energy for models with Rm = 2 × 103 and κ0 = 10-5 (∗ ) as well as with Rm = 2 × 105 and κ0 = 10-5 (solid line) for the interface dynamo model of Sect. 3.1 with dynamical α quenching. The saturation curve for zero fluxes for the model with Rm = 2 × 105 has been shown by the dashed line; b) and c) show time-latitude diagrams for the toroidal field at the depths indicated for models with Rm = 2 × 103 and κ0 = 10-5.

Fig. 11 The time-latitude plot of toroidal field a) and c), and αMb) and d) with α = 2αc for Rm = 2 × 105 and κ0 = 0.01. Meridional snapshots of e) sign(Bφ)(| Bφ | /Beq)1/2 and f) αm × 103 for the same case.

For κ = 0, the saturation curves in Fig. 4 show that the B_sat either displays underdamped oscillations for R_m = 2 × 10³ or goes through very low values for R_m ~ 2 × 10⁵ and takes a long time to relax to a steady amplitude. In this section we introduce a diffusive flux with κ(r) = κ₀η(r). Figure 10 shows the results for models with two different values of R_m and κ. For R_m = 2 × 10³ and κ₀ = 10^-5, the saturation energy is comparable to the corresponding case with κ = 0 (see dash-dotted line in Fig. 4). But now grand minima-like episodes with respect to the primary mode appear in the system, as may be seen in the oscillations in the volume-averaged energy (represented by  ∗ ’s in Fig. 10a) with a period  ~5 times the period of the equatorward-propagating mode long after saturation. Another interesting behaviour can be discerned from the butterfly diagram of the toroidal field (Fig. 10b, c). It appears that the dynamo here is governed by the competition between the equatorward-propagating primary mode and the poleward-propagating secondary mode.

For R_m = 2 × 10⁵ and κ = 0.01, we observe underdamped oscillations with a final B_sat ~ 0.1B_eq. On studying the corresponding butterfly diagrams (Fig. 11a–d) we find a poleward- propagating mode caused by the radial diffusion of the α_M into the stable layers. This was not possible for the model with the same R_m but κ = 0. Figures 11e, f show meridional snapshots of sign(B_φ)(| B_φ | /B_eq)^1/2 and α_M in order to get a clear idea of the distribution of magnetic fields.

The poleward-propagating mode is now driven by supercritical diffusive helicity fluxes, as opposed to supercritical Vishniac & Cho fluxes (see Brandenburg & Subramanian 2005, for examples of this behaviour). For the same model (R_m = 2 × 10⁵) there exists a critical κ_c ~ 10^-5 such that the secondary dynamo wave does not appear if κ₀ < κ_c. The volume-averaged magnetic energy in this case decays eventually. Note that this threshold for κ is highly dependent on R_m. For instance a model with R_m = 2 × 10³ and κ₀ = 10^-5 produces a dynamo with finite saturation magnetic energy and dynamo wave propagation governed by α_M, whereas for κ₀ = 10^-4, the dynamo shows a runaway growth.

Unlike in GCB, we used super-critical helicity fluxes in this section. Note that the critical value, κ_c, is much lower than that used in GCB where the maximum value of κ₀ used is  ~10. Interface dynamo models are a particular case of αΩ models where a strong toroidal field is produced below the convection zone. Even a low flux of α_M into the stable layers can produce very large magnetic fields and power a secondary dynamo. Direct numerical simulations of α² dynamos have established that a large-scale magnetic field is easily excited on the scale of the system i.e., ${\mathit{k}}_{\mathrm{1}}^{-1}$ for a high k_f/k₁ ratio (Archontis et al. 2003). However the length scale of the magnetic field in Figs.10b, c and 11a, c, e is comparable to ${\mathit{k}}_{\mathrm{f}}^{-1}$, which suggests that the degree of scale separation may have become insufficient to write the electromotive force as a simple multiplication, as is done in the expression $\overline{{ℰ}}\mathrm{=}\mathit{\alpha }\overline{{B}}\mathrm{-}{\mathit{\eta }}_{\mathrm{t}}\overline{{J}}$. Then it may become necessary to write the electromotive force as a convolution, which essentially corresponds to a low-pass filter (see, e.g., Brandenburg et al. 2008). However, we have not pursued this aspect any further.

3.4. Flux transport Babcock-Leighton dynamo

In this section we perform numerical experiments of the flux transport dynamo model described in Sect. 2.2. As in Sect. 3.1, we first find the critical α_BL required to have a self-excited dynamo. In this case α_c = 5.1 m s^-1 for R_m = 2 × 10³. We pursue the rest of the calculations with α_BL = 6.0 m s^-1 and including dynamical α quenching. This slightly supercritical choice of α_BL avoids the production of very large α_M leading to secondary dynamos as discussed in Sect. 3.2. We should emphasize that Eq. (9) represents a first order correction to α and should be treated with caution while in supercritical regimes.

Fig. 12 Time evolution of the volume averaged magnetic energy for the flux transport dynamo model of Sect. 3.2 for Rm = 2 × 103 with κ0 = 0.3 (dashed); Rm = 2 × 105 with κ0 = 0.3 (solid); Rm = 2 × 103 with κ = 0 (dashed-dotted); Rm = 2 × 105 with κ0 = 0 (diamond+dashed).

At first we do not consider any small-scale helicity fluxes in the model, i.e., F_α = 0 in Eq. (11). The saturation curve for R_m = 2 × 10³ in Fig. 12 is now under-damped, whereas the dynamo fails to generate a finite B_sat for R_m = 2 × 10⁵ even though it initially has the same growth rate. On increasing α_BL = 10 ms^-1 from 6 ms^-1 the saturation curve for R_m = 2 × 10⁵ also displays underdamped behaviour. This indicates that for α_BL = 6 ms^-1, the total α in the domain was simply becoming sub-critical and the dynamo was not able to sustain itself through the saturation phase. In Fig. 13a, b, we show the meridional snapshots of the magnetic fields and α_M for a saturated state. In this figure, the colour bars indicate that α_M increases while B_φ decreases with R_m for the same value of α_BL.

Now we include an advective flux of α_M caused by the meridional circulation in Eq. (11). In order to keep the system numerically stable, we also require a diffusive flux in Eq. (11). It is clear from Fig. 12 that the underdamped behaviour in the model without fluxes is suppressed due to a diffusive flux of α_M. The production of α_M is now countered by a diffusive decay in a time scale much shorter than that given by ${\mathit{t}}_{\mathit{\alpha }}\mathrm{=}{\mathit{R}}_{\mathrm{m}}\mathit{/}{\mathit{\eta }}_{\mathrm{t}}{\mathit{k}}_{\mathrm{f}}^{\mathrm{2}}$. The same figure shows that the saturation value of the magnetic energy is now almost independent of R_m. With diffusive and advective fluxes owing to the meridional circulation in Eq. (11) the small-scale helicity is distributed through out the convection zone as shown in Fig. 14a, b. It is instructive to compare the magnitudes of α_M in Figs. 13 and 14.

The diffusive flux of small-scale helicity is therefore crucial for the operation of a successful mean field αΩ dynamo. An interesting observation is the distribution of α_M in the concentrated region at the lower part of the convection zone (see Fig. 13) in contrast to Fig. 14. Even though α_BL is a surface phenomena, considerable magnetic helicity is generated in the entire convection zone when the meridional circulation sinks the poloidal field lines at high latitudes and brings them near the tachocline where the toroidal fields are generated.

Fig. 13 Meridional cross-sections showing the distribution of toroidal field and αM for a Babcock-Leighton dynamo without MC and diffusive helicity fluxes in Eq. (11) for a) Rm = 2 × 103 and b) Rm = 2 × 105. The streamlines of the positive and negative poloidal field are shown by solid and dashed lines respectively. Note that the magnetic field has decayed to very small values for Rm = 2 × 105.

Fig. 14 Meridional cross-sections showing the distribution of toroidal field and αM for a Babcock-Leighton dynamo with MC and diffusive helicity fluxes for Rm = 2 × 103 at two different epochs. The streamlines of the positive and negative poloidal field are shown by solid and dashed lines respectively.

4. Conclusions

We have performed calculations for αΩ dynamos in a spherical shell for spatially segregated α and Ω source regions. The two classes of models we studied resemble Parker’s interface dynamo and the flux transport dynamo with a Babcock-Leighton α.

In agreement with earlier work (Chatterjee et al. 2010), we found that it is not possible to escape the catastrophic α quenching by merely separating the regions of shear and α-effect. The saturation value of magnetic energy in the interface model decreases as $\mathrm{~}{\mathit{R}}_{\mathrm{m}}^{-1}$ for both dynamical quenching with g_α = 0 and the algebraic quenching with g_α = R_m (Fig. 5). Nevertheless, we found that a richer dynamical behaviour emerges for the cases with dynamical α effect, in terms of parity fluctuations (Fig. 8) and the appearance of “secondary” dynamos (Fig. 9), namely dynamo waves governed by the magnetic α effect. We observe this interesting behaviour in two different situations. Firstly, for low R_m, the toroidal field is weakly confined below the overshoot layer and therefore α_M is generated there by the term $\mathrm{-}{\mathit{\eta }}_{\mathrm{t}}\overline{{J}}\mathrm{·}\overline{{B}}$ in the expansion of the forcing term $\overline{{ℰ}}\mathrm{·}\overline{{B}}$ in Eq. (11). Secondly, for high R_m, the secondary dynamos are only observed if there is a supercritical diffusive helicity flux where the α_M diffuses into the stable layers. However, because of the lack of scale separation between the mean field and the forcing scale of the helical turbulence we refrain from interpreting this in terms of the poleward migration seen in the Sun. We have to be cautious about using the dynamical quenching equation for dynamo numbers that are not very large compared to the critical dynamo number. For highly supercritical α_K, the behaviour of the system begins to be governed by α_M. According to our knowledge, there are no direct numerical simulation of secondary dynamos or magnetic-helicity driven dynamos applied to stellar magnetic fields. This will be explored in a future paper.

We do not see any evidence for chaotic behaviour in the time series of magnetic energy because the dynamo period and the saturation energy remain fairly constant. This may not be the case for diffusive helicity fluxes, which introduce further complexity to the system. Yet the addition of a diffusive helicity flux relaxes the catastrophic ${\mathit{R}}_{\mathrm{m}}^{-1}$ dependence of the saturation magnetic energy (Figs. 10a, 12).

We would expect that the magnetic field affects all turbulent coefficients including α and η. But for this analysis we did not included an equation for the variation of η_t. This is justified for the simple two layer model with a lower η_t in the region of the production of strong toroidal fields and a higher η_t in the region of weaker poloidal fields. Note also that by quenching the diffusivity inversely with magnetic energy in a nonlinear dynamo model, Tobias (1996) was able to produce a bonafide interface model where the magnetic field was restricted to a thin layer at an interface between a layer of shear and cyclonic turbulence. None of the previous interface models have used the dynamical quenching equation.

In the flux transport dynamo model, the magnetic energy fails to reach significant saturation when both the meridional circulation and the diffusive helicity fluxes are artificially turned off in the helicity evolution equation. This is expected and is demonstrated in Fig. 12. It is interesting that the Babcock-Leighton dynamos, where α is concentrated only in a narrow layer at the surface, also produce considerable amount of magnetic helicity inside the convection zone when dynamical quenching is employed (Figs. 13, 14).

It remains to explore the role of solar wind and coronal mass ejections, which might help in disposing of small-scale helicity from the Sun and thus alleviate catastrophic quenching. An investigation of the effects of Vishniac & Cho fluxes found them to be of secondary importance compared to diffusive

helicity fluxes for αΩ mean field dynamos (Guerrero et al. 2010). Unfortunately the direct numerical simulations have not yet reached the modest Reynolds numbers used in this paper (~10⁴), which are still much lower than the astrophysical dynamos. In order to verify if the equation for dynamical quenching works for the αΩ dynamos in the same way as it does in α² dynamos, we need to embark upon systematic comparisons between DNS with shear and convection and mean field modelling.

Acknowledgments

We are grateful to Eric Blackman for comments that helped to improve the paper considerably. This work was supported in part by the European Research Council under the AstroDyn Research Project No. 227952 and the Swedish Research Council Grant No. 621-2007-4064.

References

All Figures

	Fig. 1 Profiles of radial shear ∂Ω/∂r, α and η as a function of fractional solar radius.
In the text

	Fig. 2 Critical α in terms of a fraction of ηtkf as a function of magnetic Reynolds number Rm for the interface dynamo model of Fig. 1.
In the text

	Fig. 3 Time evolution of the volume-averaged magnetic energy in the domain scaled with the equipartition energy for the case with algebraic quenching and gα = 1 for models with Rm = 1 (diamond), Rm = 20 (solid), Rm = 200 (dashed), Rm = 2 × 103 and Rm = 2 × 105 (triangles).
In the text

	Fig. 4 Same as Fig. 3 but with dynamical quenching and gα = 0 for the Rm indicated in the figure.
In the text

	Fig. 5 Saturated value of the volume-averaged magnetic energy scaled with the equipartition energy as a function of Rm for models with dynamical α quenching (triangles+solid) and algebraic quenching with gα = 1 (squares + dashed) and with gα = Rm (cross + dashed-dotted). Note that all runs included in this figure used α0 = 2αc.
In the text

	Fig. 6 Radial profiles of Aφ and Bφ at two different latitudes (λ) in the saturated phase for Rm = 2 × 103.
In the text

	Fig. 7 Time-latitude plot of αM(0.72   R⊙,θ) for a) Rm = 20 and b) Rm = 200.
In the text

	Fig. 8 a) Evolution of parity (purely dipolar  = −1 and purely quadrupolar  =  +1) for an interface dynamo model with Rm = 2 × 103 and no fluxes. b) A small part in the time-latitude diagram of Bφ indicated by dotted lines in a) where the parity is changing from quadrupolar to dipolar.
In the text

	Fig. 9 Time-latitude plot of the toroidal field a) and c) and αMb) and d) for interface dynamo models with α = 4αc and Rm = 20.
In the text

Fig. 10 a) Time evolution of the volume averaged magnetic energy in the domain scaled with the equipartition energy for models with Rm = 2 × 103 and κ0 = 10-5 (∗ ) as well as with Rm = 2 × 105 and κ0 = 10-5 (solid line) for the interface dynamo model of Sect. 3.1 with dynamical α quenching. The saturation curve for zero fluxes for the model with Rm = 2 × 105 has been shown by the dashed line; b) and c) show time-latitude diagrams for the toroidal field at the depths indicated for models with Rm = 2 × 103 and κ0 = 10-5.

In the text

	Fig. 11 The time-latitude plot of toroidal field a) and c), and αMb) and d) with α = 2αc for Rm = 2 × 105 and κ0 = 0.01. Meridional snapshots of e) sign(Bφ)(| Bφ | /Beq)1/2 and f) αm × 103 for the same case.
In the text

	Fig. 12 Time evolution of the volume averaged magnetic energy for the flux transport dynamo model of Sect. 3.2 for Rm = 2 × 103 with κ0 = 0.3 (dashed); Rm = 2 × 105 with κ0 = 0.3 (solid); Rm = 2 × 103 with κ = 0 (dashed-dotted); Rm = 2 × 105 with κ0 = 0 (diamond+dashed).
In the text

	Fig. 13 Meridional cross-sections showing the distribution of toroidal field and αM for a Babcock-Leighton dynamo without MC and diffusive helicity fluxes in Eq. (11) for a) Rm = 2 × 103 and b) Rm = 2 × 105. The streamlines of the positive and negative poloidal field are shown by solid and dashed lines respectively. Note that the magnetic field has decayed to very small values for Rm = 2 × 105.
In the text

	Fig. 14 Meridional cross-sections showing the distribution of toroidal field and αM for a Babcock-Leighton dynamo with MC and diffusive helicity fluxes for Rm = 2 × 103 at two different epochs. The streamlines of the positive and negative poloidal field are shown by solid and dashed lines respectively.
In the text