CYCLIC MAGNETIC ACTIVITY DUE TO TURBULENT CONVECTION IN SPHERICAL WEDGE GEOMETRY

The solar magnetic field exhibits a quasi-periodic cycle with a period of approximately 22 years. This cycle is manifested by the appearance of sunspots in low-latitude activity belts that migrate toward the equator as the sunspot cycle progresses. Reproducing this behavior remains a major challenge to theoreticians. Mean-field models, where small-scale turbulent effects are parameterized (e.g., Krause & Rädler 1980), have reproduced many aspects of the solar cycle, but with broadly varying assumptions for the various parameterizations (see, e.g., Dikpati & Charbonneau 1999; Ossendrijver 2003; Käpylä et al. 2006; Kitchatinov & Olemskoy 2012).

Another, computationally much more demanding, but physically more consistent route is to solve the equations of magnetohydrodynamics directly without resorting to ill-defined parameterizations for the small scales. In practice, however, realistic Reynolds and Rayleigh numbers, describing the effects of molecular diffusion with respect to advection, are not accessible to simulations (e.g., Chan & Sofia 1986; Miesch & Toomre 2009; Käpylä 2011). The usual approach is to enhance the diffusion coefficients to levels that are computationally feasible while striving to maximize the resolution.

Early spherical shell simulations were able to produce a solar-like rotation profile, i.e., one with "equatorward acceleration," and oscillatory large-scale dynamos (Gilman 1983; Glatzmaier 1985). However, the direction of propagation of the dynamo wave was toward the poles, in contradiction to the Sun. This can be qualitatively explained by a Parker dynamo wave with positive radial shear near the equator in conjunction with negative kinetic helicity density, or a positive α-effect, in the northern hemisphere (Parker 1955). More sophisticated simulations with solar rotation rate and luminosity have failed to produce strong large-scale magnetic fields (Brun et al. 2004) or clear cyclic behavior (Miesch et al. 2011). These runs omitted a stable layer below the convection zone. When such a layer is added, non-oscillatory large-scale fields are also found for solar parameters (Browning et al. 2006). Later, oscillatory solutions were obtained from similar simulations with subgrid-scale modeling (Ghizaru et al. 2010; Racine et al. 2011). These are the most solar-like solutions so far, but also in them the activity is at too high latitudes and the activity belts do not propagate toward the equator. When the rotation rate is increased from the solar value in runs without an overshoot layer, first stable wreaths of strong large-scale fields appear (Brown et al. 2010), and at even more rapid rotation, poleward migrating activity is found (Käpylä et al. 2010; Brown et al. 2011).

We report here results from simulations of turbulent convection in spherical wedge geometry with solar-like equatorward acceleration that exhibit, for the first time, equatorward migrating magnetic activity near the equator and a polar branch at high latitudes. The numerical model is the same as that in Käpylä et al. (2011a) but here we add magnetic fields.

We model a segment of a star, i.e., a "wedge," in spherical polar coordinates, where (r, θ, ϕ) denote radius, colatitude, and longitude. The radial, latitudinal, and longitudinal extents of the wedge are 0.7 R ⩽ r ⩽ R, θ₀ ⩽ θ ⩽ π − θ₀, and 0 ⩽ ϕ ⩽ ϕ₀, respectively, where R is the radius of the star. Here we take θ₀ = π/8 and ϕ₀ = π/2.

We solve the compressible hydromagnetics equations,

$$\begin{equation} \frac{\partial \bm A}{\partial t} = {\bm u}\times {\bm B} - \eta \mu _0 {\bm J}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{D \ln \rho }{Dt} = -\bm \nabla \cdot \bm {u}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{D\bm {u}}{Dt} = \bm {g} -2\bm \Omega \times \bm {u}+\frac{1}{\rho } \left(\bm {J}\times \bm {B}-\bm \nabla p +\bm \nabla \cdot 2\nu \rho \bm {\mathsf {S}}\right), \end{equation} \tag{ 3 }$$

$$\begin{equation} T\frac{D s}{Dt} = \frac{1}{\rho }\left[\bm \nabla \cdot (K \bm \nabla T + \chi _{\rm t}\rho T \bm \nabla s) + 2\nu \bm {\mathsf {S}}^2\right], \end{equation} \tag{ 4 }$$

where ${\bm A}$ is the magnetic vector potential, $\bm {u}$ is the velocity, ${\bm B} =\bm \nabla \times {\bm A}$ is the magnetic field, ${\bm J} =\mu _0^{-1}\bm \nabla \times {\bm B}$ is the current density, η is the magnetic diffusivity, μ₀ is the vacuum permeability, $D/Dt = \partial /\partial t + \bm {u} \cdot \bm \nabla$ is the advective time derivative, ν is the kinematic viscosity, K is the radiative heat conductivity, χ_t is the unresolved turbulent heat conductivity, ρ is the density, s is the specific entropy, T is the temperature, and p is the pressure. The fluid obeys the ideal gas law with p = (γ − 1)ρe, where γ = c_P/c_V = 5/3 is the ratio of specific heats at constant pressure and volume, respectively, and e = c_VT is the internal energy. The gravitational acceleration is $\bm {g}=-GM\hat{\bm {r}}/r^2$, where G is the gravitational constant, M is the mass of the star, and $\hat{\bm {r}}$ is the unit vector in the radial direction. We omit the centrifugal force (cf. Käpylä et al. 2011b). The rate of the strain tensor $\bm {\mathsf {S}}$ is given by $\mathsf {S}_{ij}={\textstyle {(1/2)}}(u_{i;j}+u_{j;i}) -{\textstyle {(1/3)}}\delta _{ij}\bm \nabla \cdot \bm {u}$, where the semicolons denote covariant differentiation (Mitra et al. 2009).

2.1. Initial and Boundary Conditions

The initial state is isentropic and the hydrostatic temperature gradient is ∂T/∂r = −g/[c_V(γ − 1)(m + 1)], where m = 1.5 is the polytropic index. We fix the value of ∂T/∂r on the lower boundary. The density profile follows from hydrostatic equilibrium. The heat conduction profile is chosen so that radiative diffusion is responsible for supplying the energy flux in the system, with K decreasing more than two orders of magnitude from bottom to top (Käpylä et al. 2011a). A weak random small-scale seed magnetic field is taken as initial condition (see below).

The radial and latitudinal boundaries are taken to be impenetrable and stress free; see Equations (14) and (15) of Käpylä et al. (2011b). For the magnetic field we assume perfect conductors at the lower radial and latitudinal boundaries, and radial field at the outer radial boundary; see Equations (15)–(17) of Käpylä et al. (2010). On the latitudinal boundaries we assume that the thermodynamic quantities have zero first derivatives, thus suppressing heat fluxes through the boundaries.

On the upper boundary we apply a blackbody condition

$$\begin{equation} \sigma T^4 = -K\frac{\partial T}{\partial r} - \chi _{\rm t} \rho T \frac{\partial s}{\partial r}, \end{equation} \tag{ 5 }$$

where σ is the Stefan–Boltzmann constant. We use a modified value for σ that takes into account that our Reynolds and Rayleigh numbers are much smaller than in reality, so K and therefore the flux are much larger than in the Sun.

2.2. Dimensionless Parameters

We obtain non-dimensional quantities by choosing R = GM = ρ₀ = c_P = μ₀ = 1, where ρ₀ is the initial density at 0.7 R. Our simulations are defined by the energy flux imposed at the bottom boundary, F_b = −(K∂T/∂r)|_{r = 0.7R}, the temperature at the top boundary, T₁ = T(r = R), as well as the values of Ω₀, ν, η, and χ_tm = χ_t(r_m = 0.85 R). The corresponding non-dimensional input parameters are the luminosity parameter $\mathcal {L} = L_0/[\rho _0 (GM)^{3/2} R^{1/2}]$, the normalized pressure scale height at the surface, ξ = [(γ − 1)c_VT₁]GM/R, the Taylor number Ta = (2ΩR²/ν)², the Prandtl number Pr = ν/χ_tm, the magnetic Prandtl number Pm = ν/η, and the non-dimensional viscosity $\tilde{\nu }=\nu /\sqrt{GMR}$. Other useful diagnostic parameters are the Reynolds number Re = u_rms/νk_f and the Coriolis number Co = 2Ω₀/u_rmsk_f, where u_rms = ((3/2)〈u²_r + u²_θ〉)^1/2 is the rms velocity. Note that for u_rms we omit the contribution from the azimuthal velocity, because its value is dominated by effects from the differential rotation (Käpylä et al. 2011b). The Taylor number can also be written as Ta = Co²Re²(k_fR)⁴, with k_fR ≈ 21. Due to the fact that the initial stratification is isentropic, we quote the (semi-) turbulent Rayleigh number Ra_t from the thermally relaxed state of the run,

$$\begin{equation} {\rm Ra}_{\rm t}\!=\!\frac{GM(\Delta r)^4}{\nu \chi _{\rm tm}R^2} \bigg (-\frac{1}{c_{\rm P}}\frac{{d}s}{{d}r} \bigg)_{r_{\rm m}}, \end{equation} \tag{ 6 }$$

where k_f = 2π/Δr is an estimate of the wavenumber of the largest eddies and Δr = 0.3 R is the thickness of the layer. The magnetic field is expressed in equipartition field strengths, $B_{\rm eq}(r)=\langle \mu _0 \rho \bm {u}^2 \rangle ^{1/2}_{\theta \phi }$, where the subscripts indicate averaging over θ and ϕ with azimuthally averaged mean flows subtracted.

The simulations were performed with the Pencil Code,⁴ which is a high-order finite difference method for solving the compressible equations of magnetohydrodynamics.

Our primary simulation (Run B4m) is continued from a thermally relaxed snapshot of a hydrodynamic Run B4 of Käpylä et al. (2011a) with $\mathcal {L}=3.8\times 10^{-5}$, ξ = 0.02, Ta ≈ 1.4 × 10¹⁰, $\tilde{\nu }=2.9\times 10^{-5}$, and Pr = 2.5, resulting in Re = 36, Co = 7.6, and Ra_t ≈ 3 × 10⁶. The discussion of the results refers to this run unless stated otherwise. We also consider two other runs with Co = 4.7 and Re = 39 (Run B3m), as well as Co = 14.8 and Re = 31 (Run B5m). The former is continued from Run B3 of Käpylä et al. (2011a) whereas the latter is run from the initial conditions stated above. Our seed magnetic field has an amplitude of ≈10⁻⁴ B_eq. As a starting point, we use Pm = 1 and a resolution of 128 × 256 × 128 mesh points, but also study the magnetic Prandtl number dependence by continuing Run B4m with values Pm = (0.25, 0.5). In Run B4m the magnetic field grows exponentially over roughly 1500 convective turnover times before reaching the saturated stage in which the total rms magnetic field is B_rms = 0.72 B_eq.

The convection pattern near the surface shows smaller scales at high latitudes and larger elongated structures or "banana cells" near the equator. Figure 1 shows the rotation profile in the saturated regime of the dynamo from Run B4m. The equator rotates faster than the high latitudes and significant radial differential rotation is present only near the equator. In the lower part of the convection zone, $\partial \overline{\Omega }/\partial r$ is negative at low latitudes (just outside the inner tangent cylinder) and positive at high latitudes. The meridional circulation shows several small cells outside the tangent cylinder on both hemispheres. The latitudinal differential rotation, measured by $\Delta _\Omega \equiv (\Omega _{\theta =\theta _0}-\Omega _{\rm eq})/\Omega _{\rm eq}$, where Ω_eq = Ω(θ = π/2), decreases from 0.08 in the kinematic regime to 0.07 in the saturated state. The rotation profiles in Runs B3m and B5m are qualitatively similar.

Zoom In Zoom Out Reset image size

We define mean quantities as averages over longitude and denote them by an overbar. In Run B4m the relative kinetic helicity density $h_{\rm rel} =\overline{{\bm u} \cdot \bm \omega }/u_{\rm rms}\omega _{\rm rms}$, with $\bm \omega = \bm \nabla \times {\bm u}$, is negative (positive) in the northern (southern) hemisphere; see Figure 1. No pronounced sign reversal with depth is seen. The maximum value of h_rel is around 0.3, which allows us to determine the dynamo number describing the strength of the α-effect as C_α = α/η_t0k₁ ≈ h_relk^(ω)_f/k₁ ≈ 2.7, where k^(ω)_f = ω_rms/u_rms is the approximate wavenumber of the energy-carrying eddies, k₁ = π/Δr is the lowest radial wavenumber in the domain, while α ≈ h_relu_rms/3 and η_t0 = u_rms/3k^(ω)_f are estimates for α-effect and turbulent magnetic diffusivity. (We note that χ_t/η_t0 varies between 1.9 near the surface and 0.15 within the convection zone.) The relevant dynamo number characterizing the radial differential rotation is $C_\Omega =\Delta \Omega /\eta _{\rm t0}k_1^2\approx 55$, where we have used ΔΩ/Ω₀ = 0.06 for the normalized radial shear (not to be confused with $\Delta _\Omega$ defined above). The ratio $C_\Omega /C_\alpha$ is well over 10. Following Roberts & Stix (1972), this suggests that we are in what is known as the αΩ regime where shear is strong enough to favor cyclic behavior.

Time series of the averaged longitudinal component of the magnetic field are shown in Figure 2 for different values of Co = (4.7, 7.6, and 14.8). For Co = 7.6, two activity belts are visible; one propagating poleward at high latitudes and another propagating equatorward between 10° and 30° latitude. The equatorward branch at mid-latitudes is also visible for Co = 4.7, but there the cycle is irregular and the magnetic field at low latitudes is weaker. The dynamo mode appears to change to a non-oscillatory one after around 3800 tu_rmsk_f. In the most rapidly rotating case with Co = 14.8 the cyclicity and equatorward migration of the field are clearly present but fluctuations from one cycle to the next are again larger than in Run B4m. In all of the runs a poleward branch with a shorter period is visible near the surface at low latitudes which looks similar to the solution obtained in the nonlinear stage in Käpylä et al. (2010). This mode can be occasionally distinguished in the saturated stage, in particular in Run B5m, but it remains subdominant to the mode with longer period exhibiting equatorward migration.

Zoom In Zoom Out Reset image size

The magnetic field is strongest at r/R ≈ 0.85 and seems to propagate from there to top and bottom of the convection zone; see Figure 3(a), which shows $\overline{B}_\phi$ as a function of r and t being regenerated in the bulk of the convection zone during each cycle. As the field approaches the surface, it propagates equatorward at low latitudes; see Figure 3(b). This mode becomes apparent in the nonlinear phase whereas in the kinematic stage the solution in the bulk of the convection zone does not oscillate; see Figure 3(c) for t u_rmsk_f < 1000. Since $\overline{B}_\phi$ is here normalized by the instantaneous average value, one sees the spatio-temporal structure, and that no reversals occur. Opposite transitions (from oscillatory to quasi-steady) have been observed in Cartesian simulations (Hubbard et al. 2011; Käpylä et al. 2012).

Zoom In Zoom Out Reset image size

On theoretical grounds, we would expect $|\overline{B}_\phi /\overline{B}_r|$ to be of the order of $|C_\Omega /C_\alpha |^{1/2}\approx 4.5$, but the actual ratio is only around unity; see Figure 4. We cannot therefore be certain that the dynamo is really in the αΩ regime, as discussed above. Interestingly, $\overline{B}_r$ shows a greater amplitude at high latitudes while $\overline{B}_\phi$ is stronger at lower latitudes. Furthermore, $\overline{B}_r$ at high latitudes changes sign approximately when $\overline{B}_\phi$ in the low-latitude activity belt changes sign.

Zoom In Zoom Out Reset image size

Visualizations of the toroidal magnetic field from Run B4m near the surface (Figure 5) show a persistent activity belt near the equator which is changing polarity with a period of roughly 400τ, where τ = (u_rmsk_f)⁻¹ is the convective turnover time. We find that decreasing Rm to 18 (Pm = 0.5) shortens the cycle period by roughly 20% to 330τ. At Rm = 9 (Pm = 0.25) the magnetic field decays, but the decay mode is oscillatory with a period of 270τ. Similar increase of the cycle period with Rm has been observed in forced turbulence simulations (Käpylä & Brandenburg 2009) and suggests that the current runs are not in a regime where the molecular diffusivities are unimportant. The same cyclic behavior is seen throughout the depth of the convection zone above r = 0.75 R. Relating the turnover time of our highest Rm model to that of the deep layers of the solar convection zone, i.e., one month, leads to a magnetic cycle period of roughly 33 years. The cycle might well be shorter if the relevant depth is shallower. On the other hand, if we used k^(ω)_f instead of k_f, our cycle period would be 4–5 times longer. It is also noteworthy that $\overline{B}_\phi$ has mixed parity about the equator, except around the time tu_rmsk_f = 2500 when the field is of odd parity; see Figure 4.

Zoom In Zoom Out Reset image size

We have reported solar-like magnetic cycles from simulations of turbulent convection in spherical wedge geometry. The magnetic activity is concentrated in two belts, a high-latitude one propagating poleward, and a low-latitude one propagating equatorward. The strongest magnetic fields, however, occur in the high-latitude activity branch. Simulations with moderately slower and faster rotation show similar behavior. These results will be discussed in more detail in forthcoming publications. Relating the convective turnover time in the simulation to that of the Sun we obtain a cycle period of 33 years which is somewhat longer than that in the Sun and half that obtained by Ghizaru et al. (2010) from quite a different model exhibiting similar solutions, but without equatorward migration. One of the main differences from our earlier work (Käpylä et al. 2010) is that we have omitted a stably stratified overshoot layer beneath. This allowed us to cover almost an order of magnitude larger density contrast within the convectively unstable layer. Furthermore, convective energy transport now dominates over radiative diffusion and a blackbody boundary condition is used for the temperature (cf. Käpylä et al. 2011a). However, compared with the Sun, our contours of differential rotation are still too cylindrical and also the banana-cell pattern of radial velocity might not be realistic. Both may be a consequence of having a large Taylor number; even the turbulent Taylor number, (ν/ν_t)²Ta = 9 Ta/Re² ≈ 10⁸, is rather large. Here, ν_t ≈ η_t0 has been used as an estimate of the turbulent viscosity. The magnetic activity in our model is distributed throughout the convection zone, in contrast to the widely accepted flux-transport dynamo mechanism (Dikpati & Charbonneau 1999) in which a one-cell counterclockwise (north) meridional circulation is crucial. In our model, meridional flows are convergent toward low latitudes (see Figure 1) and may contribute to the resulting equatorward migration. On the other hand, the negative radial differential rotation in the near-surface shear layer (as anticipated by Brandenburg 2005), which is here absent, is not the explanation for the resulting equatorward migration. Clarifying this is an important goal for future work.

The authors acknowledge the anonymous referee for constructive comments on the manuscript. The simulations were performed using the supercomputers hosted by CSC-IT Center for Science Ltd. in Espoo, Finland, who are administered by the Finnish Ministry of Education. Financial support from the Academy of Finland grants Nos. 136189, 140970 (P.J.K.) and 218159, 141017 (M.J.M.), as well as the Swedish Research Council grant 621-2007-4064, and the European Research Council under the AstroDyn Research Project 227952 are acknowledged. The authors thank NORDITA for hospitality during their visits.