Free Access
Issue
A&A
Volume 525, January 2011
Article Number A5
Number of page(s) 9
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201015073
Published online 26 November 2010

© 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 , an expression introduced by Jepps (1975), or by , where α0 is the unquenched value and Rm is the magnetic Reynolds number, is the mean magnetic field and Beq is the equipartition magnetic field. The latter expression has been discussed since the early work of Vainshtein & Cattaneo (1992). The Rm 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 Rm 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 (Bsat) with increasing magnetic Reynolds number (Rm). 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 α2 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 ≤ Rm ≤ 2 × 105. 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 < Rm < 2 × 105, 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) (1)where urms is the rms velocity of the turbulent eddies, kf 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 kf = 10k1, where k1 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 (2)where re = 0.73   R, and de = 0.025   R. Thus the magnetic diffusivity profile is given by (3)The turbulent kinetic α effect considered in this model has the form (4)where gα is a non-dimensional coefficient equal to 1 or Rm depending on the assumed form of algebraic quenching in the models and 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 α0 = τϵfωrmsurms/3, where ωrms is the rms vorticity of the turbulence and τ ~ (kfurms)-1 is the eddy correlation time scale. The prefactor ϵf, usually of the order of 0.1 or less, is used because (u·ω)rms < urmsωrms. The case ϵf = 1 means that the flow is maximally helical. These approximations give us an estimate of α0 in terms of eddy diffusivity ηt and forcing scale kf as We would consider ϵf instead of α0 as a free parameter in the model apart from ηt. Even though the helical turbulence pervades almost the entire convection zone, here we take ra = 0.77   R and da = 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 (5)where Ω0 = 14 nHz, rw = 0.68   R and dw = 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 Tcyl of the oscillatory dynamos to a reasonably small fraction of the diffusion time tdiff.

thumbnail 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 (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 (7)(8)where ζ = R/r−1, rb = 0.71R, ζb = R/rb−1, and . In the equations above, the flow is poleward at the surface with a maximum amplitude of u0 = 20 m s-1, and the equatorward return flow occurs at rb 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 rw = 0.7   R. Finally, the turbulent diffusivity has the same profile as in Eq. (3), but with ηt = 2 × 1011 cm s-1 and re = 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, (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 using the relation (10)As the equation for 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 hf 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, (11)where and 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 hf 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 (). Moreover, and are in principle gauge-dependent. However, if there is scale separation, can be written as a density of linkages and is therefore gauge-independent (Subramanian & Brandenburg 2006), while is not. This was demonstrated in Hubbard & Brandenburg (2010).

The flux term, , 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, , 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 (12)where κ = κ0η(r) is the diffusion coefficient for αM. Let us define the diffusion time in the model as . The decay time in Eq. (11) is then . 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 Rm 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 × 105 cm2 s-1 to 2 × 1010 cm2 s-1 while keeping ηt fixed at 4 × 1010 cm2 s-1.

thumbnail 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 Rm, it is important to first calculate the critical value of α0, 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 Rm are shown in Fig. 2. From this figure we observe that the model is most efficient for Rm ~ 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 α0 = 2αc, corresponding to the Rm for each model. The cycle period of these dynamo models (Tcyl) slightly increases with Rm and is found to be of the order of 2tdiff, for the investigated range of Rm. 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 Rm within the errors in the numerical determination of αc.

thumbnail 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).

thumbnail Fig. 4

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

thumbnail 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 Rm dependence, which is reminiscent of catastrophic quenching in all astrophysical dynamos, can be easily discerned from this figure. For high values of Rm, 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 (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α = Rm, as well as models with dynamical quenching (Eq. (11)). Both algebraically quenched models with gα = Rm and those with dynamical α quenching exhibit a monotonic decrease of the saturation magnetic energy with Rm.

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 Bsat on Rm is reproduced in models where the two source regions are further spatially separated by setting ra = 0.87   R instead of 0.77   R.

thumbnail 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 Tcyl. For low Rm, tα ≤ Tcyl, so that the time evolution of is flat in the saturation regime – similar to the algebraic quenching case (see Fig. 3). On the other hand, for Rm = 2 × 103, the decay time tα ≫ Tcyl 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 Rm = 2 × 105, the underdamped oscillation has a long period (tα) so that the model has to be run for more than 500 tdiff 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 tdiff. Hence, the determination of saturation magnetic energy may be inaccurate.

thumbnail 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 Rm = 20 (upper panel), αM shows strong oscillations (see the modulation in the colour of the filled contours) since tα ≪ Tcyl. For Rm = 200 (lower panel), tα ~ Tcyl, 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 Rm (≳2000) the models start showing changes in parity for t > 40tdiff. Nevertheless, the magnetic energy and the dynamo period, Tcyl, 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.

thumbnail 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.

thumbnail 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 (α0 = 4αc instead of 2αc as in Sect. 3.1) in the interface dynamo model with dynamical α quenching and Rm = 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., . Because of a low value of Rm, the toroidal field generated in the shear layer diffuses into the overshoot layer, and by virtue of the term , 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-6Beq. 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 . 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 ≤ Rm ≤ 2 × 105 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 α2 dynamos can be alleviated by introducing a Fickian diffusive flux in Eq. (11) given by (13)There was an attempt to calculate the diffusion coefficient from direct numerical simulations and it was found that κ ~ 0.3ηt for Rm ~ 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 kf/k1 ~ 10. Their Fig. 7 moreover shows that levels off at high Rm for finite diffusive fluxes even though the value of at which the curve levels off may vary with kf.

thumbnail 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.

thumbnail 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 Bsat either displays underdamped oscillations for Rm = 2 × 103 or goes through very low values for Rm ~ 2 × 105 and takes a long time to relax to a steady amplitude. In this section we introduce a diffusive flux with κ(r) = κ0η(r). Figure 10 shows the results for models with two different values of Rm and κ. For Rm = 2 × 103 and κ0 = 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 Rm = 2 × 105 and κ = 0.01, we observe underdamped oscillations with a final Bsat ~ 0.1Beq. 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 Rm but κ = 0. Figures 11e, f show meridional snapshots of sign(Bφ)(| Bφ | /Beq)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 (Rm = 2 × 105) there exists a critical κc ~ 10-5 such that the secondary dynamo wave does not appear if κ0 < κc. The volume-averaged magnetic energy in this case decays eventually. Note that this threshold for κ is highly dependent on Rm. For instance a model with Rm = 2 × 103 and κ0 = 10-5 produces a dynamo with finite saturation magnetic energy and dynamo wave propagation governed by αM, whereas for κ0 = 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 κ0 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 α2 dynamos have established that a large-scale magnetic field is easily excited on the scale of the system i.e., for a high kf/k1 ratio (Archontis et al. 2003). However the length scale of the magnetic field in Figs.10b, c and 11a, c, e is comparable to , 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 . 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 Rm = 2 × 103. 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.

thumbnail 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 Rm = 2 × 103 in Fig. 12 is now under-damped, whereas the dynamo fails to generate a finite Bsat for Rm = 2 × 105 even though it initially has the same growth rate. On increasing αBL = 10 ms-1 from 6 ms-1 the saturation curve for Rm = 2 × 105 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 Rm 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 . The same figure shows that the saturation value of the magnetic energy is now almost independent of Rm. 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.

thumbnail 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.

thumbnail 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 for both dynamical quenching with gα = 0 and the algebraic quenching with gα = Rm (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 Rm, the toroidal field is weakly confined below the overshoot layer and therefore αM is generated there by the term in the expansion of the forcing term in Eq. (11). Secondly, for high Rm, 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 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 (~104), 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 α2 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

  1. Archontis, V., Dorch, S. B. F., & Nordlund, Å. 2003, A&A, 397, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Blackman, E. G., & Brandenburg, A. 2002, ApJ, 579, 359 [Google Scholar]
  3. Blackman, E. G., & Brandenburg, A. 2003, ApJ, 584, L99 [NASA ADS] [CrossRef] [Google Scholar]
  4. Blackman, E. G., & Field, G. B. 2000, MNRAS, 318, 724 [NASA ADS] [CrossRef] [Google Scholar]
  5. Blackman, E. G., & Field, G. B. 2004, Phys. Plasmas, 11, 3264 [NASA ADS] [CrossRef] [Google Scholar]
  6. Brandenburg, A., & Dobler, W. 2001, A&A, 369, 329 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Brandenburg, A., & Subramanian, K. 2005, AN, 326, 400 [Google Scholar]
  8. Brandenburg, A., Bigazzi, A., & Subramanian, K. 2001, MNRAS, 325, 685 [NASA ADS] [CrossRef] [Google Scholar]
  9. Brandenburg, A., Rädler, K.-H., & Schrinner, M. 2008, A&A, 482, 739 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Brandenburg, A., Candelaresi, S., & Chatterjee, P. 2009, MNRAS, 398, 1414 [NASA ADS] [CrossRef] [Google Scholar]
  11. Cattaneo, F., & Hughes, D. W. 1996, PRE, 54, R4532 [Google Scholar]
  12. Charbonneau, P., 2005, Liv. Rev. Sol. Phys., 2, http://www.livingreviews.org/lrsp-2005-2 [Google Scholar]
  13. Chatterjee, P., Nandy, D., & Choudhuri, A. R. 2004, A&A, 427, 1019 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Chatterjee, P., Brandenburg, A., & Guerrero, G. 2010, GAFD, in press, [arXiv:1005.5708] [Google Scholar]
  15. Choudhuri, A. R., Schüssler, M., & Dikpati, M. 1995, A&A, 303, L29 [NASA ADS] [Google Scholar]
  16. Covas, E., Tworkowski, A., Brandenburg, A., & Tavakol, R. 1997, A&A, 317, 610 [NASA ADS] [Google Scholar]
  17. Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508 [NASA ADS] [CrossRef] [Google Scholar]
  18. Dikpati, M., & Choudhuri, A. R. 1994, A&A, 291, 975 [NASA ADS] [Google Scholar]
  19. Durney, B. R. 1995, Sol. Phys., 160, 213 [Google Scholar]
  20. Gressel, O. 2010, MNRAS, 405, 41 [NASA ADS] [CrossRef] [Google Scholar]
  21. Guerrero, G., & de Gouveia Dal Pino, E. M. 2008, A&A, 485, 267 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Guerrero, G., Chatterjee, P., & Brandenburg, A. 2010, MNRAS, in press, [arXiv:1005.4818] (GCB) [Google Scholar]
  23. Hubbard, A., & Brandenburg, A. 2010, GAFD, submitted, [arXiv:1004.4591] [Google Scholar]
  24. Jepps, S. A. 1975, JFM, 67, 625 [Google Scholar]
  25. Jouve, L., Brun, A. S., Arlt, R., et al. 2008, A&A, 483, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2008, A&A, 491, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Kleeorin, N., & Rogachevskii, I. 1999, PRE, 59, 6724 [Google Scholar]
  28. Kleeorin, N., Moss, D., Rogachevskii, I., & Sokoloff, D. 2000, A&A, 361, L5 [NASA ADS] [Google Scholar]
  29. MacGregor, K. B., & Charbonneau, P. 1997, ApJ, 486, 484 [NASA ADS] [CrossRef] [Google Scholar]
  30. Markiel, J. A., & Thomas, J. H. 1999, ApJ, 523, 827 [NASA ADS] [CrossRef] [Google Scholar]
  31. Mitra, D., Candelaresi, S., Chatterjee, P., Tavakol, R., & Brandenburg, A. 2010, Astron. Nachr., 331, 130 [Google Scholar]
  32. Ossendrijver, M., Stix, M., Brandenburg, A., & Rüdiger, G. 2002, A&A, 394, 735 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Parker, E. N. 1955, ApJ, 122, 293 [NASA ADS] [CrossRef] [Google Scholar]
  34. Parker, E. N. 1993, ApJ, 408, 707 [Google Scholar]
  35. Petrovay, K., & Kerekes, A. 2004, MNRAS, 351, L59 [NASA ADS] [CrossRef] [Google Scholar]
  36. Pouquet, A., Frisch, U., & Léorat, J. 1976, JFM, 77, 321 [NASA ADS] [CrossRef] [Google Scholar]
  37. Shukurov, A., Sokoloff, D., Subramanian, K., & Brandenburg, A. 2006, A&A, 448, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Spiegel, E. A., & Zahn, J.-P. 1992, A&A, 265, 106 [NASA ADS] [Google Scholar]
  39. Subramanian, K., & Brandenburg, A. 2006, ApJ, 648, L71 [Google Scholar]
  40. Sur, S., Shukurov, A., & Subramanian, K. 2007, MNRAS, 377, 874 [NASA ADS] [CrossRef] [Google Scholar]
  41. Sur, S., Brandenburg, A., & Subramanian, K. 2008, MNRAS, 385, L15 [NASA ADS] [CrossRef] [Google Scholar]
  42. Tobias, S. M. 1996, ApJ, 467, 870 [NASA ADS] [CrossRef] [Google Scholar]
  43. van Ballegooijen, A. A., & Choudhuri, A. R. 1988, ApJ, 333, 965 [NASA ADS] [CrossRef] [Google Scholar]
  44. Vainshtein, S. I., & Cattaneo, F. 1992, ApJ, 393, 165 [NASA ADS] [CrossRef] [Google Scholar]
  45. Vishniac, E. T., & Cho, J. 2001, ApJ, 550, 752 [NASA ADS] [CrossRef] [Google Scholar]
  46. Zhang, H., Sokoloff, D., Rogachevskii, I., et al.2006, MNRAS, 365, 276 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

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

In the text
thumbnail 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
thumbnail 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
thumbnail Fig. 4

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

In the text
thumbnail 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
thumbnail Fig. 6

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

In the text
thumbnail Fig. 7

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

In the text
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.