CONVECTIVE BABCOCK–LEIGHTON DYNAMO MODELS

The emergence of magnetic flux through the solar photosphere regulates solar variability and powers space weather. It is clear that this flux originates in the solar interior and is produced by the solar dynamo. However, it is currently unclear what role flux emergence plays in establishing the 22 yr solar activity cycle. Is it an essential ingredient or merely a superficial by-product of deeper-seated dynamics?

One of the principle means by which flux emergence may act as an essential ingredient in the operation of the solar dynamo is through the so-called Babcock–Leighton (BL) mechanism (Babcock 1961; Leighton 1964). The BL mechanism arises through the dynamics of flux emergence. As a buoyant flux tube rises through the convection zone, the Coriolis force induces a twist in the axis of the tube that is manifested upon emergence as a poleward displacement between the trailing and leading edges of a bipolar active region. When the active region subsequently fragments and disperses due to surface convection and meridional flow, the redistribution of vertical flux induces a net electromotive force (emf) that converts mean toroidal field to poloidal field. Although doubts remain about the viability of the BL mechanism as the principle source of poloidal flux, its empirical foundation and robustness have made it an integral element in many recent mean-field dynamo models of the solar cycle (reviewed by Dikpati & Gilman 2009; Charbonneau 2010).

Current simulations of global solar and stellar convection do exhibit sustained dynamo action and, depending on the parameter regime, can produce magnetic cycles (Ghizaru et al. 2010; Racine et al. 2011; Brown et al. 2011). Yet, these simulations do not have sufficient resolution to reliably capture flux emergence or the BL mechanism. Nelson et al. (2011) recently reported the first convective dynamo simulation to exhibit the spontaneous, self-consistent generation of buoyant magnetic flux structures generated by convectively driven rotational shear. However, even these cannot realistically simulate the subtle, multi-scale dynamics of flux emergence and dispersal that underlies the BL mechanism. This would require (1) very low diffusion to form concentrated flux structures, (2) very high resolution to capture the destabilization and coherent rise of those structures, and (3) a realistic depiction of surface convection, meridional flow, and radiative transfer in order to capture the emergence, coalescence, fragmentation, and dispersal of bipolar active regions (Cheung et al. 2010). Each of these is a formidable computational challenge in its own right, stretching the limits of modern supercomputers. No global model is capable of unifying convective dynamos and the BL mechanism through direct numerical simulation.

Thus, the presence of magnetic cycles in global convection simulations demonstrates that the BL mechanism is not a necessary ingredient for cyclic activity. How this may or may not apply to the solar dynamo is a complex, unresolved issue and we make no attempt at a comprehensive discussion here. Our purpose is rather to investigate how flux emergence may alter the behavior of a convective dynamo by means of the BL mechanism. The BL mechanism is modeled using a mean-field parameterization intended to mimic dynamics that cannot be explicitly captured by the simulation itself for the reasons outlined above. Our approach is described in Section 2 and simulation results are presented in Section 3, along with interpretive discussion.

Before proceeding, it is worth emphasizing up front that the simulations we consider here are of a solar-like star rotating three times more rapidly than the Sun (3Ω_☉). This is done because it is a tidy numerical experiment; without BL forcing, this dynamo builds strong large-scale fields that do not undergo cycles. Other parameter regimes, including those at the solar rotation rate, will be considered in future work.

The starting point for our investigation is the convective dynamo simulation that we refer to as Case D3 and that is described in detail by Brown et al. (2010). We refer the reader to that paper for further information on the set up and results of the simulation as well as the anelastic spherical harmonic (ASH) code that is used to solve the equations of magnetohydrodynamics (MHD) in a rotating spherical shell under the anelastic approximation. Solar values are used for the luminosity and the background stratification but the rotation rate is a factor of three faster than the Sun, as mentioned above. The spatial resolution of all simulations reported in this Letter is 96 × 256 × 512 (r, θ, ϕ). The simulation domain spans from r₁ = 0.718 R to r₂ = 0.966 R, where R is the solar radius.

The salient feature of Case D3 that we are interested in here is the presence of persistent toroidal field structures that we term magnetic wreaths. The convection exhibits an intricate, evolving small-scale structure (Figure 1(a)) but produces a substantial differential rotation (Figure 1(b)) that in turn promotes the generation of prominent magnetic wreaths (Figure 1(c)). The mean toroidal field is approximately symmetric about the equator, with one wreath in each hemisphere of opposite polarity. Although continual buffeting by convective motions induces non-axisymmetric and temporal fluctuations, the mean field is remarkably persistent, retaining its essential structure indefinitely (Figure 1(d)). After they are established, the wreaths persist for at least 60 yr (the duration of the simulation) with no sign of abating (Brown et al. 2010). This time interval is much longer than the rotation period (9.3 days), the convective turnover timescale (about 20 days), and the ohmic diffusion time (about 3.6 yr).

Zoom In Zoom Out Reset image size

The simulations described in this Letter are all restarted from the same iteration of Case D3, defined as time t = 0. This includes the unmodified D3 run shown in Figure 1, which was continued beyond the restart iteration for comparison purposes. The absence of initial transients in Figure 1 demonstrates that this simulation has reached a statistically steady state by the mutual reference time t = 0.

The anelastic MHD equations for the conservation of mass, momentum, and thermal energy are solved with no modification. The only difference between the simulations presented here and that presented in Brown et al. (2010) is in the magnetic induction equation where we add an additional poloidal source term S(r, θ) as follows:

$$\begin{equation} \frac{\partial {\bf B}}{\partial t} = \mbox{\boldmath $\nabla \times $}({\bf v}\mbox{\boldmath $\times $}{\bf B}- \eta \mbox{\boldmath $\nabla \times $}{\bf B}+ S \mbox{\boldmath $\hat{\phi }$}). \end{equation} \tag{ 1 }$$

The additional term is intended to mimic the generation of mean poloidal field by the BL mechanism. Following the mean-field BL dynamo model of Rempel (2006), we choose

$$\begin{equation} S(r,\theta) = \alpha f(r) g(\theta) \hat{B}_\phi, \end{equation} \tag{ 2 }$$

where f(r) and g(θ) are radial and latitudinal profiles,

$$\begin{equation} f(r) = \mbox{max}\left[ 0, 1 - \frac{(r-r_2)^2}{d^2}\right] \end{equation} \tag{ 3 }$$

$$\begin{equation} g(\theta) = \frac{3 \sqrt{3}}{2} \sin ^2\theta \cos \theta, \end{equation} \tag{ 4 }$$

and $\hat{B}_\phi$ is a measure of the mean toroidal flux at the base of the convection zone,

$$\begin{equation} \hat{B}_\phi = \int _{r_1}^{r_b} h(r) \langle B_\phi \rangle\, dr. \end{equation} \tag{ 5 }$$

Brackets 〈 〉 denote an average over longitude and h is an averaging kernal given by h(r) = h₀ (r − r₁) (r_b − r). The integration is confined to a region near the base of the convection zone, below r_b = 0.79 R and the normalization h₀ is defined such that $\int _{r_1}^{r_b} h(r) dr = 1$.

Note that the radial profile in Equation (3) is nonzero only near the top of the convection zone, for r > r₂ − d. We choose d = 20 Mm so the poloidal source operates above r = 0.937 R. Thus, the BL term is nonlocal in the sense that the poloidal source near the surface is proportional to the mean toroidal flux near the base of the convection zone. This is typical of BL dynamo models (e.g., Rempel 2006).

The amplitude of the BL term, α, includes an algebraic quenching of the form α = α₀(1 + B²_t/B²_q)⁻¹ where B_q = 1 MG is the quenching field strength and $B_t^2 = (1/2) \int _0^\pi \hat{B}^2 \sin \theta d\theta$. In practice, the fields generated are much less than B_q so the quenching plays little role.

The magnetic diffusivity η in all simulations is the same as in Case D3, varying from 1.56 to 7.69 × 10¹² cm² s⁻¹ from the bottom to the top of the convection zone, $\propto \overline{\rho }^{-1/2}$, where $\overline{\rho }$ is the background density.

The objective of this Letter is to vary the amplitude of the BL term α₀ in order to investigate how flux emergence may alter the nature of the dynamo. We emphasize again that the three-dimensional MHD equations are unmodified apart from the S(r, θ) term in Equation (1), so setting α₀ = 0 reproduces Case D3 as shown in Figure 1 and as described at length in Brown et al. (2010).

Figure 2 shows "butterfly" diagrams (latitude-time plots of the mean toroidal field near the base of the convection zone) for a series of simulations with progressively increasing values of the BL forcing amplitude α₀. For α₀ = 1 m s⁻¹ the BL term is too weak to significantly influence the operation of the dynamo (Figure 2(a)). The axisymmetric component of the poloidal field near the surface does increase but not enough to affect the maintenance of persistent wreaths in the lower convection zone.

Zoom In Zoom Out Reset image size

For α₀ = 10 m s⁻¹ the results are very different. The wreaths essentially vanish within a few years of simulation time (Figure 2(b)). This can be attributed to the differential rotation operating on the mean poloidal field generated by the BL mechanism via the Ω-effect. The sense of the poloidal field produced by the BL term is such that the Ω-effect generates toroidal flux in the upper convection zone that is of opposite polarity to the sense of the wreaths. Convection rapidly mixes this flux, bringing together opposite polarities that are annihilated through ohmic dissipation. The magnetic energy in the mean fields drops rapidly, decaying by a factor of 10⁶ by t  ∼ 30 yr.

Yet, beyond about 30 yr, the magnetic energy begins to rise again and the dynamo is reborn, saturating by about t  ∼ 75 yr (Figure 2(c)) with a very different structure. Persistent toroidal wreaths are again present but they are symmetric about the equator, with two wreaths per hemisphere and a somewhat lower magnetic energy (about 27% less than in Case D3).

Increasing α₀ by another order of magnitude produces prominent magnetic cycles (Figure 2(d)). Toroidal wreaths form at mid-latitudes in each hemisphere and propagate toward the equator, reminiscent of the solar butterfly diagram (Figure 2(e)). However, the cycle period is much shorter than in the Sun—about 6 months compared to 22 yr. Nonlinear modulation of the cycle is evident, with waxing and waning amplitudes and transient, weaker substructure in the butterfly diagram near the equator.

It is clear from Figure 2(e) that reversals in the two hemispheres are not precisely synchronized. A more careful analysis reveals that the phase difference shifts over time, suggesting that the two hemispheres are to some extent decoupled. For example, the northern hemisphere leads the southern hemisphere from t  ∼ 4.5–11.5 yr but the reverse is true for t  ∼ 12–14.5 yr. Similar shifts in the hemispheric phase difference (albeit less pronounced) have been reported in sunspot records; in particular, the southern hemisphere was apparently leading in cycles 18–20 while the northern was leading in cycles 21–23 (S. McIntosh et al. 2012, in preparation).

The symmetry of the dynamo about the equator can be quantified by the parity ${\cal P} = (B_s^2 - B_a^2)/(B_s^2 + B_a^2)$, where B_s and B_a are the symmetric and antisymmetric components of 〈B_ϕ〉, respectively (sampled at r = 0.84 R). The parity of Case D3 and for α₀ = 1 m s⁻¹ ranges between −0.5 and −1 while that for α₀ = 10 m s⁻¹ ranges between positive 0.5 and 1. By contrast, the cyclic dynamo (α₀ = 100 m s⁻¹) shifts between positive and negative parity, as shown in Figure 3. It does not exhibit the high synchronization of the solar cycle which exhibits a persistent negative parity.

Zoom In Zoom Out Reset image size

As expected, the addition of a BL term has a substantial influence on the amplitude and structure of the mean poloidal field. Without it, the poloidal field has a roughly octupolar structure, such that 〈B_r〉 is radially outward near the north pole, inward near the core of the northern wreath, and antisymmetric about the equator (Brown et al. 2010). As noted above, the BL term generates opposing poloidal flux at mid-latitudes near the surface, producing multi-polar structure and enhancing dissipation. This plays an essential role in establishing the cycles of Figures 2(d) and (e) and is evident in movies of the mean-field evolution. The typical amplitude of the mean poloidal field near the surface for α₀ = 100 m s⁻¹, 1–2 kG, is much larger than for α₀ = 0 (∼200 G in Case D3) but the overall (three-dimensional) magnetic energy is about a factor of two smaller.

The transition from steady to cyclic dynamos occurs when the BL term S competes with the fluctuating emf $\langle {\bf v}^\prime \mbox{\boldmath $\times $}{\bf B}^\prime \rangle$, where v' = v − 〈v〉 and B' = B − 〈B〉. This in turn occurs when α₀ becomes comparable to the velocity of the convective motions, V_c. The relevant scale for V_c is the rms value of the meridional components of v', which is about 120 m s⁻¹ near the top of the convection zone in Case D3, decreasing with depth.

A legitimate question is whether the cyclic dynamo in Figures 2(d) and (e) is operating in an essentially mean-field, axisymmetric mode. In other words, if the BL term were to dominate the generation of mean poloidal field and the differential rotation were to dominate the generation of mean toroidal field via the Ω-effect, then one might expect this three-dimensional simulation to behave very similarly to an analogous, axisymmetric mean-field model. Under this mean-field scenario, coupling between the poloidal and toroidal source regions might occur through the nonlocality of the BL term, the mean meridional flow, and the turbulent diffusion η. The primary role of the resolved convective motions would then be to maintain the mean flows.

In order to assess whether this mean-field scenario is indeed a valid interpretation of the cyclic activity shown in Figures 2(d) and (e), we have initiated another simulation in which we have artificially suppressed the fluctuating emf. More specifically, we have replaced the ${\bf v}\mbox{\boldmath $\times $}{\bf B}$ term in Equation (1) with only mean-field induction $\langle {\bf v}\rangle \mbox{\boldmath $\times $}\langle {\bf B}\rangle$. This simulation was restarted from that shown in Figures 2(d) and (e) at t = 8.25 yr with the same parameters. The only difference is the absence of the fluctuating emf.

The dynamo mode changes dramatically, as demonstrated in Figure 4. The non-axisymmetric magnetic field is quickly dissipated by ohmic diffusion, with the corresponding energy decreasing by six orders of magnitude within two years (shear promotes more rapid dissipation than the nominal 3.6 yr diffusion timescale). The total magnetic energy is about a factor of three larger than in the complementary simulation of Figures 2(d) and (e) and is dominated by a strong, axisymmetric toroidal field that is symmetric about the equator (positive parity) and reverses cyclically with a period of about 1.3 yr. The butterfly diagram in Figure 4 suggests poleward propagation but closer inspection of the mean-field evolution reveals a dynamo wave that propagates toward the rotation axis, with a cylindrical orientation for the wave front. This is consistent with the cylindrical nature of the differential rotation profile (Figure 1(b)) but is strikingly different from the complementary simulation in Figures 2(d) and (e) which exhibits virtually no mean toroidal field at the equator even when the parity is positive. In short, cycles are longer, stronger, and less solar-like without a turbulent emf.

Zoom In Zoom Out Reset image size

Thus, the cyclic dynamo in Figures 2(d) and (e) is not operating as a mean-field dynamo, or at least not in a naive sense. Convection contributes mean-field generation and transport that plays an essential role in shaping the dynamo even for α₀ as high as 100 m s⁻¹. Parity selection in cyclic dynamos is a subtle issue and space limitations preclude a thorough discussion here. We only remark that mean-field BL models have demonstrated that the relatively homogeneous nature of turbulent pumping (Guerrero & de Gouveia Dal Pino 2008), the antisymmetric structure of a supplementary α-effect operating near the base of the convection zone (Dikpati & Gilman 2001) and efficient turbulent mixing that peaks in the upper convection zone (Hotta & Yokoyama 2010) can all help promote negative (dipolar) parity. These mean-field results are loosely consistent with the three-dimensional convection simulations reported here but warrant a more careful investigation which we reserve for future work.

In summary, we have shown that flux emergence can promote cyclic magnetic activity in a convective dynamo by means of the BL mechanism. Although this result is to some extent anticipated by mean-field dynamo models, we have demonstrated it for the first time in a global convection simulation. The BL mechanism is not required to achieve cycles in convective dynamo simulations but, as we have shown here, it may help shape cycle properties such as the period, amplitude, and equatorward propagation of toroidal flux.

However, achieving cyclic activity through the BL mechanism is not easy. To make a significant impact on the dynamo, the amplitude of the BL α-effect must be comparable to the convective velocity, V_c. In particular, we find that α₀  > 10 m s⁻¹ is required to induce cyclic activity when V_c  ∼ 100 m s⁻¹ in the upper convection zone. By comparison, typical values of α used in mean-field BL dynamo models of the solar cycle are less than 1 m s⁻¹ (e.g., Dikpati & Gilman 2009; Charbonneau 2010).

The relatively large value of α₀ used here, together with the strong shear $\vert {\boldmath \nabla } \Omega \vert$ and the efficiency of convective transport, can likely account for the very short cycle period in this young, rapidly rotating Sun (Ω = 3Ω_☉). Might other parameters produce a 22 yr period comparable to the solar cycle? This remains to be seen. Whether the cycle period is regulated by a dynamo wave or flux transport, the large value of α₀ needed to promote cyclic activity (>10 m s⁻¹) and the short convective turnover timescale (∼20 days) may favor relatively short cycle periods. Longer cycle periods can generally be achieved in mean-field BL dynamo models by reducing the efficiency of poloidal flux transport (e.g., Dikpati & Gilman 2009; Charbonneau 2010) but we do not have that freedom here. Here the efficiency of transport is set by the convection, which in turn is an output of the simulation, regulated mainly by factors such as the stellar mass and luminosity that are set by observations. This implies either that the mean-field generation and transport by solar convection is much less efficient than in the models considered here (due possibly to dynamical quenching of the turbulent α and β-effects or overestimation of the convective velocity), or that the solar dynamo does not adhere to a simple BL paradigm.

We thank Mausumi Dikpati and Matthias Rempel for insightful discussions and comments on the manuscript and we thank the anonymous referee for a prompt and constructive report. This work is supported by NASA grants NNH09AK14I (SR&T) and NNX08AI57G (HTP) and computing resources were provided by the NASA High-End Computing (HEC) Program and the NSF-sponsored Teragrid resources at NICS and TACC. B. P. Brown is supported in part by NSF Astronomy and Astrophysics postdoctoral fellowship AST 09-02004. NCAR is sponsored by NSF and CMSO is supported by NSF grant PHY 08-21899.