ON THE MODE OF DYNAMO ACTION IN A GLOBAL LARGE-EDDY SIMULATION OF SOLAR CONVECTION

A detailed presentation of our simulation framework and numerical integration scheme will be given in a forthcoming publication, so that what follows is meant as an overview. The simulation solves the anelastic form of the ideal MHD equations in a thick, gravitationally stratified spherical shell (0.62 ⩽ r/R_☉ ⩽ 0.96) spanning 3.4 density scale heights and rotating at the solar rate Ω_☉ = 2.69 × 10⁻⁶ rad s⁻¹. The analytic equations governing momentum, energy, and induction evolution are written in symbolic vector form as

$$\begin{eqnarray} &&\frac{D\bm {u}}{Dt} =-\nabla \pi ^\prime - {\bf g}\frac{\Theta ^{\prime }}{\Theta _o} + 2\bm {u} \times {\mathbf \Omega } + \frac{1}{\mu \rho _o}\left(\bm {B}\cdot \nabla \right) \bm {B}, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\frac{D\Theta ^\prime }{Dt} =-\bm {u}\cdot \nabla \Theta _e + {\mathcal H} -\alpha \Theta ^\prime, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\frac{D\!\bm {B}}{Dt}= \left(\bm {B}\cdot \nabla \right) \bm {u} - \bm {B}(\nabla \cdot \bm {u}). \end{eqnarray} \tag{ 3 }$$

In the above expressions $\bm {u}$ and $\bm {B}$ are the flow velocity and magnetic field, respectively, and Θ is the potential temperature, effectively a measure of specific entropy (s = c_pln Θ). The Lagrangian derivative is defined in the usual manner

$$\begin{eqnarray} &&\frac{D}{Dt}\equiv \frac{\partial }{\partial t}+(\bm {u}\cdot \nabla). \end{eqnarray} \tag{ 4 }$$

These evolution equations are also subject to mass continuity and solenoidality constraint on $\bm {B}$:

$$\begin{eqnarray} &&\nabla \cdot \left(\rho _o\bm {u}\right) =0, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\nabla \cdot \bm {B} =0. \end{eqnarray} \tag{ 6 }$$

In the momentum equation (1), π' is a density-normalized pressure perturbation in which magnetic pressure and centrifugal forces have been subsumed. The basic state defining the background stratification is unmagnetized, rigidly rotating, isentropic, and satisfies hydrostatic balance with g∝r⁻². Subscripts "o" refer to this basic state, with primes denoting deviations from a prescribed ambient state, assumed to be in global thermodynamic equilibrium, consistent with helioseismically calibrated solar structural models (e.g., Christensen-Dalsgaard 2002). In our simulation this ambient state (subscript "e") is defined through a combination of polytropes, with a polytropic index varying linearly with r in the stable layer from n = 3.0 at the base of the domain to n = 1.5 at the base of the unstable layer at r/R_☉ = 0.718, and set to n = 1.49995 within the model's convection zone. Radiative diffusion, symbolized by ${\mathcal H}$ on the right-hand side of Equation (2), is expressed in terms of the potential temperature perturbation as

$$\begin{eqnarray} &&{\mathcal H} ={\Theta _o\over \rho _oT_o}\nabla \cdot \left(\kappa _r{\rho _oT_o\over \Theta _o}\nabla \Theta ^\prime\right), \end{eqnarray} \tag{ 7 }$$

with κ_r the coefficient of radiative diffusivity, and T_o(r) the temperature profile associated with the isentropic basic state. Convection is driven through the Newtonian cooling term αΘ' in Equation (2), which forces the system toward the ambient state on a timescale given by the inverse of α (α = 2 × 10⁻⁸ s⁻¹ for the simulation considered here). This procedure, borrowed from simulation of an idealized terrestrial climate (Held & Suarez 1994; Smolarkiewicz et al. 2001), is physically equivalent to imposing a heat flux through the fluid layer. The simulation of Ghizaru et al. (2010), subjected to the residual (perturbation) solar heat flux, is only weakly forced from the convection point of view. Ongoing simulations operating under different forcing regimes indicate that simulated large-scale dynamo behavior can be sensitive to details of the forcing. Such sensitivity is not at all unusual in a turbulent fluid system (Piotrowski et al. 2009).

For numerical implementation, the system of evolutionary Equations (1)–(3) is written as a set of Eulerian conservation laws—using mass continuity and solenoidality, Equations (5) and (6)—and recast in the geospherical coordinate system (Prusa & Smolarkiewicz 2003):

$$\begin{equation} {\partial (\rho ^\ast {\bf \Psi }) \over \partial t} + \nabla \cdot ({\bf V}^* {\bf \Psi }) = \rho ^\ast {\bf R}, \end{equation} \tag{ 8 }$$

where

$$\begin{equation} {\bf \Psi }=\lbrace \bm {u}, \Theta ^\prime, \bm {B}\rbrace \end{equation} \tag{ 9 }$$

denotes the vector of prognosed physical (viz., measurable) dependent variables, and R symbolizes all right-hand-side terms in Equations (1)–(3), including metric forces. The quantity ρ* = Gρ_o combines the anelastic density and the Jacobian of coordinate transformation, and ${\bf V}^\ast =\rho ^\ast {\dot{\bf x}}$ is an effective flux transport term constructed using contravariant velocity ${\dot{\bf x}}$ associated with the actual curvilinear coordinates used. The Lorentz force term in Equation (1) as well as the induction terms in Equation (3) are written in fully conservative form that incorporate $\bm {u}$ and $\bm {B}$ under the divergence operator; cf. Section III in Bhattacharyya et al. (2010). Furthermore, an auxiliary term of the form −∇π* is added to the right-hand side of the induction Equation (3), in order to compensate for the truncation error-level departures from $\nabla \cdot \bm {B}=0$, through the solution, at every time step, of the associated elliptical problem for the auxiliary potential π*. The solution of the resulting conservation laws proceeds as outlined in Ghizaru et al. (2010), following the established EULAG procedures (Prusa & Smolarkiewicz 2003; Prusa et al. 2008).

Viscous, thermal, and Ohmic dissipation are conspicuously absent on the right-hand side of Equations (1)–(3); indeed, with the exception of the radiative diffusion term appearing explicitly on the right-hand side of Equation (2), all dissipation is delegated to the non-oscillatory advection scheme MPDATA (Smolarkiewicz 2006) at the core of EULAG. In essence, the higher-order truncation terms of MPDATA provide an implicit turbulence model (Domaradzki et al. 2003; Margolin et al. 2006; Margolin & Rider 2007). The simulation analyzed in what follows is computed at relatively low spatial resolution (N_r × N_θ × N_ϕ = 47 × 64 × 128), which permits long temporal integrations, yet the low dissipative properties of the underlying computational framework still yield a reasonably turbulent regime, with estimated Reynolds numbers of order 10² and magnetic Prandtl number of order unity.

Simulations typically use a static state as initial condition, with small velocity and magnetic perturbations introduced at t = 0. The integration proceeds at a fixed-size time step, specifically half an hour, for a grand total of some 6 × 10⁶ time steps for the simulation discussed in what follows.

Figure 1 shows a temporal sequence of the toroidal (zonal) magnetic component in the simulation, extracted on a spherical shell corresponding to the core–envelope interface (r/R_☉ = 0.718 in the simulation), plotted in latitude–longitude Mollweide projection. From top to bottom, the sequence runs from one magnetic maximum to the next and is temporally centered on a polarity reversal (middle panel). Despite strong fluctuations in the magnetic field, a large-scale axisymmetric component antisymmetric about the equatorial plane is clearly apparent except at the time of polarity reversal. In fact, except for the obvious reversal of magnetic polarity, the magnetic field distributions at the top and bottom are remarkably similar, showing concentration at mid latitudes and comparable peak strengths (∼0.3 T) in both hemispheres. Polarity reversals occur through a gradual weakening of the large-scale axisymmetric magnetic component, with antisymmetry about the equatorial plane maintained reasonably well as the time of reversal is approached (t = 202.5 years), and establishing itself rapidly again once the next half-cycle starts to build up (t = 221.9 years).

Zoom In Zoom Out Reset image size

The spatiotemporal evolution of the large-scale axisymmetric component is best viewed by zonally averaging the total magnetic field present at any given time in the simulation. Working in spherical coordinates (r, θ, ϕ), where −π/2 ⩽ θ ⩽ π/2 is the latitude (rather than the polar angle), such a zonal average is defined as

$$\begin{equation} \langle \bm {B} \rangle (r,\theta,t)= {1\over 2\pi }\int _0^{2\pi }\bm {B}(r,\theta,\phi,t){d}\phi. \end{equation} \tag{ 10 }$$

The top panel in Figure 2 shows a time–latitude diagram of the toroidal magnetic component (〈B_ϕ〉) averaged in this manner, constructed at a depth corresponding to the core–envelope interface in the model (r/R_☉ = 0.718). If the toroidal magnetic flux ropes assumed to give rise to bipolar active regions are indeed stored at this depth, as suggested by stability analyses (e.g., Ferriz-Mas et al. 1994; Fan 2009), and rise radially to the photosphere, then this is the simulation's equivalent to the sunspot butterfly diagram (e.g., Hathaway 2010). The toroidal magnetic component is concentrated at mid latitudes (30° ⩽ |θ| ⩽ 70°), as opposed to the low latitudes (5° ⩽ |θ| ⩽ 40°) suggested by the butterfly diagram, but does show a tendency for equatorial migration as each half-cycle unfolds. On a time–latitude diagram such as Figure 2 (top panel), this is seen in the strong toroidal field concentrations, which take an elongated, elliptical shape, with the "major axis" tilted toward the equator as the cycle is followed in time. In other words, throughout a cycle the latitude of the peak large-scale toroidal magnetic field occurs at decreasing latitudes, until the cycle terminates and the next one begins anew at higher latitudes.

Zoom In Zoom Out Reset image size

The middle panel of Figure 2 shows the time–radius slice of the same simulation run, extracted in the southern hemisphere at latitude −45°, where the toroidal field is strongest at the core–envelope interface (cf. top panel). Note how polarity reversals begin well within the convection zone, at depth r/R_☉ ≃ 0.8, with radial drift and concentration of the magnetic field both upward as well as downward all the way to the core–envelope interface, where the field reaches its peak strength.

The bottom panel of Figure 2 shows the corresponding time evolution of the zonally averaged radial surface magnetic component, again in a time–latitude diagram. The surface field is characterized by a well-defined dipole moment closely aligned with the rotational axis, with transport of surface fields taking place from lower latitudes, and possibly contributing to the reversal of the dipole moment. Comparing the three panels in Figure 2 reveals that the dipole moment and deep-seated toroidal component reach their peaks at approximately the same time, indicating that they oscillate in phase, without significant temporal lag.

Overall, the magnetic cycle characterizing the large-scale, zonally averaged magnetic component is quite regular, here with a half-period of 30 years, almost thrice the 11 years of the solar cycle. Note the good long-term synchrony maintained between the northern and southern hemispheres, persisting despite significant fluctuations in the amplitude and duration of cycles in each hemisphere. In keeping with solar tradition, we delineate the cycles from one toroidal field polarity reversal to the next. Remember therefore that what we henceforth refer to as "cycle," as numbered on the bottom panel of Figure 2, spans in fact one-half of a complete magnetic cycle. While the average period of our cycles so-defined is almost exactly 30 years for our 11 cycles, this period can become as low as 25 years (simulated cycle 11) or as high as 35 (simulated cycle 4) from one cycle to the next.

Stratification and rotation, the latter through the agency of the Coriolis force, introduce cross-correlation between turbulent velocity components that lead to significant Reynolds stresses powering large-scale flows throughout the solar convection zone. This phenomenon also arises naturally in the type of rotating, stratified convection simulations considered here, in either the hydrodynamical or magnetohydrodynamical regimes (e.g., Miesch et al. 2000; Brun et al. 2004; Browning et al. 2006, and references therein).

The zonal averaging procedure defined in Equation (10) can be used to extract a mean flow $\langle \bm {u}\rangle (r,\theta,t)$ from the simulation output. The result of this procedure is shown in Figure 3, in the form of isocontours of zonally averaged rotational frequency (angular velocity divided by 2π) and meridional flow vectors superimposed on a color coding of the latitudinal mean-flow component 〈u_θ〉. This is a purely hydrodynamical parent simulation subjected to the same thermal forcing as the MHD simulation described above. Differential rotation is reasonably solar-like, with equatorial acceleration and near-radial isocontours at mid to high latitudes, transiting rapidly to rigid rotation across a thin tachocline-like shear layer located immediately beneath the convective layers. The sharpness of this rotational transition layer is noteworthy for these types of simulations; it is a direct reflection of the low dissipation embodied in our numerical scheme and of the rapidly increasing subadiabaticity with depth characterizing our ambient stratification. The pole-to-equator contrast in surface and subsurface rotation rate, ∼100 nHz, is also reasonably solar-like. This differential rotation pattern remains fairly steady over the time span of the simulation. The large-scale meridional flow, in contrast, is much more strongly fluctuating, but when averaged over a sufficient long time span, away from the equatorial regions the net flow is found to be poleward in the surface and subsurface layers, and equatorward near the base of the convection zone. The surface meridional flow, peaking at ∼1 m s⁻¹ at mid to high latitudes, is about a factor of 10 slower than observed on the solar surface. At lower latitudes, the elongated flow structures parallel to the rotation axis are the residual signature of persistent elongated convective rolls stretching across the equator, typical of these types of simulations running in this turbulent regime (see, e.g., Miesch et al. 2000, Figure 5; Brun et al. 2004, Figure 1; Käpylä et al. 2010, Figure 2).

Zoom In Zoom Out Reset image size

Turning now to our MHD simulation, these large-scale flows are now found to carry the signature of the magnetic cycle. This is shown in Figure 4, displaying again the rotational frequency and meridional flow in meridional planes, in the same format and using the same color scales and ranges as in Figure 3, together now with the zonally averaged toroidal magnetic component (rightmost plots). These were all averaged in time over as temporal window of width Δt = 3.3 years, about one-tenth of a magnetic half-cycle, centered over the peak of simulated cycle 1 (t = 16 years, top row) and following cycle minimum (t = 32 years, bottom row). Differential rotation (left) is again characterized by equatorial acceleration, but there now remains little latitudinal gradients at mid latitudes as compared to the parent hydrodynamical simulation. The pole-to-equator contrast in rotation rates is down to 40 nHz at epochs of cycle minimum, and even less at maximum, which is now a far less solar-like situation. Differential rotation spreads somewhat deeper into the stable layer, especially at high latitudes, but a modest tachocline-like shear layer still persists. Significant cycle-driven torsional oscillations materialize throughout the simulation. In Figure 4, this is mostly apparent at high latitudes, with the polar region speeding up at time of cycle maxima, but lower-amplitude torsional oscillations are actually present also in equatorial regions and at all depths.

Zoom In Zoom Out Reset image size

The surface and subsurface meridional flow remains generally poleward at mid to high latitudes, although secondary, counterrotating flow cells develop in the descending phase of the cycles, their onset visible here on the 〈u_θ〉 plot at maximum phase. At the base of the convecting layers, the meridional flow is equatorward at mid to high latitudes at all cycle phases, but does show significant cyclic variations in speed, generally being more vigorous near cycle maximum, reaching half a meter per second at ±60° latitude. The flow pattern near times of "cycle minimum" more closely resembles what is observed in the purely hydrodynamical simulation (cf. Figure 3).

With zonal means of the total flow $\bm {u}$ and magnetic field $\bm {B}$ defined through Equation (10), subtracting these means from the total flow and magnetic field then yields the non-axisymmetric contributions:

$$\begin{equation} \bm {u}^\prime (r,\theta,\phi,t)= \bm {u}(r,\theta,\phi,t)-\langle \bm {u}\rangle (r,\theta,t), \end{equation} \tag{ 11 }$$

$$\begin{equation} \bm {B}^\prime (r,\theta,\phi,t)= \bm {B}(r,\theta,\phi,t)-\langle \bm {B}\rangle (r,\theta,t). \end{equation} \tag{ 12 }$$

In the following section, we shall demonstrate that these non-axisymmetric contributions can be reasonably associated with small spatial scales and are henceforth referred to as such, while the zonal means can be associated with large spatial scales, commensurate with the size of the simulation domain. It is then straightforward to compute the contributions of these various flow and field components to the total specific energy budget, through volume averages over the full simulation domain:

$$\begin{equation} E(\langle u\rangle)=\frac{1}{2}\int _V \langle u\rangle ^2{d}V, \end{equation} \tag{ 13 }$$

$$\begin{equation} E(u^\prime)=\frac{1}{2}\int _V (u^\prime)^2{d}V, \end{equation} \tag{ 14 }$$

$$\begin{equation} E(\langle B\rangle)=\frac{1}{\mu _0}\int _V \rho ^{-1}\langle B\rangle ^2{d}V, \end{equation} \tag{ 15 }$$

$$\begin{equation} E(B^\prime)=\frac{1}{\mu _0}\int _V \rho ^{-1}(B^\prime)^2{d}V, \end{equation} \tag{ 16 }$$

with $\langle u\rangle ^2=\langle \bm {u}\rangle \cdot \langle \bm {u}\rangle$, etc. The comparison between these various contributions is carried out pictorially in Figure 5, in the form of time series spanning the full simulation run. All four contributions are similar in magnitude, and the cycle stands out as a prominent, smooth cyclic variation in the energy contribution of the axisymmetric component of the magnetic field, which varies by a factor of 2–3 between epochs of cycle "maxima" and "minima." The small-scale magnetic energy E(B') does show some in-phase variations with E(〈B〉), but the energy contributions of the non-axisymmetric flow component, dominated by turbulent convection, remain remarkably steady in comparison, showing only short timescale fluctuations about a well-defined mean value. The energy contribution of the axisymmetric flow, dominated by differential rotation, shows variations on timescale commensurate with the cycle of the axisymmetric part of the magnetic field and reflects the signature of torsional oscillations.

Zoom In Zoom Out Reset image size

Mean-field electrodynamics is predicated on the assumption that the fluid flow ($\bm {u}$) and magnetic field ($\bm {B}$) can be separated into "mean" (usually spatially large-scale and slowly varying in time) and "fluctuating" (usually small-scale and rapidly varying) components:

$$\begin{equation} \bm {u}=\langle \bm {u}\rangle +\bm {u}^\prime,\qquad \bm {B}=\langle \bm {B}\rangle +\bm {B}^\prime, \end{equation} \tag{ 17 }$$

where the angular brackets denote an intermediate averaging scale for which

$$\begin{equation} \langle \bm {u}^\prime \rangle =0,\qquad \langle \bm {B}^\prime \rangle =0. \end{equation} \tag{ 18 }$$

Given the well-defined large-scale axisymmetric component present in our simulation, it becomes natural to associate the averaging scale with the zonal average defined through Equation (10), so that the fluctuating components become defined as the non-axisymmetric contributions through Equations (11) and (12). In practice, it is of course possible to define fluctuating flows and fields in this way for any global MHD simulation run. However, for the mean-field electrodynamics approach to be mathematically well-posed and physically meaningful, it is essential for a good separation of scales to hold. Accordingly, we first examine in some detail whether this is indeed the case for our simulation.

Figure 6 shows the mean and fluctuating toroidal magnetic field components resulting from this decomposition applied to the toroidal field distribution plotted on the bottom panel of Figure 1 (simulation time t = 230 years, very near the peak of simulated cycle 8). At the base of the convection zone, both components have similar strength, but the fluctuating component in Figure 6 shows little or no azimuthal structuring on scales comparable to the solar radius, nor any clear hemispheric pattern. In particular, no sign of a significant non-axisymmetric dipolar or quadrupolar components is visible.

Zoom In Zoom Out Reset image size

This visual impression is confirmed by a modal decomposition in spherical harmonics, which reveals significant power concentrated in the axisymmetric (m = 0) mode of low, odd angular degree ℓ. This is illustrated in Figure 7, showing the result of a modal decomposition of the toroidal magnetic component at the core–envelope interface (r/R_☉ = 0.718). The decomposition was carried out at every time step and then averaged for two sets of disjoint data blocks of temporal width 900 days centered either on the epochs of polarity reversals (panel A, for "cycle minima") or peak toroidal field (panel B, for "cycle maxima"). The top row of Figure 7 displays the result of this procedure in the form of color coding of the absolute values of the modal coefficients, in the plane defined by the angular degree ℓ (horizontal) and azimuthal degree m (vertical). The diamond shape results from the truncation of azimuthal modes, chosen so as to ensure comparable peak spatial resolution in the latitudinal and azimuthal directions.

Zoom In Zoom Out Reset image size

Comparing panels A and B of Figure 7 reveals significant power in the axisymmetric (m = 0) angular modes of odd-ℓ angular degree at cycle maximum that all but vanish at times of cycle minimum. This is particularly prominent for the dipolar (ℓ, m) = (1, 0) mode. At high-ℓ values, on the other hand, power is broadly and more evenly distributed in modes of all azimuthal orders m and shows no clear variations with the phase of the cycle for the large-scale field. Indeed, subtracting panel B from panel A, yielding panel C in Figure 7, reveals very little residual power at high wavenumbers.

Figure 7(D) shows the variation of the same (unsigned) modal coefficients as a function of the latitudinal wavenumber ℓ for the axisymmetric modes only (m = 0, solid lines), and the quadratic sum of all coefficients for non-axisymmetric modes ($m\not=0$, dotted lines). The latter are seen to dominate over the former for ℓ ⩾ 10, but are of comparable amplitude for smaller ℓ, except at maximum epochs (red lines) where the modal amplitudes are dominated by the ℓ = 1 and ℓ = 5 axisymmetric modes. Note also how the summed non-axisymmetric coefficients at maximum and minimum epochs (red and green dotted lines) are almost exactly coincident at all ℓ-values. For axisymmetric modes (solid lines), this holds only for ℓ ⩾ 10. In addition, the peaks at ℓ = 1 and l = 5 remain present at almost the same amplitude in the subtracted distribution (black), while the amplitude for other modes drops significantly, by over an order of magnitude at high ℓ-values. In fact, for axisymmetric modes the only significant difference between maximum and minimum phases is suppression of the ℓ = 1 and ℓ = 5 modal amplitudes. All this indicates that the cycle in the large-scale axisymmetric magnetic field does not alter significantly the spectral properties of the small-scale, "turbulent" magnetic component.

Examining time series of individual modal coefficients reveals that the cycle shows up prominently in the (ℓ, m) = (1, 0) and (5, 0) modes, but the time series for $m\not=0$ modes, as well as for even-ℓ axisymmetric modes, have much lower amplitudes and show no sign of persistent cyclic variation. This has been further verified by computing Fourier transforms of these various time series, which show clear peaks at periods of ∼29 and ∼58 years in the odd-ℓ axisymmetric modes, but no peaks and overall power levels inferior by over an order of magnitude for even-ℓ and $m\not= 0$ modes. Similar temporal patterns are observed within the convection zone and/or for the other magnetic components. Overall, the above analyses then suggest that the decomposition embodied in Equations (11) and (12) represents a viable ansatz, at least for this simulation run.

With the assumption of scale separation vindicated at least to some degree, our purpose is now to identify and characterize the dynamo mechanism underlying the observed magnetic cycle described above. The starting point of the analysis is the magnetohydrodynamical induction equation (e.g., Davidson 2001), describing the evolution of a magnetic field $\bm {B}$ subjected to the inductive action of a flow field $\bm {u}$ in addition to Ohmic dissipation of the associated electrical current density:

$$\begin{equation} \frac{\partial \bm {B}}{\partial t} = \bm {\nabla } \times [\bm {u} \times \bm {B} - \eta \bm {\nabla } \times \bm {B}]\,, \end{equation} \tag{ 19 }$$

where η = (μ₀σ_e)⁻¹ is the magnetic diffusivity, inversely proportional to the electrical conductivity σ_e (SI units are used throughout). Note that explicit Ohmic dissipation is now retained, unlike in the induction equation solved by the simulation (cf. Equation (3)). Inserting Equations (17) and performing the zonal average introduced earlier, and remembering that the averaging operator commutes with spatial derivatives pertaining to the large spatial scales, one readily obtains

$$\begin{equation} \frac{\partial \langle \bm {B} \rangle }{\partial t} = \bm {\nabla } \times [\langle \bm {u} \rangle \times \langle \bm {B} \rangle + \bm {\mathcal {E}} - \eta \bm {\nabla } \times \langle \bm {B} \rangle]\,, \end{equation} \tag{ 20 }$$

where

$$\begin{equation} \bm {\mathcal {E}} = \langle \bm {u}^\prime \times \bm {B}^\prime \rangle \end{equation} \tag{ 21 }$$

is the mean EMF due to the fluctuations about the large-scale magnetic field. Note that this mean turbulent EMF is generally nonzero, because the correlation between the fluctuating flow and magnetic field does not necessarily vanish upon averaging, even though $\bm {u}^\prime$ and $\bm {B}^\prime$ individually do by definition. For more details on some of the many subtleties involved, see, e.g., Moffatt (1978), Rüdiger & Hollerbach (2004), Hoyng (2003), and Ossendrijver (2003).

From the simulation results, it is straightforward to reconstruct $\bm {u}^\prime$ and $\bm {B}^\prime$ via Equations (11) and (12), and then compute the EMF by performing the required cross product and zonal averaging at every time step, as per Equation (21). The result of this calculation is presented in Figure 8 for the r (top), θ (middle), and ϕ (bottom) components of the EMF, in the form of time–latitude slices at depth r/R_☉ = 0.85 near the middle of the convection zone (left panels), where polarity reversals are initiated (see Figure 2), and meridional slices extracted at simulation time t = 112 years, at the peak of the fourth magnetic half-cycle (right panels). Note already how the EMF vanishes rapidly as one moves below the core–envelope interface, as expected since convectively driven turbulence disappears there. Otherwise the EMF pervades most of the convection zone. Despite being a very noisy "signal," on the largest spatial scales all EMF components exhibit a well-defined symmetry ($\mathcal {E}_r$, $\mathcal {E}_\phi$) or antisymmetry ($\mathcal {E}_\theta$) about the equatorial plane. Note in particular that a ϕ-component of the EMF positive (negative) in both hemispheres is precisely what is required to support a positive (negative) dipole moment (cf. bottom panels in Figures 2 and 8), as per the right-hand rule, since under the MHD approximation the EMF drives a parallel electrical current density. Moreover, reversing the direction of this EMF every half-cycle then amounts to reversing the sign of that dipole moment, as observed in the simulation.

Zoom In Zoom Out Reset image size

All EMF components show strongly reduced amplitudes in the equatorial regions. This can be traced to a change in the topological character of convective fluid motions, which at low latitudes tend to organize themselves in longitudinal stacks of large, latitudinally elongated convective cells approximately aligned with the rotational axis (see Figure 1(A) in Ghizaru et al. 2010). These are quite typical of these types of global convection simulations (see, e.g., Browning et al. 2006, Figure 1; Brown et al. 2010, Figure 1; Käpylä et al. 2010, Figure 2). Evidence of a quenching of the small-scale flow components can also be found in the spatial distribution of the small-scale surface magnetic field, which shows a significant decrease in both coverage and magnitude at very low latitudes (see, e.g., Figure 1(B) in Ghizaru et al. 2010).

The most striking global spatiotemporal feature visible in Figure 8 is certainly the cyclic variation of all EMF components, with a period identical to that of the cycle observed in the large-scale magnetic field (cf. Figure 2). Here the large-scale velocity field is roughly stationary over the entire simulation, and the small-scale turbulent velocity field is also stationary in a statistical sense (cf. Figure 5). Since the EMF term depends only on the small-scale fluctuations about the large-scale magnetic fields, it is then not at all obvious a priori that the EMF should exhibit the same cyclic evolution as the mean field; the fact that it does already provides us with important clues as to the nature of the large-scale dynamo operating in our simulation. More specifically, it suggests that the turbulent EMF actually plays a key role in the production of the large-scale magnetic component. We now turn to this question, using the mathematical machinery of mean-field electrodynamics.

The essence of mean-field electrodynamics is to provide an expression for the EMF in terms of the large-scale magnetic field, so that the small scales (fluctuations) are effectively removed from the problem. The procedure consists of a simple linear expansion of the EMF as a power series about the large-scale magnetic field and its derivatives, i.e.,

$$\begin{equation} \mathcal {E}_i = \alpha _{ij} \langle B_j \rangle + \beta _{ijk} \partial _j \langle B_k \rangle + \,\, {\rm higher\,\, order \,\, derivatives}, \end{equation} \tag{ 22 }$$

with summation implied over repeated indices. The tensors appearing on the right-hand side of this expression and relating the EMF to the mean field can depend on the properties of small-scale flow and field fluctuations, but not on the mean field itself, a situation expected to hold only if the latter is too weak to impact the dynamics of the fluctuating flow. Even though it is defined globally, the steadiness of the time series of the energy contribution of the small-scale flow on the temporal scale of the cycle in the large-scale field, as seen in Figure 5, is consistent with this working hypothesis. Note also here that since this expansion is linear in the large-scale field, if that field exhibits cyclic properties, then so should the EMF, which is what we observe (cf. Figure 8). This is our central motivation for appealing to mean-field electrodynamics to analyze the dynamo operating deep in the convection zone of our simulation.

The leading-order contribution in the right-hand side of Equation (22) is usually named the "α-effect." The second term in the expansion is responsible, among other effects, for turbulent diffusion. For our first analysis, however, we retain only the α-effect for simplicity. Owing to the rough stationarity of the large-scale flows in the simulation, we also assume that α_ij is time independent, an assumption that will prove itself well-verified a posteriori. We thus write as a first approximation

$$\begin{equation} \mathcal {E}_i(t,r,\theta) = \alpha _{ij}(r,\theta) \langle B_j \rangle (t,r,\theta)\,. \end{equation} \tag{ 23 }$$

One can further decompose α_ij into symmetric and antisymmetric parts, such that the EMF may be rewritten as

$$\begin{equation} \mathcal {E}_i = \alpha _{(ij)} \langle B_j \rangle + (\bm {\gamma } \times \langle \bm {B} \rangle)_i\,, \end{equation} \tag{ 24 }$$

where

$$\begin{equation} \gamma _i = -\frac{1}{2}\epsilon _{ijk} \alpha _{jk} \end{equation} \tag{ 25 }$$

is called the general turbulent pumping velocity, because it effectively provides an additive contribution to the mean, large-scale inductive flow in Equation (20), although its physical origin lies with the fluctuating, small-scale flow and magnetic field. The flow-like vector field $\bm {\gamma }$ so-defined is in general non-solenoidal ($\nabla \cdot \bm {\gamma }\not=0$) and acts on the total magnetic field, with variations between magnetic component subsumed into the off-diagonal terms of the symmetric part of the α-tensor (see Ossendrijver et al. 2002 for a more thorough discussion).

A number of methods have been designed to measure the components of the α-tensor in turbulent MHD simulations (e.g., Ossendrijver et al. 2001, 2002; Käpylä et al. 2006, 2009; Hubbard et al. 2009, and references therein). These method have been designed in the context of MHD simulations that do not produce a well-defined large-scale magnetic component, so that the latter must be imposed externally—and artificially—with the consequence that the α-tensor being measured is not necessarily characterizing the simulation prior to the application of the external field. For more on these—and other—difficulties, see Brandenburg (2009), Cattaneo & Hughes (2009), and references therein. In contrast, we are in the very advantageous position to have in hand a simulation that does generate a large-scale magnetic field, in a manner consistent with the dynamical interaction of flow and field on all numerically resolved spatial scales.

The procedure we employ here to extract α_ij from the simulation data differs from the approaches found in the literature, and we shall therefore first detail it. Essentially, we attack this problem from an experimentalist's point of view, i.e., we have (numerical) data in our possession, namely the EMF and the average magnetic field, and we wish to verify if these data match a specific model, namely the α-effect parameterization of the EMF embodied in Equation (23). This task therefore amounts to a fitting problem, in other words the components of the tensor α_ij are calculated a posteriori by minimizing the difference between the left-hand side and right-hand side of Equation (23). We opt here for a standard least-squares fit, based upon singular value decomposition (hereafter SVD; see Press et al. 1992, Section 15.4). We state results of various theorems without proof.

For a given component, say m (=r, θ or ϕ) of the EMF at a fixed grid point⁴ (r_b, θ_c), we define the time-dependent functions y(t) and X_k(t) as follows:

$$\begin{eqnarray} y(t) &=& \mathcal {E}_m(t,r_b,\theta _c)\,, \end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray} X_k(t) &=& \langle B_k \rangle (t,r_b,\theta _c)\,. \end{eqnarray} \tag{ 26b }$$

We also define

$$\begin{equation} a_k = \alpha _{mk}(r_b,\theta _c)\,, \end{equation} \tag{ 27 }$$

which then allows us to rewrite the parameterization of the EMF as

$$\begin{equation} y(t) = \sum _{k=1}^3 a_k X_k(t) \,. \end{equation} \tag{ 28 }$$

Denoting by N_t the number of time steps of our simulation, we define a merit function as

$$\begin{equation} \chi ^2 = \sum _{i=1}^{N_t} \left[y(t_i) - \sum _{k=1}^3 a_kX_k(t_i)\right]^2\,, \end{equation} \tag{ 29 }$$

where t_i is the value of the time coordinate at the ith step. The goals are to find the three parameters a_k that minimize the merit function and to obtain an estimate of the goodness of fit. The method we employ to carry out these two tasks simultaneously is the SVD of the "design matrix" $\bm {{\rm A}}$, defined as

$$\begin{equation} A_{ij} = X_j(t_i)\,. \end{equation} \tag{ 30 }$$

The design matrix therefore depends on the large-scale magnetic field components. This matrix has N_t rows and three columns, and by virtue of a theorem of linear algebra, can always be decomposed as follows:

$$\begin{equation} \bm {{\rm A}} = \bm {{\rm U}} \cdot \bm {{\rm w}} \cdot \bm {{\rm V}}^T\,, \end{equation} \tag{ 31 }$$

where the matrix $\bm {{\rm U}}$ is an n_t-by-three column orthogonal matrix, $\bm {{\rm w}}$ is a three-by-three diagonal matrix containing the so-called singular values, and $\bm {{\rm V}}$ is a three-by-three orthogonal matrix. The key element is that the solution to the minimization of the merit function (29) is directly obtained in terms of the three matrices $\bm {{\rm U}}$, $\bm {{\rm V}},$ and $\bm {{\rm w}}$. The solution vector $\bm {a} = (a_1,a_2,a_3)$ is given by

$$\begin{equation} \bm {a} = \bm {{\rm V}} \cdot \bm {{\rm w}}^{-1} \cdot \bm {{\rm U}}^T \cdot \bm {y}\,, \end{equation} \tag{ 32 }$$

where the vector $\bm {y} = (y_1, y_2, \ldots, y_{N_t})$ and therefore depends on the selected EMF component. To construct the SVD of our design matrix, we use the algorithm and routines of Press et al. (1992), as implemented in the IDL programming language. By repeating this procedure for all three components of the EMF and for each of the N_r × N_θ grid points (r_b, θ_c), we can calculate the complete tensor α_ij(r, θ).

Figure 9 shows the results of this fitting procedure over the complete 337 year simulation interval. Each meridional slice encodes one component of the α-tensor or combinations thereof. The 3 × 3 structure of the subfigures reflects the 3 × 3 components of the tensor, with the top-left to bottom-right diagonal corresponding to the diagonal components α_rr, α_θθ, and α_ϕϕ. The three meridional slices above this diagonal (top right) correspond to the off-diagonal elements of the symmetric part of the α-tensor (indices within parentheses), while the three plots below the diagonal represent the three independent components of the antisymmetric part of the tensor, plotted here as components of the turbulent pumping velocity ${\bm \gamma }$ defined through Equation (25). To facilitate comparison, the same color scale is used on all meridional slices.

Zoom In Zoom Out Reset image size

Although all components of the α-tensor are roughly of the same order of magnitude, the strongest components are found to be α_rr and α_ϕϕ, both antisymmetric about the equatorial plane. The former is highly structured spatially, peaking in subsurface layers and in the equatorial regions and with numerous sign changes in each hemisphere. The latter is spatially more homogeneous, with sign changes only across the equatorial plane and across a near-spherical surface located above the core–envelope interface, at r/R ≃ 0.76. Significant amplitudes are obtained in the off-diagonal contributions, particularly α_(rθ) and its associated pumping velocity γ_ϕ. This can be traced to the development of persistent non-axisymmetric flow structures in the low-latitude portions of the outer convection zone (see Figure 4), which contribute a strong signal to the "small-scale" flow component even after the azimuthally averaged mean-flow is subtracted (viz., Equation (11)). These flow structures are also responsible for the strong low-latitude signal seen in α_rr.

The α_ϕϕ component is of particular interest here, as it is the primary contributor to the production of the large-scale poloidal magnetic component (cf. bottom panel on Figure 2). Its first important structural property is antisymmetry with respect to the equatorial plane, as expected for cyclonic turbulence, where reflectional symmetry is broken by Coriolis forces (Parker 1955). Note that the fitting method we use to measure the components of the α-tensor treats each hemisphere independently, so that the high degree of antisymmetry observed here is a real characteristic of the simulation, rather than a fitting artifact. In the northern hemisphere, the α_ϕϕ component is positive in the bulk of the convection zone and peaks at high latitudes, as expected in view of the sense of twist imparted by Coriolis forces on diverging convective updrafts and converging downdrafts (Parker 1955). The sign change near the base of the convecting layer has been observed before in measurements of the α-tensor in other MHD simulations of convection (e.g., Ossendrijver et al. 2001) and is associated with the rapid decrease of the turbulent intensity as one moves downward toward and into the convectively stable fluid layer underlying the convection zone.

The SVD method has the additional advantage to automatically provide the variance in the estimate of a given parameter a_k. Denoting this variance by σ²_k, it is given by

$$\begin{equation} \sigma ^2_k = \sum _{i=1}^3 \left(\frac{V_{ki}}{w_{ii}}\right)^2\,. \end{equation} \tag{ 33 }$$

The standard deviation (square root of variance) is in the range of 1–2 m s⁻¹ in most of the convective layers for all α-components, but reaches much larger values in the underlying stable layers. This is inconsequential, as it arises simply because the EMF and large-scale magnetic field are both very weak below the convective zone, reducing greatly the quality of the fit in that region.⁵ Examination of the standard deviation distribution in meridional slices indicates that even though the stronger α-tensor components peak at high (α_ϕϕ) or low (α_rr) latitudes, the absolute quality of the fit is actually best at mid latitudes, where the EMF signal is strongest (see Figure 8), as expected. Nonetheless, at mid-convection zone depth the inferred values of α_ϕϕ in polar regions deviate from zero by more than three standard deviations, leaving no doubt that these values are physically meaningful. On the other hand, in the equatorial region delimited by |θ| ≲ 30° and 0.7 ≲ r/R_☉ ≲ 0.8 the standard deviation is of the same order of magnitude as α_ϕϕ, indicating a poorer fit in that area. This is simply due to the fact that both the EMF and the toroidal large-scale field are small in that region, as can be seen in Figures 2 and 8.

While the variance of the fit provided by the SVD is a useful, well-defined tool to characterize the accuracy of the fit, other methods can provide complementary information regarding the parameterization of the EMF. Here we demonstrate the quality of the fit by looking at the ability of the α-effect parameterization to reproduce the observed EMF from the α-tensor and the large-scale magnetic field. In Figure 10(A), we plot a time–latitude diagram of the residual between the observed EMF and the reconstructed EMF, i.e., the quantity

$$\begin{equation} \Delta _\phi = \mathcal {E}_\phi - \alpha _{\phi j}\langle B_j \rangle, \end{equation} \tag{ 34 }$$

computed at depth r/R_☉ = 0.85. It is quite remarkable to see that essentially no cyclic features remain in Δ_ϕ. Apart from the equatorial band where the EMF barely manifests itself, there are in fact no clearly discernible structures in Δ_ϕ. This implies that the α-effect parameterization performs very well at capturing the origin of the cyclic features of the EMF, the remainder of the EMF being essentially turbulent noise. This visual impression is confirmed by averaging latitudinally and performing a Fourier transform of the resulting time series to produce the power spectrum plotted in Figure 10(B). The low wavenumber portion of this spectrum is tolerably well represented by a 1/f power law and does not show any particular features or changes in slope at frequencies corresponding to the cycle period of the large-scale field or its first few harmonics.

Zoom In Zoom Out Reset image size

The three meridional slices under the diagonal in Figure 9 represent the components of the turbulent pumping speed, as defined through Equation (25). The vertical component γ_r is found to be negative through most of the convective envelope, with upward pumping taking place only in the subsurface layers. This is better viewed in Figure 11 showing radial cuts of γ_r extracted at latitude +45° (panel A, solid line) and latitudinal cuts extracted at depth r/R_☉ = 0.75 (panel B). This is qualitatively similar to the results obtained by Ossendrijver et al. (2002, see their Figure 5) in their local Cartesian simulations including rotation.⁶ This downward pumping reaches significant speeds, ∼1 m s⁻¹ in the bottom half of the convective layers, except at very low latitudes where it all but vanishes. Upward pumping occurs at all latitudes for r/R_☉ ⩾ 0.9, except at very low latitudes |θ| ⩽ 15°), where the γ_r > 0 region reaches almost to the base of the convecting layer (see Figure 9, central bottom meridional slice).

Zoom In Zoom Out Reset image size

We also detect in our simulation a significant equatorward latitudinal pumping in the outer two thirds of the convective envelope (see the bottom-left meridional slice in Figure 9 and dashed lines in Figures 11(A) and (B)), particularly prominent at mid latitudes where it reaches ∼1 m s⁻¹ in the middle of the convection zone. Both the direction and magnitude of this latitudinal pumping velocity are similar to those measured in the local simulations of Ossendrijver et al. (2002). Again, like these authors, we also observe significant poleward latitudinal pumping in the subsurface layers of the simulations, with speeds approaching 2 m s⁻¹ at latitude |θ| ≃ 60°. Given the ∼30 year duration of our cycles, it is therefore quite possible that the poleward drift and intensification of the surface magnetic field seen on the bottom panel of Figure 2 is driven at least in part by turbulent pumping, as proposed already by Ossendrijver et al. (2002).

It is also interesting to compare the turbulent pumping velocity to the drift speed of the large-scale toroidal magnetic field visible on the top and middle panels of Figure 2. The net drift speed of large-scale magnetic structures is influenced by other physical processes, notably turbulent diffusion. Indeed, if turbulent pumping plays a role in this drift its speed should still have values of the same overall order of magnitude as to the observed drift speeds. For most cycles, the equatorial drift of the toroidal component at the core–envelope interface (Figure 2, top panel), as measured by tracking the latitude of peak mean toroidal magnetic field at any given time, spans 25°–30° in 20–25 years, depending on individual cycles. This yields drift speeds in the range 0.3–0.5 m s⁻¹. Both the magnitude and direction of this drift compare well to the mid-latitude latitudinal turbulent pumping speed |γ_θ| ≃ 0.2 m s⁻¹ at r/R_☉ = 0.718, as plotted in Figure 11(A). Likewise, measuring the slope from r/R_☉ = 0.8 to 0.7 of the polarity reversal line on the middle panel of Figure 2 leads to a radial drift speed ∼−0.1 m s⁻¹, only a factor of five smaller than the mid-latitude values of γ_r in this range of depth, as plotted in Figure 11(A). All of this suggests that turbulent pumping may well play an important role in the spatiotemporal evolution of the large-scale, mean magnetic field building up in the simulation.

We next turn our attention to the possible time dependence of the tensor α_ij. Our entire analysis so far relies on the assumption that α_ij depends only on position and not time. As we have just shown, this hypothesis works fairly well in practice, in the sense that it provides an accurate parameterization of the EMF in terms of the large-scale magnetic field. At a more physical level, there are however reasons to believe that the α-effect could vary in time via a nonlinear dependence on the magnetic field strength.

In the solar context, the α-effect can arise via the systematic twist imparted by cyclonic convection on a pre-existing large-scale magnetic field. This idea was introduced by Parker (1955), who argued that it could circumvent Cowling's theorem by providing a viable regeneration mechanism for the Sun's large-scale poloidal magnetic component. Because magnetic tension should resist any twisting by the flow, it has been argued that Parker's mechanism can only operate up to the point where the large-scale magnetic field reaches local equipartition with the turbulent fluid motions. This has led to the introduction of a number of simple so-called α-quenching algebraic formulae, e.g.,

$$\begin{equation} \alpha _{ij} = {\alpha _{ij}^0\over 1+(\langle B_j\rangle /B_{\rm eq})^2}, \end{equation} \tag{ 35 }$$

where B_eq is the equipartition magnetic field strength (see also Rüdiger & Kitchatinov 1993).

In order to explore if the α-effect varies from one cycle to the next, in a manner related to the strength of the mean magnetic field, we refine our analysis by breaking up the simulation data into blocks each spanning 100 solar days (∼8 years). One set of blocks is centered on times of cycle maximum, when the large-scale magnetic field reaches peak strength, and a second set of interweaved blocks is centered on times of polarity reversals, where the large-scale magnetic component reaches its minimum. Because of the regularity of the cycles in the simulations (see Figure 2), and good hemispheric synchrony, such a partition of the simulation can be unambiguously defined and readily carried out. We then apply the fitting procedure described previously to each block of data independently, which yields a sequence of functions {α_ij(r, θ, t_c)}, where t_c is the time around which a given block is centered.

With the SVD fit now carried out over a much reduced time span, the resulting α-tensor is much noisier, with the standard deviation associated with the fit now often exceeding the fitted values themselves. In order to extract a statistically meaningful signal from these results, we spatially average the α-tensor components over meridional planes for each t_c, yielding a single measure α_c for each time block; for example, applying this procedure to the α_ϕϕ component

$$\begin{eqnarray} \alpha _c &=& \frac{1}{S} \int _{r_{{\rm min}}}^{r_{{\rm max}}} \int _{hem.} {\alpha _{\phi \phi }(r,\theta,t_c)\over \sigma _{\phi \phi }(r,\theta,t_c)} \, r \, {d}r\,{d}\theta \,, \end{eqnarray} \tag{ 36a }$$

$$\begin{eqnarray} S &=& \int _{r_{{\rm min}}}^{r_{{\rm max}}} \int _{hem.} {1\over \sigma _{\phi \phi }(r,\theta,t_c)} \, r \, {d}r\,{d}\theta \,. \end{eqnarray} \tag{ 36b }$$

Note that we use the standard deviation returned by the SVD algorithm as a weighting factor on the values being averaged. The radial integration bounds r_min and r_max are chosen to include the convective zone, but not the stable layers (r/R_☉ ⩽ 0.718) nor the surface of the simulation, in order to capture the contributions from regions where the α-effect is expected to operate substantially. The latitudinal integration is restricted to either the north or south hemisphere, to avoid uninteresting cancellations as the sign of α_ϕϕ flips from one hemisphere to the other. The standard deviation on the spatially averaged quantity is computed as

$$\begin{equation} \sigma _c = A/S, \end{equation} \tag{ 37 }$$

where A is the surface over which the above averaging integrals are defined. This amounts to treating the integral as a weighted sum of statistically independent measurements, probably a debatable hypothesis here, but nonetheless sufficient for the foregoing analysis.

The results of those calculations, as applied to the α_ϕϕ component, are shown in Figure 12, where each histogram-type column represents one temporal block. The total width of the error bars is given by 2σ_c. Viewed in this way, the actual values of α_c do not vary substantially in both hemispheres throughout the entire simulation and do not appear to show any clear, systematic variations between maximum and minimum phases. Similar results are obtained for the other components of the α-tensor. The α_ϕϕ component does show a very weak (r = −0.27) anticorrelation with the mean magnetic field strength spatially averaged in each time block, but even there the range of variation is only about one-half of the standard deviation.

Zoom In Zoom Out Reset image size

In an attempt to further suppress the noise, we then used the SVD fitting algorithm to do a single fit on all time blocks associated with maximum phases, and a separate fit for all time blocks centered on minimum phases. The resulting α-tensor components were then averaged spatially in the same manner as previously. The resulting α_c values, again for the ϕϕ component, are listed in Table 1 with the corresponding spatial average for the α_ϕϕ component extracted from the full simulation run, as plotted in Figure 9 (bottom right). Here we finally obtain a min/max difference that exceeds the range of the estimated error bars, and which moreover runs in the direction expected from the idea of α-quenching, namely a weaker α_c when the mean field is strongest. However, the associated variation in α_c between maximum and minimum phases remains rather modest, at the 20% level for the southern hemisphere, down to the 10% level in the northern.

The results of Table 1 reflect the variation of the spatially integrated α-profiles and could of course hide more substantial local variations. This is examined in Figure 13, showing now the spatial distribution of the α_ϕϕ component extracted from the set of concatenated data blocks centered on cycle minima (top row) and maxima (bottom) row, together with the corresponding spatial distributions of the standard deviation returned by our SVD fitting scheme (at right). At first glance both distributions appear quite similar, but a closer look reveals a marked reduction of α_ϕϕ at the base of the convection zone at maximum epochs, as well as at high latitudes. While these amplitudes of these variations (by ∼1–2 m s⁻¹) are well within the local one standard deviation at high latitude (∼2–4 m s⁻¹), at lower latitude they are comparable or slightly larger than one standard deviation, suggesting a statistically significant local variation. Note also that this is also where the toroidal large-scale component reaches its peak value (see Figure 2), so this is in fact where α-quenching would be expected to be most important.

Zoom In Zoom Out Reset image size

In contrast to the α_ϕϕ component, α_rr, α_θθ as well as the radial pumping velocity γ_r shows no significant local variation whatsoever between maximum and minimum epochs. Other components show high-latitude variations, but all at or within the associated one standard deviation range. In fact, the only other noteworthy pattern is a significant enhancement of the poleward latitudinal turbulent pumping in the surface layers at maximum epochs, by a factor of nearly two, which is slightly larger than the local one standard deviation range for γ_θ.

Taken as a whole, the results presented above indicate that little or no significant α-quenching occurs in our simulation, i.e., the magnitude of the α-tensor components is for the most part insensitive to the value of the large-scale magnetic field. This justifies, a posteriori, our use of the tensorial expansion of the turbulent EMF in terms of the mean magnetic field (viz., Equation (22)), as the starting point of the analysis carried out throughout this section.

To which (if any) of these three possible mean-field model categories could we identify our simulation? The answer clearly hinges on the relative importance of the $\langle \bm {u}\rangle \times \langle \bm {B}\rangle$ and $\langle \bm {u}^\prime \times \bm {B}^\prime \rangle$ terms on the right-hand side of the mean-field induction Equation (20). In keeping with our "experimental" approach, we simply use the mean and fluctuating flow and magnetic field components defined previously do directly calculate, using simple centered finite differences, the azimuthal components of both $\nabla \times (\langle \bm {u}\rangle \times \langle \bm {B}\rangle)$ and $\nabla \times \langle \bm {u}^\prime \times \bm {B}^\prime \rangle$, which then measures the rate of change of the magnetic field associated with each inductive contribution.

Figure 14 shows the results of this exercise, in the form of time–latitude diagrams for the azimuthal component of each of these two contributions (top and middle rows), as well as their sum (bottom row), extracted at the base of the convection zone (r/R_☉ = 0.718, left column) and higher up in the convection zone (r/R_☉ = 0.803, right column). The time span covered corresponds to the second half of the simulation run. It must be remembered that computed in this manner, the plotted tendencies reflect not only the action of source terms but also transport by turbulent pumping and large-scale flows.

Zoom In Zoom Out Reset image size

At any phase of the cycle, both contributions are here clearly of the same overall order of magnitude. Furthermore, constructing similar diagrams at other depths indicates that this remains the case throughout the bulk of the convection zone. In the terminology of mean-field electrodynamics, one would then conclude that in this one specific simulation, the large-scale magnetic field is sustained by an α²Ω dynamo, in the sense that both source terms on the right-hand side of Equation (41) contribute to the evolution of the toroidal field. However, it is also clear from Figure 14 that these two contributions to toroidal field production act in opposition to each other at most latitudes during most phases of the cycle. Interestingly, a similar pattern was observed in the MHD simulations of Brown et al. (2010), which produced steady large-scale magnetic fields. In their simulations D3 (see their Section 5.1), production of toroidal magnetic fields by the large-scale differential rotation shear was found to dominate over an opposing turbulent induction and Ohmic dissipation. Here, in the absence of explicit Ohmic dissipation and with both contributions to toroidal field production of comparable strengths, sign notwithstanding, the net tendency (bottom row) shows a spatiotemporal evolution that is markedly different from either contribution, and overall much smaller in magnitude.

The solid contours in Figures 14(C) and (F) show the spatiotemporal evolution of the large-scale toroidal magnetic component at the corresponding depths. At the base of the envelope the 〈B_ϕ〉 component and its temporal variation reconstructed from the upper two panels, shown as color levels, are seen to be offset in time, with d〈B_ϕ〉/dt leading 〈B_ϕ〉 in a manner consistent with the development of the observed polarity reversals at that depth. This phase offset, in turn, results from the residual tendency left from the near cancellation of the turbulent EMF and large-scale shear contributions plotted in panels A, B, and D, E. A similar temporal offset is observed at mid-convection zone depth at the mid latitudes, but not in the vicinity of the equatorial plane.

These results, especially when considered jointly with those of Brown et al. (2010, 2011), suggest that the development of cycles arises through a rather delicate balance between the magnitude and relative phasing of turbulent induction and large-scale shear. One could then legitimately suspect that the cycle period in these simulations is also sensitively dependent on this balance, which is itself likely dependent on various simulation parameters such as rotation rate, forcing, stratification, and so on. How and when the transition from steady to cyclic dynamo action takes place requires a detailed parameter study, which goes well beyond the scope of the present paper. Nonetheless, the demonstrated difficulty in "finding" regular cyclic large-scale dynamo action in these kinds of turbulent simulations, at the solar rotation rate at least, may well be a reflection of the need for such a fine balance to materialize.

The mathematical machinery of mean-field electrodynamics also involves the characterization of the turbulent EMF in terms of the mean magnetic field. In particular, for turbulent flows satisfying certain statistical properties and operating in certain specific physical regimes, it offers the mean of calculating the form of the α-tensor (see, e.g., Moffatt 1978). In the case of a turbulent flow that is isotropic and homogeneous but lacks reflectional symmetry, the α-tensor reduces to the diagonal form αδ_jk, with

$$\begin{equation} \alpha =-{\tau \over 3} h_v, \end{equation} \tag{ 42 }$$

where τ is the correlation time of the turbulence and h_v is the mean kinetic helicity associated with this turbulent flow:

$$\begin{equation} h_v=\langle \bm {u}^\prime \cdot \nabla \times \bm {u}^\prime \rangle. \end{equation} \tag{ 43 }$$

The key result here is that albeit for the sign difference, the α-effect coefficient is directly proportional to the kinetic helicity, a quantity readily computed a posteriori from the simulation output. The result of such a calculation is plotted in Figure 15(C), with the kinetic helicity having been averaged temporally as well as azimuthally over the complete time span of the simulation. Except for the sign difference, the degree of resemblance with the α_ϕϕ tensor component plotted on panel A is quite striking. Both functionals are strongly peaked at very high latitudes, and in each hemisphere only show a sign change near the base of the convection zone. The most prominent difference is found at low latitudes, where h_v shows a secondary peak in the middle of the convection zone that is absent in the α_ϕϕ profile. The α_θθ component (see Figure 9) is also tolerably well represented by the negative of the mean kinetic helicity, the most prominent discrepancy being that α_θθ peaks at mid latitudes. The α_rr diagonal component is that which shows the least similarity with the mean kinetic helicity, presumably because it is more sensitive to the break of homogeneity on which Equation (42) is predicated, induced by stratification of the background state and/or upper boundary condition. Another significant departure relates to the temporal behavior of the kinetic helicity, which shows a clear signature of the cycle in the large-scale magnetic field, while most α-tensor components extracted from the simulation do not, at least not at a level that could be deemed statistically significant; in the simulation, the cyclicity of the EMF is well captured by the cyclicity of the large-scale magnetic component (viz., Figure 10).

Zoom In Zoom Out Reset image size

Equation (42) is also predicated on a number of other strong assumptions, including the requirement that the small-scale turbulent flow remains entirely unaffected by the small-scale turbulent magnetic field. Once the latter becomes dynamically significant, it is expected that it will alter the small-scale flow giving rise to the α-effect. Based on a set of local Cartesian incompressible MHD simulations operating in the strongly nonlinear regime, Pouquet et al. (1976) suggested that Equation (43) should be replaced by

$$\begin{equation} \alpha ^\ast =-{\tau \over 3}\left(h_v-h_B\right), \end{equation} \tag{ 44 }$$

where h_B is the mean current helicity, generalized to our anelastic MHD formulation as

$$\begin{equation} h_B={1\over \mu _0\rho _0} \langle \bm {B}^\prime \cdot \nabla \times \bm {B}^\prime \rangle. \end{equation} \tag{ 45 }$$

The physical link between kinetic and current helicity becomes more transparent upon noting that ${{\bm B}^\prime }/\sqrt{\mu _0\rho _0}$ is the Alfvén velocity, measuring the propagation speed of transverse MHD waves for which magnetic tension provides the restoring force; in essence, what Equation (44) expresses is simply that magnetic tension tends to oppose any twisting of magnetic field lines by the small-scale turbulent flow. The mean current helicity in our simulation is plotted in Figure 15(D), again as a combined azimuthal and temporal average spanning the full simulation duration. This magnetic contribution to the α-coefficient is found to be significantly smaller than its kinetic counterpart, by a factor of about 10. The strong peaks straddling the base of the convection zone at mid latitude in both hemispheres are associated with the strong mean magnetic field building up in these regions, which, acted upon by the small-scale turbulent flow, feeds the production of small-scale current helicity.

Figure 15(B) shows the α*(r, θ) profile reconstructed according to Equation (44). The correlation time within the unstable layers is estimated in a similar manner as in Brown et al. (2010), i.e.,

$$\begin{equation} \tau (r)= {H_\rho \over u^\prime }, \end{equation} \tag{ 46 }$$

where H_ρ is the density scale height of the background stratification and u' is the rms average of the small-scale flow component, the averaging being carried out zonally, latitudinally, and temporally over the full time span of the simulation. In the stable layers, τ is artificially set to zero, since Equation (44) is not expected to hold there, pertaining as it does to fully developed turbulence.

With the kinetic helicity about an order of magnitude larger than the magnetic helicity, the structure of this estimated α*(r, θ) reflects primarily the spatial structure of the kinetic helicity, except in the immediate vicinity of the core–envelope interface where the kinetic and magnetic contributions have comparable magnitudes. While the overall amplitude of this reconstructed α profile is about a factor of three larger than the measured α_ϕϕ (Figure 15(A)), in terms of the spatial variations Figures 15(A) and (B) are remarkably similar. Interestingly, the overall magnitude of the reconstructed α, ∼20 m s⁻¹, is quite comparable to the α reconstructed similarly in Brown et al. (2010, see their Figure 8) in a simulation rotating at thrice the solar rate, and sustaining a temporally steady, rather than cyclic, large-scale magnetic field component.

The limited applicability of Equation (44) to our simulation notwithstanding, the weak contribution of the specific current helicity as compared to the kinetic helicity term suggests that the dynamical impact of the small-scale magnetic component on the small-scale flow is correspondingly weak. Recall that we also demonstrated, in Section 3.6, that the α-tensor components show little or no significant variation with the strength of the large-scale magnetic field. This would suggest that the α-effect in our simulation operates in the quasi-linear regime and that the saturation of the dynamo amplitude must be achieved via a different mechanism.