A mean-field Babcock-Leighton solar dynamo model with long-term variability

Dynamo models relying on the Babcock-Leighton mechanism are successful in reproducing most of the solar magnetic field dynamical characteristics. However, considering that such models operate only above a lower magnetic field threshold, they do not provide an appropriate magnetic field regeneration process characterizing a self-sustainable dynamo. In this work we consider the existence of an additional \alpha-effect to the Babcock-Leighton scenario in a mean-field axisymmetric kinematic numerical model. Both poloidal field regeneration mechanisms are treated with two different strength-limiting factors. Apart from the solar anti-symmetric parity behavior, the main solar features are reproduced: cyclic polarity reversals, mid-latitudinal equatorward migration of strong toroidal field, poleward migration of polar surface radial fields, and the quadrature phase shift between both. Long-term variability of the solutions exhibits lengthy periods of minimum activity followed by posterior recovery, akin to the observed Maunder Minimum. Based on the analysis of the residual activity during periods of minimum activity, we suggest that these are caused by a predominance of the \alpha-effect over the Babcock-Leighton mechanism in regenerating the poloidal field.


Introduction
The Sun is a magnetic active star, which undergoes periods of high and low magnetic activity approximately each 11 years. Its dynamical behavior imposes important consequences to Figure 1: Representative scheme of the main processes thought to occur during the solar cycle, departing from (a), an initial poloidal field. (b) and (c) represent the generation of the toroidal field by differential rotation -the Ω-effect. (d) and (e) show the effect of cyclonic turbulence on former toroidal fields, creating small-scale secondary poloidal magnetic fields -the α-effect. Averaged, they result in a net electromotive force generating a new large-scale poloidal field (f), closing the first half part of the magnetic cycle with a new poloidal field (g), with opposite polarity than the initial one. (h) represents the beginning of the Babcock-Leighton mechanism: toroidal flux tubes buoyantly rise to the surface forming sunspots, tilted bipolar regions. In (i), the field from the bipolar region diffuses and reconnects with the correspondents of each hemisphere and with the polar fields. The resulting poloidal flux is advected by meridional circulation to the poles (j), generating the final large-scale poloidal field in (g). poloidal field and thus maintaining the solar cycle (Parker, 1955). This latter process is known as the mean-field α-effect, and unlike the Ω-effect, is far from being totally understood, as well as the preferred place for its action. It is also known that the Lorentz force back-reaction of strong magnetic field inhibits turbulence; so the conventional α-effect would not lead to great effectiveness in regenerating the poloidal field in a dynamo within a strong magnetic field regime (Cattaneo and Hughes, 1996). Other mechanisms, like interface dynamos, would overcome this problem, for the αΩ process would occur in the stably stratified layer comprising the tachocline (Parker, 1993).
Alternatively, the Babcock-Leighton mechanism, may have an important role on the poloidal field regeneration. Differently from the mean α-effect, the Babcock-Leighton mechanism relies on the diffusion of the sunspots magnetic field, operating at the solar surface (Babcock, 1961;Leighton, 1969). However, the mechanism behind sunspot formation remains elusive (e.g. Guerrero and Käpylä (2011)), making it difficult to specify the exact nature of the Babcock-Leighton process in a model. In addition, magnetic pumping at the base of the convection zone provides an interesting complement to the Babcock-Leighton scenario (e.g. Guerrero and de Gouveia Dal Pino (2008)).
Axisymmetric numerical models of the solar dynamo are widely used as a tool to investi-gate the relevance of the main effects supposed to govern the solar cycle (Charbonneau, 2010). Many of the solar features are well reproduced by numerical models, but no agreement is achieved considering the particular causes. The variability of the solar cycle, for example, is generally explained either by stochastic forcing (Choudhuri, 1992;Charbonneau and Dikpati, 2000) or by dynamical nonlinearities (Bushby, 2006). Alternatively, simple form of timedelays arising from the spatial decoupling of the α and Ω-effects operation place in the solar convection zone (like in the Babcock-Leighton case), are known to also yield long-term modulation of the solar cycles (Wilmot-Smith et al., 2006;Jouve et al., 2010), even reaching a chaotic behavior . In this paper we use a 2D kinematic solar dynamo model merging concepts of the mean-field theory and the Babcock-Leighton mechanism, accounting for different magnetic field strength-limiting thresholds, in order to achieve solar-like long-term variability.

Model formulation
To access the variability of the magnetic field in the kinematic context of the solar dynamo, in which the flow field is steady and prescribed, the problem reduces to solving the MHD induction equation for the magnetic field B where U is the flow field and η m the magnetic diffusivity. To a first approximation, the largescale magnetic and flow field can be represented as contributions of their large-scale mean and small-scale fluctuating parts, B = B + b ′ and U = U + u ′ , respectively. Upon substitution of these quantities in equation (1), and proper averaging, we get the mean-field induction equation (Moffatt, 1978), given by The term u ′ × b ′ corresponds to a mean electromotive force E arising from the interactions of turbulent motions with the small-scale magnetic field. It can thus be written in terms of a parameterization of the turbulent effects on the mean-magnetic field as E = αB − β∇ × B . Substitution in the mean-induction equation (2) gives in which the averaging brackets have been omitted, and will remain so throughout. The α term represents the turbulent magnetic helicity, due to cyclonic motions oriented by the Coriolis force and η = η m +β is now the effective magnetic diffusivity, covering both magnetic diffusion at the microscopic level and turbulent diffusion, respectively. An additional important issue arises from working within the kinematic context: how should one deal with the feedback of the magnetic field on fluid motions (the Lorentz force)?
Since the Navier-Stokes equation is not solved, this effect ought to be parameterized in the induction equation. This must be done in a way which enforces the saturation of the magnetic field, based on the equipartition of energy between its small-scale magnetic and turbulent kinetic components. Generally, in the mean-field context, the resulting quenching of the magnetic field is heuristically formulated as a decreasing function of its intensity, B eq being the equipartition magnetic field and α 0 a typical value characterizing the α-effect. Studies of magnetoconvection on the solar interior give B eq ∼ 10 4 G (Fan, 2009), and lead to conjectures that a turbulent α-effect would not represent an effective dynamo mechanism (Cattaneo and Hughes, 1996).
Alternatively to the conventional α-effect, Babcock (1961) and Leighton (1969) suggested that the diffusion of the magnetic field of sunspot pairs played a crucial role on the process of poloidal field regeneration. The bipolar sunspot regions are tilted in a way as the magnetic field of the leading spot, nearer the equator, diffuses and reconnects with the field of the leading spot at the other hemisphere, which has nearly always an opposite polarity. At the same time, the magnetic field of the following spot will also decay and connect with the polar magnetic field, which has opposite polarity. This process will at some point annihilate the flux in the polar region, causing a poloidal polarity reversal. In this case the poloidal and toroidal fields regeneration processes are spatially separated (recall Figure 1h to Figure 1j); a mechanism transporting the new polar surface magnetic flux generated by the Babcock-Leighton mechanism to the bottom of the convection zone is necessary. Facing this requirement, initially a meridional circulation flow was assumed in Babcock-Leighton dynamo modeling. As a matter of fact, a poleward meridional flow is indeed observed at the solar surface (Duvall Jr, 1979). Numerical models comprising this additional physical ingredient are termed flux transport dynamos and they are usually successful in reproducing the equatorward tendency of the toroidal field (Küker et al., 2001). In summary, the meridional flow not only acts as a means to transport the polar surface magnetic field down to the base of the convection zone where it is transformed into a toroidal field, but it drags the toroidal field towards the equator as well, yielding a solar-like behavior of the magnetic field. However, the meridional flow profile used in such model is questionable, for it usually consists of a single meridional cell per hemisphere, with a return flow penetrating deep into the convection zone. From full MHD simulations, the resulting pattern of meridional circulation is thought to be far more complex, being divided into smaller cells in each hemisphere (e.g. Brun et al. (2004)).
Simulations of the rising of thin magnetic toroidal flux tubes from the base of the convection zone suggest that for matching the observed tilts, the flux tubes generating sunspots would need to have an initial magnetic field strength ranging from approximately 10 4 to 10 5 G. Stronger toroidal flux tubes tend to reach the solar surface with almost no tilt, while weaker flux tubes arise at too high latitudes (D'silva and Choudhuri, 1993), hence in contradiction with Joy's law. The Babcock-Leighton mechanism should therefore depend upon the initial intensity of the toroidal flux tube at the base of the convection zone, just above the tachocline.
Despite the observations supporting the Babcock-Leighton mechanism, much discussion exists on whether it is the unique responsible for the poloidal field regeneration, an active but not unique component of this part of the cycle, or yet a minor contributor to the whole process. The first possibility is unlikely, as a dynamo based only on the Babcock-Leighton mechanism is not self-sustainable -it needs the formation of sunspots, and therefore operate only above critical toroidal magnetic field buoyancy values. A probable scenario is that an additional αeffect operates on the convection zone, and the resulting magnetic field is the combined product of both poloidal field regeneration processes.
Numerical simulations based on the Babcock-Leighton mechanism have the tendency to produce equatorially symmetric solutions, in opposition to what is observed in the Sun (Chatterjee et al. 2004). Recent results have shown that an α-effect located within a thin layer just above the tachocline is more successful at yielding equatorially antisymmetric solutions (Bonanno et al., 2002;Dikpati and Gilman, 2001). In this study, we focus on this scenario for the location of the α-effect (although other options are possible, see for instance Käpylä et al. (2009).). Residual turbulent motions acting on toroidal flux tubes right above the tachocline, before these reach a critical strength and become buoyant, support such a disposition of the α-effect. This concept is applied in the present model, in addition to the Babcock-Leighton effect at the solar surface.

Mathematical description
Under the assumption of axisymmetry, the magnetic and flow fields are written in terms of their poloidal and toroidal components in spherical coordinates (r, θ, φ) as in which A φ is the poloidal potential and B φ is the toroidal field. The steady flow profile U is given by the meridional circulation u p and the differential rotation Ω. The poloidal-toroidal decomposition of the magnetic field enables a separation of the mean induction equation (3) into two partial differential equations for A φ and B φ , The appearance of the three numbers quantifying the strength of the processes discussed in the introduction, namely and of the magnetic Reynolds number, results from the nondimensionalization of the equations using the solar radius R as the characteristic length scale and the effective magnetic diffusion time R 2 /η s as the characteristic time scale. B tc φ is the toroidal field just above the tachocline; Ω eq , α o , S o and u o are the rotation rate at the equator, the typical magnitudes of the poloidal source terms, and the peak velocity of the meridional flow at the surface.η p andη t are the normalized effective magnetic diffusivities for the poloidal and toroidal components, respectively. Note that the extra α-term that would arise in equation (7) has been neglected, as it is common in the solar case approximation to suppose that C Ω ≫ C α . The flow specifications, u p and Ω are the same as the ones used in the reference work of Dikpati and Charbonneau (1999), and are shown in the left panels of Figure 2. In this paper, the tachocline will comprise the region from the top of the radiative zone r r = 0.6R, to the top of the region with the largest radial angular velocity gradients, r tc = 0.7R.
It is worth noticing that the S term, representing the Babcock-Leighton poloidal regeneration process at the surface, has been added in an ad-hoc manner to equation (7). Unlike the α-term, it is non-local -it depends on the toroidal field at the tachocline. This comes from the latitudinal tilt given by Joy's law being rather dependent on the initial magnetic field strength of the rising flux tube, therefore at the tachocline (D'silva and Choudhuri, 1993).
The α-effect in equation (7) must be suppressed for magnetic field strengths above a limiting value associated with energy equipartition. Therefore, the α-term will be written following equation (4), as where B eq = 10 4 G and f α (r, θ) is a function of spatial coordinates given by where r 1 = 0.675R, r 2 = 0.725R and d 1 = d 2 = 0.01R. This function constrains the αeffect to a thin layer at the base of the convection zone, just above the tachocline, and to mid-latitudes. As mentionned above, other options are possible, but the detailed exploration of these is beyond the scope of the present study. (In passing, we tried the end-member case of a quasi-homogeneous distribution of α, which can produce a solar-like dynamo, but over a limited range of C α , 1 C α 3.) Similar conjectures apply to the Babcock-Leighton S source term in equation (7), with the difference that it operates between lower and upper limiting values, Here B tc φ min = 10 4 G, B tc φ max = 10 5 G and the radial and latitudinal distribution f S (r, θ) is given by where r 3 = 0.95R, r 4 = 1.0R and d 3 = d 4 = 0.01R, restricting the Babcock-Leighton mechanism to the near-surface layers. The radial profiles of the α and S poloidal source terms are shown in the top right panel of Figure 2. The effective diffusivity follows the concept of Chatterjee et al. (2004), who parameterize the suppression of turbulent diffusion by separating the diffusivity into a poloidal component η p and a toroidal component η t .
in which r 5 = 0.7R, r 6 = 0.72R, r 7 = 0.95R and d 5 = d 6 = d 7 = 0.025R. Such separation comes from the fact that the toroidal field, at least in the deeper layers, tends to be much more intense than the poloidal field, and therefore more effective in suppressing turbulent diffusion. η r is the diffusivity near the radiative zone, η cz the diffusivity in the turbulent convection zone associated with the toroidal magnetic field, and η s the diffusivity in the radial surface layers and for the weaker poloidal field. They are set in this study to η r = 5 × 10 8 cm 2 /s, η cz = 1 × 10 10 cm 2 /s and η s = 3 × 10 11 cm 2 /s. The corresponding radial profiles are shown in the bottom right panel of Figure 2.

Numerical simulation
Given the physical ingredients described in the previous section, equations (7) and (8) are solved numerically over a regular grid in an annular meridional plane covering θ ∈ [0, π] and r ∈ [0.6, 1], that is, from slightly below the tachocline up to the solar surface. The boundary and initial conditions are the following: the inner boundary condition is set as to represent the radiative zone as a perfect conductor. An approximation of this condition gives that at r r = 0.6, A φ = B φ = 0 (Chatterjee et al., 2004). The outer boundary condition is that of a vacuum region, requesting the magnetic field to connect with a potential field in the exterior region (Dikpati and Charbonneau, 1999). As for the initial condition, we use a dipolar field permeating the convective envelope. In this case, A φ = sin θ/r 2 for R ≥ r ≥ 0.7R and zero elsewhere, whereas B φ = 0 everywhere. The solution procedure was performed by an adaptation of the Parody code (Dormy, 1997;Dormy et al., 1998;Aubert et al., 2008), based on a pseudo-spectral method. It rests on a spherical harmonic expansion of the angular dependence of the poloidal and toroidal scalars and finite differences in the radial direction. More details on the code, including benchmarks with published numerical solutions (Jouve et al., 2008), are presented in the Appendix. The results presented in the following were obtained with spectral truncation L max = 65, number of radial points N r = 65 and a constant, non-dimensional time stepping size ∆t = 5 × 10 −6 .

Results
On the basis of helioseismic data (Ω eq = 2π × 460.7 nHz), we fix C Ω = 4.7 × 10 4 ; in this case, variations of the free parameters C α , C S and Rm allow for a broad range of solutions. Meridional circulation measured at the solar surface at mid-latitudes displays an average value of 15 m/s. Considering Rm varying around this value, from 318 to 378, the minimum configuration for a proper dynamo solution consists of C α ≃ 2 and C S ≃ 0.5. In order to access the solar representativeness of the solution, some observable aspects should be matched (Charbonneau 2010): 1. Cyclic polarity reversals with approximately 11 years periodicity; 2. Strong deep toroidal fields (∼ 10 4 − 10 5 G) at a 30 o latitudinal belt migrating equatorward; 3. Poleward migration of the polar radial field (∼ 10 − 100 G); 4. Phase lag of π/2 between the deep mid-latitudinal toroidal and surface polar fields; 5. Antisymmetric coupling of the magnetic fields between the hemispheres; 6. Long-term variability of the solar cycle. and Rm = 318. The gray scale map represents the toroidal field at the tachocline, whereas high latitude black and white contours represent the radial field at the solar surface, with a maximum value of 1,700 G.
A useful way to analyze the solar semblance of simulated results is to display the magnetic field in a time-latitude map and compare it with proper synoptic magnetograms and sunspot butterfly diagrams (Hathaway, 2010). A reference solution of the model is displayed in this form in Figure 3, in which the surface polar field contours are superimposed with the grayscale map of the toroidal field just above the tachocline (r tc = 0.7R). This case was chosen because it met most of the obserational requirements listed above. Criterion 1, for example: the periodicity associated with the solar cycle is of 10.95 years. The main dependence of the periodicity resides on the strength of meridional circulation, a well-known characteristic of flux transport dynamos (Dikpati and Charbonneau, 1999). Using a similar model as the present one, Charbonneau et al. (2005) observed the persistent key part played by the meridional circulation in setting the cycle periodicity.
The magnetic field morphology presented in criteria 2 and 3 is also achieved, even though the polar strength of most of the solutions (peaking at 1, 700 G of the reference case in Figure 3) is an order of magnitude higher than the observed one. A phase lag of approximately π/2 between poloidal and toroidal fields cited in criterion 4 is also a general property of flux transport Babcock-Leighton models (Dikpati and Charbonneau, 1999).
On the other hand, the parity requirement 5 is still an issue. Although accounting for an αeffect at a thin layer above the tachocline was thought to help yielding anti-symmetric solutions (Dikpati and Gilman, 2001;Bonanno et al., 2002), the parity coupling was not straightforward in this model case. Actually, there does not seem to exist a clear preferred mode for the solutions: they vary between periods of symmetric, anti-symmetric and out of phase modes. In the reference case, as it is noticeable in Figure 3, the activities within each hemisphere are slightly out of phase. This probably originates from the chaotic nature of the solutions. For higher C S and C α , there is a tendency for the magnetic fields to evolve independently in each hemisphere, with no stable phase lag. Charbonneau et al. (2005) analyzed the general chaotic behavior in a Babcock-Leighton dynamo and ascribed its cause to time-delays connected with the spatial segregation of the toroidal and poloidal field regeneration processes. Long-term variability (recall item 6 above) is a consequence of a chaotic behavior. Figure 4 displays the evolution of the toroidal magnetic field energy at a certain latitude (toroidal magnetic field energy is generally used as a proxy for activity cycle amplitude). We observe frequent short periods of minimum activity with a duration of approximately 3 solar cycles -we typically get 8 of these every 1,000 yr. Moreover, the model also reveals periods of extended minimum activity, reminiscent of the Maunder Minimum, in which the cycle is apparently not fully developed, but persists with a kind of residual activity, lasting for approximately 500 years.
The residual activity episode in Figure 4a suggests that there are two different regimes for the large-scale magnetic field behavior. In the same figure, we show the evolution of the spatial distribution of the poloidal and toroidal magnetic fields in the meridional plane at different epochs, one during regular activity ( Figure 4b) and another during the episode of solar quiescence (Figure 4c). During the normal activity period, the magnetic field is mainly largescale. At the tachocline, a strong toroidal field is created by the shearing effect of differential rotation, from this toroidal field, a poloidal field is created at the solar surface by means of the Babcock-Leighton mechanism. The role of the meridional circulation in setting the timing of the solar cycle is clear: the flow transports the toroidal field at the botton of the convection zone towards the equator, which generates the solar-like orientation of the butterfly diagram (recall Figure 3). On the other hand, during the quiescent period (recall Figure 4c), the magnetic field is small-scale and the solar cycle is confined to the bottom of the convection zone. These differences points to a drastic change in the underlying dynamo mechanism. To investigate this further, we now analyze separatedly the α-effect and the Babcock-Leighton mechanism in two distinct dynamo models. Figure 5 shows the toroidal magnetic energy time series of those two models. The poloidal and toroidal fields during half a solar cycle are also shown in the meridional planes. In the Babcock-Leighton case (Figure 5a), the solar cycle evolution is smooth and tends to display a persistent weak-strong amplitude bundling configuration (therefore with twice the cycle period). This feature, resembling the observed solar pattern known as the Gnevyshev-Ohl rule, is also a consequence of the time delays inherent to the Babcock-Leighton mechanism (Charbonneau et al., 2007). The meridional plots of the magnetic field in Figure 5b show that B is large scale, which matches the overall description of the normal activity period of the reference model, depicted in Figure 4b. In the case of a pure Babcock-Leighton scenario, the toroidal field shows a moderate level of antisymmetry about the equator.
The situation is different in the pure α-effect case, shown in Figures 5c and 5d: the typical amplitude of any given cycle is much lower than in the pure Babcock-Leighton scenario, the cycle period is about twice as long, and there are additional high frequency oscillations superimposed to the solar cycle. Figure 5d shows that in contrast with the situation of a Babcock-Leighton dynamo, the field is mostly small-scale, in agreement with the appearance of higher frequencies in the α-effect model time-series (see the power spectra in 5c).  Aware of the different behaviors concerning the different poloidal field regeneration processes, it is possible that the period of minimum activity involves the preponderance of the tachocline α-effect over the Babcock-Leighton mechanism. In such a case, lower initial poloidal fields would result in the generation of toroidal field below the buoyancy instability limit, leading to few sunspot formations. The long-term solar cycle variability would then be mainly driven by the tachocline α-effect, generating weak toroidal fields during periods of minumum activity, but being able to restart the Babcock-Leighton mechanism when the upper threshold of toroidal field strength is achieved. This situation is reminiscent of observed minimum activity periods, such as the Maunder Minimum. In addition, note that this long-lasting quiescent situation is rather rare in our model: only twice did we observe quiescent periods lasting for more than 500 years in our 20,000 year long integration (short-lived quiescent periods are more frequent, see above). It is also worth mentioning that because of the decoupling of the magnetic field between hemispheres, the minimum activity episode of Figure 4a neither starts nor ends simultaneously in the North and in the South. In this case, the southern hemisphere enters the minimum activity phase approximately 200 years after the northern hemisphere.
Most of the mean-field kinematic solar dynamo simulations able to reproduce the minimum activity periods rely either on the introduction of stochastic forcings (Charbonneau, 2005), or on the somehow arbitrary manipulation of the meridional flow and/or Babcock-Leighton poloidal sources (e.g. Karak (2010)). Here, on the account of the results we presented, we may argue that the α-effect located at the tachocline effectively replaces the stochastic forcing in producing long-lasting phases of minimum activity.

Summary and Conclusion
In view of the not self-sustainable character of a solar dynamo relying only on the Babcock-Leighton mechanism for poloidal field regeneration, we have considered an additional α-effect operating in a thin layer above the tachocline, originating from the turbulent effects on magnetic flux tubes just above a critical buoyancy level. Accounting for different limiting ranges of magnetic field on the operation of each effect, concerning their different natures, we have obtained a dynamo solution reproducing the basic solar magnetic field dynamic features, namely: cyclic reversals with a ∼ 11 year periodicity, equatorward migration of the activity belt, poleward migration of the polar radial field, proper phase lag between both, and long-term variability resembling the solar one.
Appropriate antisymmetric magnetic coupling between the hemispheres remains an issue, even if the location of the α-effect at the tachocline had been suggested as a way to solve the parity problem. In fact, the hemispheres appear to behave in a rather dissociated way, not showing any preferred relaxation mode. The decoupling may originate from the chaotic nature of the solution, making the magnetic field B evolve independently in each hemisphere. Further investigations on the parity topic are needed, possibly relying upon the joint spherical harmonic analysis of the modelled B and that of the B observed at the surface of the Sun (see e.g. Stenflo and Vogel (1986) and DeRosa et al. (2012)).
Our study spontaneously presents a Maunder-like grand minimum which, in comparison with other studies based on similar mean-field kinematic dynamo models, does not require the addition of a stochastic forcing to the right-hand side of the dynamo equations. We conclude by suggesting that grand minima periods could be caused by an intermittent phase of the solar dynamo, during which the sole deep and weak α-effect is responsible for the regeneration of the poloidal field. Why the transition from a Backbock-Leighton dominated regime to an α-effect dominated regime occurs remains a matter of investigation.

Acknowledgements
Sabrina Sanchez would like to express her gratitude to Dr. Oscar Matsuura for his valuable suggestions and fruitful discussions. We thank A. S. Brun for discussions and comments on an earlier version of the manuscript. This work is a result of a M.Sc. project in collaboration between Observatório Nacional and IPGP. We thank both institutions for the available facilities and the Brazilian agency CAPES for funding this research. The contribution of AF and JA is IPGP contribution number NNNN.

Appendix Parody Code -Mean Field Benchmarking
The Parody code used in this work was originally proposed for full 3D MHD dynamo simulations (ACD code, benchmarked in Christensen et al. (2001), see Dormy et al. (1998) and Aubert et al. (2008)). The magnetic field is decomposed according to and the scalar potentials further expanded on a spherical harmonic basis The radial dependence is treated by a second-order finite-differencing scheme. Time integration uses a Crank-Nicholson scheme for diffusion terms and an Adams-Bashforth scheme of order 2 for the nonlinear terms. The actual equations to be solved are the radial curl and radial curl of the curled induction equation (1). Further adaptation to the axisymmetric (M max = 0) mean-field scenario included a change of the magnetic boundary conditon at the inner boundary, the specification of a steady flow in terms of a differential rotation and a meridional circulation, the construction of the proper diffusivity profiles and the addition of the α and S source terms to the right-hand side of the induction equation. The perfectly conducting inner boundary condition implies setting B p = 0 and ∂(rB t )/∂r = 0 at the inner boundary. The code was tested by comparing its predictions with the published reference solutions of a community mean-field benchmark (Jouve et al., 2008). The goal here is to compute critical dynamo numbers and cycle frequencies for different dynamo models, involving either an αΩ scenario (cases A and B, differing only in the prescribed diffusivity) or a Babcock-Leighton scenario (case C). Table 1 displays the values obtained with our code against the reference ones. The butterfly diagrams for the supercritical cases SB and SC (which incorporate an α-quenching) are displayed in Figure 6. Table 1: Comparison of the critical dynamo numbers C crit α,S and frequency of the solar cycle ω within the benchmark cases A, B and C from Jouve et al. (2008). The spatial and temporal resolutions are given in terms of radial points and harmonic degree truncation (N r × L max ) and time-step size ∆t.