Shear-driven and diffusive helicity fluxes in alpha-Omega dynamos

We present nonlinear mean-field alpha-Omega dynamo simulations in spherical geometry with simplified profiles of kinematic alpha effect and shear. We take magnetic helicity evolution into account by solving a dynamical equation for the magnetic alpha effect. This gives a consistent description of the quenching mechanism in mean-field dynamo models. The main goal of this work is to explore the effects of this quenching mechanism in solar-like geometry, and in particular to investigate the role of magnetic helicity fluxes, specifically diffusive and Vishniac-Cho (VC) fluxes, at large magnetic Reynolds numbers (Rm). For models with negative radial shear or positive latitudinal shear, the magnetic alpha effect has predominantly negative (positive) sign in the northern (southern) hemisphere. In the absence of fluxes, we find that the magnetic energy follows an Rm^-1 dependence, as found in previous works. This catastrophic quenching is alleviated in models with diffusive magnetic helicity fluxes resulting in magnetic fields comparable to the equipartition value even for Rm=10^7. On the other hand, models with a shear-driven Vishniac-Cho flux show an increase of the amplitude of the magnetic field with respect to models without fluxes, but only for Rm<10^4. This is mainly a consequence of assuming a vacuum outside the Sun which cannot support a significant VC flux across the boundary. However, in contrast with the diffusive flux, the VC flux modifies the distribution of the magnetic field. In addition, if an ill-determined scaling factor in the expression for the VC flux is large enough, subcritical dynamo action is possible that is driven by the action of shear and the divergence of current helicity flux.


INTRODUCTION
A crucial point in the study of astrophysical dynamos is to understand the mechanism by which they saturate. Nevertheless, a consistent description of this process has rarely been considered in mean-field dynamo (MFD) modeling and only a heuristic description is often used. An important phenomenon happens when the dynamo operates in closed or periodic domains: the turbulent contribution to the dynamo equation, i.e., the α effect, decreases for large values of the magnetic Reynolds number. This process is known as catastrophic quenching and can pose a problem in explaining the generation of magnetic field in late type stars like the Sun or the Galaxy, where Rm could be of the order of 10 9 or 10 15 , respectively.
In the last few years the nature of the catastrophic quenching has been identified as a consequence of magnetic helicity conservation (for a review see Brandenburg & Subramanian 2005a). It has been found that ⋆ E-mail: guerrero@nordita.org (GG) in the nonlinear phase of the dynamo process, conservation of magnetic helicity gives rise to a magnetic α effect (αM) with a sign opposite to the inductive contribution due to the helical motions, i.e., the kinematic α effect. As the production of αM depends on Rm, the final value of the magnetic field should also follow the same dependence. However, real astrophysical bodies are not closed systems, but they have open boundaries that may allow a flux of magnetic helicity. The shedding of magnetic helicity may mitigate the catastrophic α quenching.
These ideas have been tested in direct numerical simulations (DNS) in both local Cartesian and global spherical domains. In the former (Brandenburg 2005;Käpylä, Korpi and Brandenburg 2008) it has been clearly shown that open boundaries (e.g. vertical field boundary conditions) lead to a faster saturation of a large-scale magnetic field compared with cases in closed domains (perfect conductor or triple-periodic boundary conditions). In the latter, it has been found that it is possible to build up large-scale magnetic fields either with forced turbulence (Brandenburg 2005;Mitra et al. 2010b) or with convectively c 2002 RAS driven turbulence (e.g., Brown et al. 2010;Käpylä et al. 2010). These models generally used vertical field boundary conditions.
In flux-transport dynamos (Dikpati & Charbonneau 1999;Guerrero & de Gouveia Dal Pino 2008) as well as in interface dynamos of the solar cycle (e.g. MacGregor ) the quenching mechanism has been considered either through an ad hoc algebraic equation or by phenomenological considerations (Chatterjee, Nandy & Choudhuri 2004), but most of the time the models do not consider the effects of magnetic helicity conservation. An exception is the recent paper by , where these effects have been considered in the context of an interface dynamo.
In general the magnetic helicity depends on time, so it is necessary to solve an additional dynamical equation for the contribution of the small-scale field to the magnetic helicity together with the induction equation for the magnetic field. In the past few years, some effort has already been made to consider this dynamical saturation mechanism in MFD models like in the 1D α 2 dynamo models presented in Brandenburg, Candelaresi & Chatterjee (2009), in axisymmetric models in cylindrical geometry for the galactic αΩ dynamo (Shukurov et al. 2006), and also in models with spherical geometry for an α 2 dynamo (Brandenburg et al. 2007). The role of various kinds of magnetic helicity fluxes have been explored in several papers (Brandenburg, Candelaresi & Chatterjee 2009;Zhang et al. 2006;Shukurov et al. 2006).
Our ultimate goal is to develop a self-consistent MFD model of the solar dynamo, with observed velocity profiles and turbulent dynamo coefficients computed from the DNS. This is a task that requires intensive efforts. Hence we shall proceed step by step, starting with simple models and then including more realistic physics on the way. In this work we will study the effects of magnetic helicity conservation in simplified αΩ dynamo models for a considerable number of cases. More importantly, we shall perform our calculations in spherical geometry, which is appropriate for describing stellar dynamos, with suitable boundary conditions, and considering shear profiles which are a simplified version of the observed solar differential rotation. We shall also explore how magnetic helicity fluxes affect the properties of the solution. Two classes of fluxes are considered in this paper: a diffusive flux and a shear-driven or Vishniac-Cho (hereafter VC) flux (Vishniac & Cho 2001). We consider models with either radial or latitudinal shear. The effects of meridional circulation will be investigated in detail in a companion paper ). This paper is organized as follows: in Section 2 we describe the basic mathematical formalism of the αΩ dynamo, give the formulation of the equation for αM and also justify the fluxes included. In Section 3 we describe the numerical method and then, we present our results in Section 4 starting from a dynamo model with algebraic quenching to models with dynamical α quenching and different fluxes. Finally, we provide a summary of this work in Section 5.

THE αΩ DYNAMO MODEL
In mean-field dynamo theory, the evolution of the magnetic field is described by the mean-field induction equation, where B and U represent the mean magnetic and velocity fields, respectively, ηm is the molecular diffusivity, E = αB − ηtµ0J is the mean electromotive force obtained using a closure theory like the first order smoothing approximation, where E gives the contribution of the small-scale components on the large-scale field, α is the non-diffusive contribution of the turbulence, ηt is the turbulent magnetic diffusivity, J = ∇ × B/µ0 is the mean current density, and µ0 is the vacuum permeability.
In spherical coordinates and under the assumption of axisymmetry, it is possible to split the magnetic and the velocity fields into their azimuthal and poloidal components, B = Bê φ + ∇ × (Aê φ ) and U = r sin θΩê φ + up, respectively. For the sake of simplicity we shall not consider the meridional component of the flow, i.e. up = 0. Then, the toroidal and poloidal components of equation (1) may be written as where D 2 = ∇ 2 − s −2 is the diffusion operator, η = ηm + ηt, s = r sin θ is the distance from the axis, and Bp = ∇×(Aê φ ) is the poloidal field. The two source terms in equations (2) and (3), sBp ·∇Ω and αB, express the inductive effects of shear and turbulence, respectively. The relative importance of these two effects may be quantified through the non-dimensional dynamo numbers: CΩ = ∆ΩL 2 /ηt and Cα = α0L/ηt, where ∆Ω is the angular velocity different between top and bottom of the domain. Note that equations (2) and (3) are valid only in the limit CΩ ≫ Cα, known as αΩ dynamo.
The inductive effects of the shear may be understood as the stretching of the magnetic field lines due to the change in the angular velocity between two adjacent points. On the other hand, the kinematic α-effect is the consequence of helical motions of the plasma which produce screw-like motions in the rising blobs of the magnetic field. Using the first order smoothing approximation it may be expressed as: where, τ is the correlation time of the turbulent motions and ω = ∇ × u is the small-scale vorticity. The saturation value of the magnetic field may be obtained by multiplying αK by the quenching function fq = 1 + B 2 /B 2 eq −1 , which saturates the exponential growth of the magnetic field at values close to the equipartition field strength given by Beq = (µ0ρu 2 ). This form of algebraic quenching was introduced heuristically (see, e.g. Stix 1972) and has been often used as the standard quenching mechanism in many dynamo simulations. However, it does not give information about the back reaction process and is independent of any parameter of the system like the magnetic Reynolds number. A consistent description of the quenching mechanism will be presented in the following section.

Dynamical α effect
Recently, it has been demonstrated that when the amplitude of the magnetic field reaches values near the equipartition, the α-effect is modified by a magnetic contribution, the so called magnetic α effect, denoted by αM. It is usually the case that αM has a sign opposite to αK resulting thus in the saturation of the magnetic field. Pouquet, Frisch & Léorat (1976) have shown that αM is proportional to the smallscale current helicity of the system, hence it is possible to write α as a sum of two contributions, one from the fluid turbulence and other from the magnetic field, as follows: where ρ is the mean density of the medium, assumed here as a constant, and j = ∇ × b/µ0 is the current density of the fluctuating field. The mathematical expression that describes the evolution of αM may be obtained by taking into account the magnetic helicity evolution (Blackman & Brandenburg 2002), which leads to: where k f = 2π/(L − rc) with rc = 0.7L0 is a suitable choice for the wave number of the forcing scale, the magnetic Reynolds number RM = ηt/ηm and F α is the flux of the magnetic α effect related to the flux of the small-scale magnetic helicity, F f through: According to previous authors αM has a finite value in the interior of the domain in absence of fluxes (F α = 0), and its sign is usually opposite to the sign of αK in such a way that the final amplitude of the total α-effect decreases, and so does the final value of the magnetic energy.

Magnetic helicity fluxes
Recently it has been pointed out that the catastrophic quenching could be alleviated by allowing the flux of smallscale magnetic (or current) helicity out of the domain, so that the total magnetic helicity inside need not be conserved any longer. Alternately, we may introduce those fluxes in the equation for αM; see equation (7). Several candidates have been proposed for the helicity fluxes in the past (Kleeorin & Rogachevskii 1999;Vishniac & Cho 2001;Subramanian & Brandenburg 2004). Amongst them are the flux of magnetic helicity across the iso-rotation contours, advective and diffusive fluxes and also the explicit removal of magnetic helicity in processes like coronal mass ejections or galactic fountain flows, for the case of the galactic dynamo. From the mathematical point of view, the nature of the flux terms in the equation for αM has not been demonstrated with sufficient rigor. However, several DNS have pointed to its existence.
Firstly, the shearing box convection simulations of Käpylä, Korpi and Brandenburg (2008) showed that in the presence of open boundaries, the large-scale magnetic field grows on temporal scales much shorter than the dissipative time scale. They concluded from this that open boundaries may allow the magnetic helicity to escape out of the system. These experiments seem to be compatible with the flux proposed by Vishniac & Cho (2001), whose functional form may be expressed as (see Subramanian & Brandenburg 2004;Brandenburg & Subramanian 2005b, for further details): where S lk = 1 2 (U l,k + U k,l ) is the mean rate of strain tensor and CVC is a non-dimensional scaling factor. As we assume up = 0, this flux has the following three components: with S φr = S rφ = r sin θ(∂Ω/∂r)/2 and S θφ = S φθ = sin θ(∂Ω/∂θ)/2. Secondly, Mitra et al. (2010a) performed α 2 dynamo simulations driven by forced turbulence in a box with an equator. They found that the diffusive flux of αM across the equator can be fitted to a Fickian diffusion law given by, They also computed the numerical value of this diffusion coefficient, and found it to be of the order of turbulent diffusion coefficient. They also found that the time averaged flux is gauge independent. Both results were later corroborated by simulations without equator, but with a decline of kinetic helicity toward the boundaries (Hubbard & Brandenburg 2010). Additionally, magnetic helicity may be advected by the mean velocity with a flux given by F ad = αMU , or it may be expelled from the solar interior by coronal mass ejections (CMEs) or by the solar wind. This flux, FCME, may account for ∼ 10% of the total helicity generated by the solar differential rotation, as estimated by Berger & Ruzmaikin (2000). It can be modeled by artificially removing a small amount of αM every τ time (Brandenburg, Candelaresi & Chatterjee 2009), or also by a radial velocity field that mimics the solar wind.
The total flux of magnetic helicity may be written as the sum of these contributions, Since in this dynamo model we do not include any component of the velocity field other than the differential rotation, in this study we will consider only the first two terms on the rhs of equation (13).

THE MODEL
We solve equation (2), (3) and (6) for A, B and αM in the meridional plane in the range 0.6L r L and 0 θ π. We consider two different layers inside the spherical shell. In the inner one the dynamo production terms are zero and go smoothly to a finite value in the external layer. The magnetic diffusivity changes from a molecular to a turbulent value from the bottom to the top of the domain. This is achieved by considering error function profiles for the magnetic diffusivity, the differential rotation, and the kinetic α effect, respectively (see Fig. 1): where Θ(r, r1,2, w) = 1 2 [1 + erf {(r − r1,2)/w1}], with r1 = 0.7L0, r2 = 0.72L0 and w1 = 0.025L0. We fix CΩ = −10 4 and vary Cα.
The boundary conditions are chosen as follows: at the poles, θ = 0, π, we impose A = B = 0; at the base of the domain, we impose a perfect conductor boundary condition, i.e. A = ∂(rB)/∂r = 0. Unless noted otherwise, we use at the top a vacuum condition by coupling the magnetic field inside with an external potential field, i.e., (∇ 2 − s −2 )A = 0. A good description of the numerical implementation of this boundary condition may be found in Dikpati & Choudhuri (1994).
The equations for A and B are solved using a secondorder Lax-Wendroff scheme for the first derivatives, and centered finite differences for the second-order derivatives. The temporal evolution is computed by using a modified version of the ADI method of Peaceman & Rachford (1955) as explained in Dikpati & Charbonneau (1999). This numerical scheme has been used previously in several works on the flux-transport dynamo and the results were found to be in good agreement with those using other numerical techniques (Guerrero & de Gouveia Dal Pino 2007Guerrero, Dikpati & de Gouveia Dal Pino 2009).
In the absence of magnetic helicity fluxes, equation (6) for αM corresponds to an initial value problem that can be computed explicitly. However, as we are going to include a diffusive flux, we use for αM the same numerical technique used for A and B. All the source terms on the right hand side of equation (6) are computed explicitly. We have tested the convergence of the solution for 64 2 , 128 2 , and 256 2 grid points. For cases with small Rm, there are no significant differences between different resolutions, but for high Rm, 64 2 grid points is insufficient to properly resolve the sharp diffusivity gradient. A resolution of 128 2 grid points is a good compromise between accuracy and speed.

αΩ dynamos with algebraic quenching
In order to characterize our αΩ dynamo model we start by exploring the properties of the system when the saturation is controlled by algebraic quenching with fq = (1+B 2 /B 2 eq ) −1 . We found that, with the profiles given by equations (14)-(16), Fig. 1, the critical dynamo number is around 2 × 10 4 (i.e., C C α = 1.975). The solution for the model is a dynamo wave traveling towards the equator since it obeys the Parker-Yoshimura sign rule (see Fig. 2). In this case, the maximum amplitude of the magnetic field depends only on the dynamo number of the system, CαCΩ, as can be seen in the bifurcation diagram in Fig. 3. The quenching formula is here independent of Rm, so the saturation amplitude is also independent on Rm.

αΩ dynamos with dynamical quenching
In this section we consider dynamo saturation through the dynamical equation for αM described in Section 2.1. In this models we distinguish three different stages in the time evolution of the magnetic field: a growing phase, a saturation phase and a final relaxation stage (see panels a, b and c of Fig. 5). The magnetic field is amplified from its initial value, 5 × 10 −4 Beq, following an exponential growth. From the ear-  liest stages of the evolution we notice the growth of αM with values that are predominantly negative in the northern hemisphere and positive in the southern hemisphere. The latitudinal distribution of αM is fairly uniform in the active dynamo region, spanning from the equator to ∼ 60 • latitude. The radial distribution exhibits two narrow layers where the sign of αM is opposite to the dominant one developing at each hemisphere. These are located at the base of the dynamo region (r ∼ 0.7L0) and at a thin layer near to the surface (r > 0.95). In the equation for the magnetic α effect, equation (6), the production term is proportional to E · B = αB 2 − ηtµ0J · B. The first component of this term has the same sign as αK, which in general is positive in the northern and negative in the southern part of the domain. The minus sign in front of the right hand side of equation (6) defines then the sign of αM. However, at the base and at the top of the dynamo region, αK → 0 and B → 0, respec-tively. The term ηtJ · B is the only source of αM and leads to the formation of these two thin layers. The space-time evolution of αM depends on the value of the magnetic Reynolds number. For small Rm, the decay term in equation 6 (i.e. the second term in the parenthesis) becomes important, so that there is a competition between the production and decay terms resulting in an oscillatory behavior in the amplitude of the magnetic α effect, as is indicated by the vertical bars in the middle panel of Fig. 7. The period of these oscillations is the half the period of the magnetic cycle. With increasing Rm, the amplitude of the oscillations decreases such that for Rm 10 3 , αM is almost steady.
The morphology of the magnetic field corresponds to a multi-lobed pattern of alternating polarity (left panels of Fig. 5). These lobes are radially distributed in the whole dynamo region with maximum amplitude at the base of this layer. The poloidal magnetic field follows a similar pattern with lines that are open at the top of the domain due to the potential field boundary condition. There is a phase shift between toroidal and poloidal components which we have estimated to be ∼ 0.4π. The model preserves the initial dipolar parity during the entire evolution.
The evolution of αM traces the growth of the magnetic field, but its final value depends on the magnetic Reynolds number. For small Rm, after saturation, αM reaches a steady state, but for large Rm, its relaxation is modulated by overdamped oscillations. The relaxation time is proportional to Rm, which means that for Rm ≫ 1 the simulation must run for many diffusion times. The differences in the relaxation time observed for αM reflects the evolution of the magnetic field, as is shown in Fig. 4.
We observe that the rms value of the magnetic field remains steady during the saturation phase for Rm < 10 2 . For 10 2 < Rm < 10 3 , a bump appears in the curve of magnetic field evolution, followed by the relaxation to a steady value, whereas for Rm > 10 3 , the magnetic energy shows over-damped relaxations with a final energy proportional to Rm −1 as has been previously reported (Brandenburg et al. 2007). These oscillations in the time evolution plot of the averaged magnetic field have been reported in mean field dynamo simulations including the dynamical α-effect (Brandenburg & Subramanian 2005b).
Not many DNS of αΩ dynamo exist so far in the literature with Rm 100 in order to compare with our results. However, in the local αΩ dynamo simulations of Käpylä, Korpi and Brandenburg (2008), a rapid decay of the magnetic field seems to occur after the initial saturation for moderate values of Rm. This decay forms a bump in the curve of the averaged magnetic field (see their Fig. 14), similar to the bump that we obtain for 10 2 < Rm < 10 3 .
For reasons of clarity in the Fig. 4 we do not show the entire time evolution of each simulation with Rm > 10 3 . The total evolution time as well as the final value of the magnetic field of each simulation are shown in the Table 1. For magnetic Reynolds numbers above 2 × 10 4 , the initial kinematic phase is followed by a decay phase during which the total α effect goes through subcritical values and then the dynamo fails to start again.
In Fig. 5 we present the meridional distribution of the magnetic field (left panel), αM (middle panel) and the total α (right panel), in normalized units, for the three dif- Note that for Rm > 10 3 , we have allowed the simulations to evolve more than 4 diffusion times, as indicated in Table 1. ferent stages of evolution corresponding to the early kinematic phase, the late kinematic phase and the saturated phase. These snapshots correspond to the simulation with Rm = 10 3 (Run Rm1e3 in Table 1). The multi-lobed pattern of the toroidal field represented with filled contours remains unchanged during the evolution even though its amplitude increases. The same occurs for the poloidal component, shown by continuous and dashed streamlines for positive and negative values, respectively.
The magnetic α effect (middle panels) is formed first at latitudes between ±30 • and then it amplifies and expands to latitudes up to ∼ ±60 • . This makes the total α effect, initially similar to αK ( Fig. 1 and top panel of Fig. 5a), smaller at lower latitudes in the central area of the dynamo region. At the bottom and at the top of the domain αM and αK have the same sign making the total α larger. However, the global effect is a decrease of the dynamo efficiency.

Diffusive flux for αM
In this section we consider a Fickian diffusion term in equation (12) for αM. We consider a diffusion coefficient varying from 5 × 10 −3 ηt to 10 ηt in the dynamo region and with κα = ηm in the bottom layer. In these cases, the initial evolution of αM is similar to the cases presented in the previous section: negative (positive) values for αM in the northern (southern) hemisphere, with narrow regions of opposite values nearby the regions where αK = 0 or B = 0. However, at the later stages, αM is much more diffuse in the entire domain and has only one sign in each hemisphere. This is the result of cancellation of αM with opposite signs occurring in each hemisphere due to radial diffusion. Contrary to the cases without fluxes, we now obtain finite values of Bsat for large values of Rm, as can be seen in Fig. 6. All the cases depicted in this figure correspond to κα = 0.005ηt. We notice that the final value of the magnetic field still remains small compared to the equipartition ( 0.1Beq), but it is clear that even this very modest diffusion prevents the α effect from being catastrophically quenched. This is also evident from the top panel of Fig. 7, where we plot the final strength of B as a function of Rm, for the cases with and without dissipative flux. In the middle and bottom panels of the Fig.  7 we compare the behavior of the normalized αM, at a given point inside the dynamo region, and also the time period, T , of the dynamo for models with and without fluxes. In both panels it is clear that for Rm above ∼ 10 3 , αM and T reach a saturated value.
Besides its dependence on Rm, the evolution of αM depends also on κα. For models with κα ≪ ηt, the evolution of αM relies on Rm, but for κα 0.1ηt, the dissipation time of αM becomes comparable to, or even shorter, than the period of the dynamo cycle. This results in αM becoming oscillatory, as shown in the bottom panel of Fig. 8. The amplitude and the period of these oscillations depend on the value of κα.
In the top panel of Fig. 8 we show the final value of the averaged mean magnetic field as a function of κα. We observe that for κα in the range (0.1-1) ηt, the value of Brms remains between 20% and 60% of the equipartition, a value similar to the one obtained in the simulations using algebraic α quenching (Section 4.1, Fig. 3). For κα > ηt, superequipartition values of the magnetic field may be reached. This is because larger values of κα result in oscillations of αM with larger amplitude, such αM may locally change its sign, increasing the value of the total α in each hemisphere and thereby enhancing the dynamo action. Such high values of the diffusion of the magnetic helicity are unlikely in nature.

The Vishniac-Cho flux
Our next step is to explore the magnetic helicity flux proposed by Vishniac & Cho (2001) in the form given by equation (8). For the moment we set κα = 0. In a previous study on the effects of the VC flux in a MFD model in Cartesian coordinates, Brandenburg & Subramanian (2005b) found that there exist a critical value for the parameter CVC above which there is a runaway growth of the magnetic field that can only be stopped using an additional algebraic quenching similar to the one used in Section 4.1. They found that this critical value, CVC * , diminishes with increasing the amount of shear. Since we have used a strong shear (CΩ = −10 4 ) we use nominal values of CVC = 10 −3 , but without any algebraic quenching.
The term ∇·F VC develops a multi-lobed pattern which travels in the same direction as the dynamo wave, this confirms that the VC flux follows the lines of iso-rotation. From equation (8), we see that the VC flux is proportional to the magnetic energy density. In the present case, with CΩ ≫ Cα, the spatial distribution of ∇ · FVC/B 2 eq is dominated by the terms involving B 2 φ in equations. 9-11 (this may be inferred from the left hand panels of Fig. 10a). This results in a new distribution of αM, with concentrated regions of positive (negative) sign at low latitudes in the northern (southern) hemisphere, and a broad region of negative (positive) sign in latitudes between 20 • and 60 • latitude (see middle panels of Fig. 10). Surprisingly we find that the general effect of this flux is to decrease the final amplitude of the magnetic field with respect to the case without any fluxes as can be seen in Fig. 11. Note that we have until now used only the potential field boundary condition for the poloidal field.  When we consider both diffusive as well as VC fluxes, with κα = 0.1ηt and CVC = 10 −3 , we obtain a magnetic field of slightly larger amplitude compared to the case with only the diffusive flux (compare the value of Brms in Runs DRm1e3 and VCD in Table 1). However we may say from the butterfly diagram of Fig. 12 that the toroidal magnetic field appears to be more concentrated at lower latitudes, where the sign of αM is same as that of αK.
With negative values of CVC, it was found that the resulting profile of αM is only weakly modified from cases without fluxes, though its value is reduced marginally such that the final amplitude of Brms is slightly larger. But even this contribution does not help in alleviating catastrophic quenching in models with large Rm (see Fig. 11).
Since VC fluxes transport helicity along lines of constant shear, it may be expected that they are more important in models with latitudinal shear, since in this case the magnetic helicity flux can travel either towards the bottom or the top boundaries, from where magnetic helicity can be expelled. For testing this possibility, we turn off the radial shear profile and consider a purely latitudinal solar-like differential rotation: where Ωeq/2π = 460.7 nHz is the angular velocity at the equator, and Ωs(θ) = Ωeq + a2 cos 2 θ + a4 cos 4 θ gives the latitudinal profile, with a2/2π = −62.9 nHz and a4/2π = −67.13 nHz.
In order for the dynamo to be slightly supercritical, as in the previous cases, we consider CΩ = 5 × 10 4 . This dynamo solution corresponds to a dynamo wave produced at mid latitudes (∼ 45 • ) that travels upwards (since CΩ now is positive). As in the previous cases with radial shear, the distribution of ∇ · F VC/B 2 eq is similar to that of the divergence of magnetic energy density (left hand panels of Fig. 10 b,c and d). If no fluxes are considered, the final amplitude of the mean magnetic field is ∼ 0.03% of the equipartition value. In presence of VC fluxes, starting with CVC = 10 −3 for a model with Rm = 10 3 , we notice that the final magnetic field is twice as large as in the case with CVC = 0.
Our model becomes numerically unstable beyond CVC = 10 −2 due to appearance of concentrated regions of strong αM. When VC and diffusive fluxes are considered simultaneously, with CVC = 10 −3 and κα = 0.1ηt, the relaxed value of Brms is only slightly below the value reached at the end of the kinematic phase (Fig. 11b). In this case αM spreads out in the convection zone, as shown in Fig. 10c, indicating that the effects of the VC flux are not important when compared with the diffusive flux.
We repeated the calculation by considering the vertical field (VF) boundary condition, ∂(rB θ )/∂θ = 0, for the top boundary, instead of the potential field (PF) condition used throughout the rest of this work. Furthermore, in the models with VF conditions the presence of the VC flux leads to an increase of Bsat by a factor of ∼ 2 compared to the case without VC flux (see Fig. 11c). It may be noted that αM shows regions of both positive and negative signs in each hemisphere (see Fig. 10d). Thus, the total α effect is increased locally to values well above the kinematic one. This implies that in the region around ±45 • the dynamo action is driven by the magnetic α effect. A similar secondary dynamo is found to be working for a different distribution of shear and αK . As with PF boundary condition, large values of CVC result in a numerical instability of the magnetic field in the simulation with VF. Figure 11. Time evolution of the averaged mean magnetic field for different values of C VC : a) Radial shear, b) latitudinal shear with potential field boundary conditions and c) latitudinal shear with vertical field boundary conditions. The width of the different bands reflects the range over which the magnetic field varies during one cycle. Note that the cycle period is short compared with the resistive time scale on which the magnetic field reaches its final saturation. If not indicated, in all models Rm = 10 3 . The two dashed lines in the panel a) corresponds to Cvc = −0.002 for Rm = 10 3 and Rm = 10 4 .
The main result of this section is that the VC flux does not alleviate catastrophic quenching of the dynamo for large values of Rm (see the dashed lines in Fig. 11 a and c). The reason for this may be related to the fact that the radial flux has components that are either proportional to B θ or to B φ (equation 9). As B φ vanishes on the top boundary, and B θ is small, the VC flux is not able to dispose of αM across

CONCLUSIONS
We have developed αΩ dynamo models in spherical geometry with relatively simple profiles of αK and shear (∂Ω/∂r and ∂Ω/∂θ). We choose potential field (also vertical field in some cases) and perfect conductor boundary conditions for the top and bottom boundaries, respectively. We estimate the critical dynamo number by fixing CΩ = −10 4 and varying Cα while using algebraic quenching. Using a dynamo number, CΩCα, that is slightly supercritical, we solve the induction equations for B and A together with an equation for the dynamical evolution of the magnetic α effect or αM. We find that for positive (negative) values of Cα in the northern (southern) hemisphere, αM is mainly negative (positive), with narrow fractions of opposite sign in regions where αK or B are equal to zero.
We find that the kinematic phase is independent of Rm. However for Rm > 10 2 there exists a phase of relaxation post saturation in which the averaged magnetic field oscillates about a certain mean. The larger the Rm, the more pronounced are the damped oscillations and the longer is the relaxation time (Fig. 4). The final value of the magnetic energy obeys a R −1 m dependency (R −0.5 m for magnetic field, Fig. 7), which is in agreement with earlier work (Brandenburg & Subramanian 2005b;Brandenburg, Candelaresi & Chatterjee 2009).
We argue that including equation (6) in MFD models is appropriate for describing the quenching of the magnetic field in the dynamo process. Since we observe large-scale magnetic fields at high magnetic Reynolds numbers in astrophysical objects, there must exist a mechanism to prevent the magnetic field from catastrophic quenching.
We have studied the role that diffusive and VC fluxes may play in this sense. Their contribution may be summarized as follows: (i) In the presence of diffusive fluxes, αM has only one sign in each hemisphere (negative in the northern hemisphere and positive in southern) and is evenly distributed across the dynamo region (Fig. 9).
(ii) For Rm < 10 2 the mean values of αM are similar to models without diffusive fluxes, whereas for Rm 10 2 , αM has smaller values that seem to be independent of Rm (see Fig. 7, middle).
(iii) Even a very low diffusion coefficient, e.g. κα = 0.001ηt, causes Brms to depart from the R −0.5 m tendency and converge to a constant value which is then around 5% of the equipartition value for large values of Rm, but below the value of 10 7 used in this study (dashed line in Fig. 7, top).
(iv) Larger values of κα result in larger final field strengths.
(v) In models with only radial shear the Vishniac-Cho flux contributes to αM with a component that travels in the same direction as the dynamo wave. This produces a different radial and latitudinal distribution of the magnetic α effect that also affects the distribution of the magnetic fields. However, it does not help in alleviating the quenching at high Rm. On the contrary, the larger the coefficient CVC, the smaller is the resultant magnetic field.
(vi) In models with only latitudinal shear the VC flux travels radially outward but it remains concentrated at the center of the dynamo region. In a given hemisphere the resultant distribution of αM has both positive and negative signs. The part of αM that has the same sign as αK enhances dynamo action. This effect is more evident in models with vertical field boundary conditions .
(vii) In models with vacuum and vertical field boundary conditions and Rm = 10 3 , the VC flux increases the final value of the magnetic field by a factor of two compared to the case without any fluxes.
(viii) The magnetic field in models with Rm 10 4 and with non-zero VC flux decays after the kinematic phase since the total α effect becomes subcritical (see dashed lines in Fig. 11 a and c).
(ix) Larger values of CVC produce narrow bands of αM which drives intense dynamo action in these regions. This positive feedback between the magnetic field and αM causes the simulation to become numerically unstable in the absence of any other quenching effect.
From the above results it is clear that diffusive fluxes are much more important in alleviating catastrophic quenching when compared to the Vishniac & Cho fluxes (in the form of equation 8) for a large range of Rm. This is somehow intriguing since it is known from DNS that shear in do-mains with open boundaries does indeed help in alleviating the catastrophic quenching. It may be understood as a result of the large value of CΩ compared with Cα and also to the top boundary condition for the azimuthal magnetic field (Brandenburg 2005;Käpylä, Korpi and Brandenburg 2008).
The results presented above indicate that considerable work is still necessary in order to understand the role of larger-scale shear in transporting and shedding small-scale magnetic helicity from the domain.
In snapshots of the meridional plane as well as in butterfly diagrams we notice that the diffusive fluxes do not significantly modify the morphology and the distribution of the magnetic field when compared with cases without fluxes or even with simulations with algebraic α quenching. On the other hand, for models with VC flux the distribution of αM becomes different and so does the magnetic field. This is clear from the butterfly diagram shown in Fig. 12b, which exhibits a magnetic field confined to equatorial latitudes reminiscent of the observed butterfly diagram of the solar cycle. Even though this result corresponds to a simplified model, it illustrates the importance of considering the dynamical α quenching mechanism for modeling the solar dynamo. Similar changes in the distribution of αM and B are expected to happen when advection terms are included in the governing equations.
In the simulations presented here, Ω and α effects are present in the same layers. An interesting question is whether the quenching of the dynamo is catastrophic when both layers are segregated, as in the Parker's interface dynamo or the flux-transport dynamo models. We address this question in detail in two companion papers .
We should notice that the back reaction of the magnetic field affects not only the α effect, but also the other dynamo coefficients, including the turbulent diffusivity. Contrary to quenching of α, the quenching of ηt may be considered through an algebraic quenching function (see e.g. Yousef, Brandenburg & Rüdiger 2003;Käpylä & Brandenburg 2009). Guerrero, Dikpati & de Gouveia Dal Pino (2009) have shown that in a flux-transport model these effects could affect properties of the models such as the final magnetic field strength and its distribution in radius and latitude. We leave the study of models with simultaneous dynamical α and η quenchings for a future paper. Solar-like profiles of differential rotation and meridional circulation along with dynamical α quenching will also be considered in a forthcoming paper.

ACKNOWLEDGMENTS
This work started during the NORDITA program solar and stellar dynamos and cycles and is supported by the European Research Council under the AstroDyn research project 227952.