ON THE CAUSE OF SOLAR-LIKE EQUATORWARD MIGRATION IN GLOBAL CONVECTIVE DYNAMO SIMULATIONS

Just over 50 yr after the paper by Maunder (1904), in which he showed for the first time the equatorward migration (EM) of sunspot activity in a time–latitude (or butterfly) diagram, Parker (1955) proposed a possible solution: migration of an αΩ dynamo wave along lines of constant angular velocity Ω (Yoshimura 1975). Here, α is related to kinetic helicity and is positive (negative) in the northern (southern) hemisphere (Steenbeck et al. 1966). To explain EM, ${\boldsymbol {\nabla }} {}\Omega$ must point in the negative radial direction. However, application to the Sun became problematic with the advent of helioseismology showing that ∇_rΩ is actually positive at low latitudes where sunspots occur (Schou et al. 1998), implying poleward migration (PM). This ignores the near-surface shear layer where a negative ∇_rΩ (Thompson et al. 1996) could cause EM (Brandenburg 2005). An alternative solution was offered by Choudhuri et al. (1995), who found that in αΩ dynamo models with spatially separated induction layers the direction of migration can also be controlled by the direction of meridional circulation at the bottom of the convection zone, where the observed poleward flow at the surface must lead to an equatorward return flow. Finally, even with just uniform rotation, i.e., in an α² dynamo as opposed to the aforementioned αΩ dynamos, it may be possible to obtain EM due to the fact that α changes sign at the equator (Baryshnikova & Shukurov 1987; Rädler & Bräuer 1987; Mitra et al. 2010; Warnecke et al. 2011).

Meanwhile, global dynamo simulations driven by rotating convection in spherical shells have demonstrated not only the production of large-scale magnetic fields, but, in some cases, also EM (Käpylä et al. 2012, 2013; Warnecke et al. 2013b; Augustson et al. 2013). Although this seemed to be successful in reproducing Maunder's observation of EM, the reason remained unclear. Noting the agreement between their simulation and the α² dynamo of Mitra et al. (2010) in terms of the π/2 phase shift between poloidal and toroidal fields near the surface, as well as their similar amplitudes, Käpylä et al. (2013) suggested such an α² dynamo as a possible underlying mechanism. Yet another possibility is that α can change sign if the second term in the estimate for α (Pouquet et al. 1976),

$$\begin{equation} \alpha ={\tau _{\rm c}\over 3}\left({-} \overline{{\boldsymbol {\omega }} {}\cdot {\boldsymbol {u}} {}} +{\overline{{\boldsymbol {j}} {}\cdot {\boldsymbol {b}} {}}\over \overline{\rho }}\right), \end{equation} \tag{ 1 }$$

becomes dominant near the surface, where the mean density $\overline{\rho }$ becomes small. Here, ${\boldsymbol {\omega }} {}={\boldsymbol {\nabla }} {}\times {\boldsymbol {u}} {}$ is the vorticity, ${\boldsymbol {u}} {}$ is the small-scale velocity, ${\boldsymbol {j}} {}={\boldsymbol {\nabla }} {}\times {\boldsymbol {b}} {}/\mu _0$ is the current density, ${\boldsymbol {b}} {}$ is the small-scale magnetic field, μ₀ is the vacuum permeability, τ_c is the correlation time of the turbulence, and overbars denote suitable (e.g., longitudinal) averaging. However, an earlier examination by Warnecke et al. (2013a) showed that the data do not support this idea, i.e., the contribution from the second term is not large enough. Furthermore, the theoretical justification for Equation (1) is questionable (Brandenburg et al. 2008).

A potentially important difference between the models of Käpylä et al. (2012, 2013) and those of other groups (Ghizaru et al. 2010; Racine et al. 2011; Brown et al. 2011; Augustson et al. 2012, 2013; Nelson et al. 2013) is the use of a blackbody condition for the entropy and a radial magnetic field on the outer radial boundary. The latter may be more realistic for the solar surface (Cameron et al. 2012).

It should be noted that a near-surface negative shear layer similar to the Sun was either not resolved in the simulations of Käpylä et al. (2012, 2013), or, in the case of Warnecke et al. (2013b), such a layer did not coincide with the location of EM. Instead, most of these simulations show a strong tendency for the contours of angular velocity to be constant on cylinders. Some of them even show a local minimum of angular velocity at mid-latitudes. Indeed, Augustson et al. (2013) identified EM with the location of the greatest latitudinal shear at a given point in the cycle and find that weak negative radial shear also plays a role.

In this Letter, we show through detailed comparison among four models that it is this local minimum, where ∇_rΩ < 0 and α > 0, which explains the EM as a Parker dynamo wave traveling equatorward. While we do not expect this to apply to Maunder's observed EM in the Sun, it does clarify the outstanding question regarding the origin of EM in the simulations. A clear understanding of these numerical experiments is a prerequisite for a better understanding of the processes causing EM in the Sun.

To isolate effects arising from changes in the α effect and the differential rotation, we consider four models. Our reference model, Run I, is the same model presented in Käpylä et al. (2012) as their Run B4m and in Käpylä et al. (2013) as their Run C1. Run II is a run in which the subgrid scale (SGS) Prandtl number, ${\rm Pr}_{\rm SGS}=\nu /\overline{\chi }_{\rm SGS}$, is reduced from 2.5 to 0.5, and the magnetic Prandtl number,  Pr_M = ν/η, is reduced from 1 to 0.5. Here ν is the viscosity, η is the magnetic diffusivity, and $\overline{\chi }_{\rm SGS}$ is the mean SGS heat diffusivity. We keep ν fixed so the effect of lowering  Pr_SGS is that the SGS diffusion is more efficiently smoothing out entropy variations. For stars, the relevant value of  Pr_SGS is well below unity, but such cases are difficult to simulate numerically. In Runs III and IV, we have replaced the outer radiative boundary condition by a cooling layer (see Warnecke et al. 2013b, for the implementation and the profile) above fractional radii r/R = 0.985 and 1.0, respectively; see Figures 1(a) and (b) for the radial temperature and density profiles. Here, R is the solar radius. The cooling profile of Run III leads to a stronger density decrease and suppression of u_rms than in the other runs; see Figures 1(b) and (c). Besides the differences in the fluid and magnetic Prandtl numbers (Run II) and in the upper thermal boundary condition (Runs III and IV), the setups are equal. We can therefore isolate the origin of the difference in the migration pattern of the toroidal field.

Zoom In Zoom Out Reset image size

Our simulations are done in a wedge |90° − θ| ⩽ 75°, 0 < ϕ < 90°, and R − ΔR ⩽ r ⩽ R + δR_C, where θ is colatitude, ϕ is longitude, r is radius, ΔR = 0.3R, and δR_C = 0.01R is the extension by the cooling layer in Runs III and IV. We solve the equations of compressible magnetohydrodynamics using the Pencil Code.⁶ The basic setup of these four models is identical to previous work; the details can be found in Käpylä et al. (2013) and Warnecke et al. (2013b). We scale our results to physical units following Käpylä et al. (2013, 2014) and choose a rotation rate of Ω₀ = 5Ω_☉, where Ω_☉ = 2.7 × 10⁻⁶ s is the solar value.

We begin by comparing the evolution of the mean toroidal field $\overline{B}_\phi$ using time–latitude and time–depth diagrams; see Figure 2. In Run I, $\overline{B}_\phi$ migrates equatorward between ±10° and ±40° latitude, and becomes strongly concentrated around r = 0.8–0.9 R. The cycle period is around five years.⁷ In Run IV, the evolution of $\overline{B}_\phi$ is similar to Run I. Therefore, the blackbody boundary condition is not a necessity for EM. In both runs, a poleward migrating high-frequency dynamo wave is superimposed on the EM, as already seen in Käpylä et al. (2012, 2013). By contrast, in Run II, $\overline{B}_\phi$ migrates poleward between ±10° and ±45° latitude. The field is strongest at r = 0.85–0.98 R and the cycle period is about 1.5 yr. This is clearly shorter than the cycle in Runs I and IV, but significantly longer than the superimposed poleward dynamo wave in those runs. In Run III, $\overline{B}_\phi$ has two superimposed field patterns: one with PM similar in frequency and location to that of Runs I and IV, and a quasi-stationary pattern with unchanged field polarity for roughly 20 yr. The poleward migrating field appears in the upper 10% of the convection zone, whereas the non-migrating field is dominant in the lower half of the convection zone.

Zoom In Zoom Out Reset image size

The distribution of $\overline{B}_\phi$ in the meridional plane can be seen in the top row of Figure 3, where we plot the rms of the mean toroidal magnetic field, time-averaged over the saturated stage, $\overline{B}_\phi ^{\rm rms}\equiv \langle \overline{B}_\phi ^2 \rangle ^{1/2}_t$. In Run I, it reaches 4 kG and is concentrated at mid-latitudes and mid-depths. The field structures are aligned with the rotation axis. Additionally, there is a slightly weaker (≈3 kG) field concentration closer to the equator and surface. A similar field pattern can be found in Run IV, but the field concentrations are somewhat weaker. In Run II, $\overline{B}_\phi ^{\rm rms}$ is concentrated closer to the surface with a larger latitudinal extent than in Run I. The shape of the field structure is predominantly aligned with the latitudinal direction. In Run III, there is some near-surface field enhancement similar to Run I, but closer to the equator. However, the maximum of $\overline{B}_\phi ^{\rm rms}$ is near the bottom of the convection zone, although at higher latitudes it occupies nearly the entire convection zone.

Zoom In Zoom Out Reset image size

Next, we compare the differential rotation profiles of the runs; see the bottom row of Figure 3. All runs develop cylindrical contours of constant rotation as a dominant pattern. However, Runs I, III, and IV possess a local minimum of angular velocity, implying the existence of a negative ∇_rΩ, between ±15° and ±40° latitude, which is the same latitude range where EM was found in Runs I and IV. In Run II, the contours of constant angular velocity are nearly cylindrical, but with a slight radial inclination, which is more than in Run I. This is expected due to the enhanced diffusive heat transport and is also seen in other global simulations (e.g., Brun & Toomre 2002; Brown et al. 2008), where  Pr_SGS is closer to or below unity. Unlike in Runs I, III, and IV, there is no local minimum of Ω. This can be attributed to the higher value of the SGS heat diffusivity in Run II, which smoothes out entropy variations, leading to a smoother rotation profile via the baroclinic term in the thermal wind balance (see corresponding plots and discussion in Warnecke et al. 2013b).

Furthermore, we calculate the local dynamo numbers

$$\begin{equation} C_\alpha ={\alpha \, \Delta R\over \eta _{\rm t0}},\quad C_\Omega ={\nabla _{\!r}\Omega \, \Delta R^3\over \eta _{\rm t0}}, \end{equation} \tag{ 2 }$$

where η_t0 = α_MLTH_pu_rms(r, θ)/3 is the estimated turbulent diffusivity with the mixing length parameter α_MLT = 5/3, the pressure scale height H_p, the turbulent rms velocity u_rms(r, θ), and α(r, θ) is estimated using Equation (1); see also Käpylä et al. (2013). In Figure 4, we plot C_α and $C_\Omega$ as functions of radius for 25° latitude for Runs I–IV. The C_α profiles in all the runs are similar: the quantity is almost always positive, except for a narrow and weak dip to negative values at the very bottom of the simulation domain. The only two exceptions are Runs III and IV, where the cooling layer causes C_α to decrease already below (Run III) or just above (Run IV) the surface, becoming even weakly negative there. The reason is a sign change of the kinetic helicity caused by the sign change of entropy gradient. The $C_\Omega$ profiles are similar for Runs I, III, and IV. There are two regions of negative values in the lower and middle part of the convection zone, with positive values near the surface. In the middle of the convection zone, these profiles coincide with clearly positive values of C_α, as required for EM. For Run II, the profiles of $C_\Omega$ are markedly different: despite the negative dip at the bottom of the convection zone, the values of $C_\Omega$ are generally positive and larger in magnitude than for Runs I, III, and IV. This suggests PM throughout most of the convection zone.

Zoom In Zoom Out Reset image size

To investigate this in more detail, we calculate the migration direction ${\boldsymbol {\xi }} {}_{\rm mig}$ as (Yoshimura 1975)

$$\begin{equation} {\boldsymbol {\xi }} {}_{\rm mig}(r,\theta)=-\alpha \hat{\boldsymbol {e}} {}_\phi \times {\boldsymbol {\nabla }} {}\Omega, \end{equation} \tag{ 3 }$$

where $\hat{\boldsymbol {e}} {}_\phi$ is the unit vector in the ϕ-direction. Note that this and our estimated α(r, θ) using Equation (1) is a strong amplification, in general, of the tensorial properties. In all of our runs, α is on average positive (negative) in the northern (southern) hemisphere.

The migration direction ${\boldsymbol {\xi }} {}_{\rm mig}$ is plotted in the top row of Figure 3 for the northern hemispheres of Runs I–IV. The white arrows show the calculated normalized migration direction on top of the color coded $\overline{B}_\phi ^{\rm rms}$ with black contours indicating $\overline{B}_\phi ^{\rm rms}$ = 2.5 kG. In Runs I and IV, Equation (3) predicts EM in the region where the mean toroidal field is the strongest. This is exactly how the toroidal field is observed to behave in the simulation at these latitudes and depths, as seen from Figure 2. The predicted EM in this region is due to α > 0 and ∇_rΩ < 0. Additionally, in a smaller region of strong field closer to the surface and at lower latitudes, the calculated migration direction is poleward. This coincides with the high-frequency poleward migrating field shown in Figure 2. In Run II, due to the absence of a negative ∇_rΩ (see the bottom row of Figure 3), ${\boldsymbol {\xi }} {}_{\rm mig}$ points toward the poles in most of the convection zone, in particular in the region where the field is strongest; see the top row of Figure 3. Here the calculated migration direction agrees with the actual migration in the simulation; see Figure 2. In Run III, there exists a negative ∇_rΩ, but in the region where the toroidal field is strongest, the calculated migration direction is inconclusive. There are parts with equatorward, poleward, and even radial migration. This can be related to the quasi-stationary toroidal field seen in Figure 2. However, in the smaller field concentration closer to the surface and at lower latitudes, the calculated migration direction is also poleward, which seems to explain the rapidly poleward migrating $\overline{B}_{\phi }$ of Run III (Figure 2). This agreement between calculated and actual migration directions of the toroidal field implies that the EM in the runs of Käpylä et al. (2012, 2013) and in Runs I and IV can be ascribed to an αΩ dynamo wave traveling equatorward due to a local minimum of Ω.

To support our case, we compute a two-dimensional histogram of $|\overline{B}_{\phi }|$ and ∇_rΩ in a band from ±15° to ±40° latitude for Runs I and II; see Figures 5(a)–(b). For Run I, the strong (>5 kG) fields correlate markedly with negative ∇_rΩ < 0. For Run II, the strong fields are clearly correlated with positive ∇_rΩ < 0. These correlations have two implications: first, strong fields in these latitudes are related to and most likely generated by radial shear rather than an α effect. Second, the negative shear in Run I is related to and probably the cause of the toroidal field migrating equatorward and the positive shear in Run II is responsible for PM.

Zoom In Zoom Out Reset image size

These indications resulting from the comparison of four different simulation models lead us to conclude that the dominant dynamo mode of all models is of αΩ type, and not, as suggested by Käpylä et al. (2013), an oscillatory α² dynamo. They based their conclusion on the following three indications. (1) The two local dynamo numbers, C_α and $C_\Omega$, had similar values; see Figures 11 and 12 of Käpylä et al. (2013). However, due to an error, a one-third factor was missing in the calculation of C_α, so our values are now three times smaller; see Figure 4. (2) The phase difference of ≈π/2 between $\overline{B}_\phi$ and $\overline{B}_r$ was observed, which agrees with that of an α² dynamo, as demonstrated in Figure 15 of Käpylä et al. (2013). As shown in Figures 5(c) and (d), this is only true close to the surface (r = 0.98 R). At mid-depth (r = 0.84 R), where $\overline{B}_\phi ^{\rm rms}$ is strong, the phase difference is close to 3π/4, as expected for an αΩ dynamo with negative shear; see Figure 15(e) of Käpylä et al. (2013). (3) Poloidal and toroidal fields had similar strengths, as was shown in Figures 15(a) and (b) of Käpylä et al. (2013). Again, however, this is only true near the surface (r = 0.98 R), where B_ϕ has to decrease due to the radial field boundary condition. As shown in Figure 5(e), $\overline{B}_\phi ^{\rm rms}$ and $\overline{B}_r^{\rm rms}$ are comparable only near r = 0.98 R, whereas in the rest of the convection zone, $\overline{B}_\phi ^{\rm rms}$ is around five times larger than $\overline{B}_r^{\rm rms}$. It is still possible that there is a subdominant α² dynamo operating close the surface causing the phase and strength relation found in Käpylä et al. (2013).

Comparing our results with Augustson et al. (2013), their differential rotation profile possesses a similar local minimum of Ω as our Runs I, III, and IV; see their Figure 2(b). This supports the interpretation that an αΩ dynamo wave is the cause of EM also in their case.

Even though the input parameters are similar to those of Runs I and IV, in Run III $\overline{B}_\phi$ does not migrate toward the equator. The only difference between Runs III and IV is the higher surface temperature in the former (Figure 1(a)). As seen from Figure 1(c), this leads to a suppression of turbulent velocities and a sign change of α close to the surface in those latitudes, where EM occurs in Runs I and IV; see also Figure 4. One of the reasons might be the fact that the sign changes. This suppresses the dynamo cycle and causes a quasi-stationary field. Another reason could be $C_\Omega$ in Run III being stronger at the bottom of the convection zone than in the middle (in contrast to Runs I and IV; see Figure 4), which implies a preferred toroidal field generation near the bottom, where the migration direction is not equatorward; see the top row of Figure 3.

By comparing four models of convectively driven dynamos, we have shown that the EM found in the work of Käpylä et al. (2012) and in Run IV of this Letter as well as the PM in Runs II and III can be explained by the Parker–Yoshimura rule. Using the estimated α and determined Ω profiles to compute the migration direction predicted by this rule, we obtain qualitative agreement with the actual simulation in the regions where the toroidal magnetic field is strongest. This result and the phase difference between the toroidal and poloidal fields imply that the mean field evolution in these global convective dynamo simulations can well be described by an αΩ dynamo with a propagating dynamo wave. We found that the radiative blackbody boundary condition is not necessary for obtaining an equatorward propagating field. Even though the parameter regime of our simulations might be far away from the real Sun, analyzing these simulations, and comparing them with, e.g., mean-field dynamo models, will lead to a better understanding of solar and stellar dynamos and their cycles.

We thank Matthias Rheinhardt for useful comments on the manuscript. The simulations have been carried out on supercomputers at GWDG, on the Max Planck supercomputer at RZG in Garching, and in the facilities hosted by the CSC—IT Center for Science in Espoo, Finland, which are financed by the Finnish ministry of education. This work was partially funded by the Max-Planck/Princeton Center for Plasma Physics (J.W) and supported in part by the Swedish Research Council grant Nos. 621-2011-5076 and 2012-5797 (A.B.), and the Academy of Finland Centre of Excellence ReSoLVE 272157 (M.J.K., P.J.K. and J.W.), and grants 136189 and 140970 (P.J.K).