Where is the solar dynamo?

In this paper I review results from recent global magnetohydrodynamical numerical simulations of solar convection, as a springboard to address the question "Where is the solar dynamo". I first describe and contrast similarities and differences in the large-scale flows and magnetic fields such simulations can produce, with emphasis on polarity reversals (or lack thereof) in the large-scale magnetic components they generate. On the basis of these simulation results, I argue that some of the significant differences in the spatiotemporal evolution of the large-scale magnetic field can be traced to the competing effects of turbulent electromotive forces and induction by large-scale flows, whose mutual near-cancellation in the nonlinearly saturated regime leads to a high sensitivity to the numerical/physical treatment of small scales. Some of these recent simulation results also reopen the possibility that dynamo action driving the solar activity cycle may reside entirely within the convection zone, with the tachocline perhaps playing a lesser role than has been assumed in the last two decades. On the other hand, other subsets of simulations suggest that magnetohydrodynamical processes taking place within the tachocline may have a significant impact on timescales comparable to or longer than the primary cycle.


Where the solar dynamo was
Helioseismology has changed forever the way we go about modelling the solar dynamo, as it provides us with vital information on the flows pervading the solar interior. These flows are the dominant energy source for the magnetohydrodynamical (MHD) dynamo process powering all of solar activity. Interestingly, helioseismology's initial impact on solar dynamo modelling was to stir trouble in what was then emerging as an almost consensual view. Instead of rotational isocontours constant on cylinders matching the observed surface latitudinal differential rotation, as expected on theoretical grounds, the first helioseismic inversions of the deep solar rotation (Brown & Morrow 1987;Brown et al. 1989;Dziembowski et al. 1989) revealed a shear that remains largely latitudinal from the subsurface layers down to the base of the solar convection zone (0.71 ≤ r/R ⊙ ≤ 1), vanishing abruptly below its base across a thin shear layer now called the tachocline. At about the same time, serious doubts were beginning to be raised concerning the viability of the other cornerstone of the solar dynamo models developed in the 1970 and 1980, namely the turbulent electromotive force associated with cyclonic turbulence, as embodied in the so-called α-effect (see Cattaneo et al. 1996;Blackman & Brandenburg 2002; and references therein). Finally, the first pioneering global magnetohydrodynamical simulations of solar convection (Gilman 1983;Glatzmaier 1984Glatzmaier , 1985 also challenged many prevailing theoretical tenets of the solar cycle as modeled through the dynamically-and geometricallysimplified dynamo models relying on mean-field electrodynamics. This state of affair led to the development of dynamo models specifically designed to bypass some of these difficulties, e.g., the interface dynamo picture proposed by Parker (1993), as well as the revival of other classes of solar cycle models formerly eclipsed by the rise of meanfield electrodynamics. Most notable among these, perhaps, are the so-called Babcock-Leighton models (Wang & Sheeley 1991;Choudhuri et al. 1995;Durney 1995;Dikpati & Charbonneau 1999), in which the sun's large-scale poloidal field is regenerated through the decay of active regions. The study of the stability and rise of thin magnetic flux tubes, presumably producing sunspots upon emergence through the photosphere, also got underway and produced a number of key insights, most important perhaps the notion that such sunspot-producing magnetic flux ropes could only be stored in the weakly subadiabatic overshoot layer underlying the solar convection zone (Ferriz-Mas et al. 1994). The strong rotational shear therein was also seen as conducive to the production of intense toroidal magnetic fields, from which strongly magnetized flux ropes would presumably form. Various instabilities were later identified, which could drive the production of a poloidal magnetic component from a strong tachocline toroidal magnetic field, thus opening the possibility of closing the dynamo loop wholly within the tachocline. For more details on all these solar cycle models, see Charbonneau (2010).
This very brief historical survey, incomplete as it may be, nonetheless illustrates the main point I want to make here: as recently as a few years ago, if one were to be asked "where is the solar dynamo", the short answer would have likely been along the line "mostly in the tachocline". But remarkably, just as this almost consensual new picture was emerging, something else came around to stir trouble once again.

Magnetic cycles in MHD simulations of solar convection
It was clearly in the air; in the years 2009-2011, various groups working independently and using distinct numerical and modelling approaches managed to achieve turbulent MHD simulations of solar/stellar convection producing magnetic fields well-organized on large spatial scales, and undergoing polarity reversals which, in some cases, showed remarkable similarities with the solar cycle. In this section I review three such classes of numerical simulations, highlighting similarities and differences in the flows and magnetic fields they generate. In principle, all simulate the same thing: turbulent convection and MHD dynamo action in a thick, rotating, stratified, convectively unstable shell of electricaly conducting fluid subjected to thermal forcing sustaining convection.
The PENCIL code (www.nordita.org/pencil-code) is a public domain parallel high order finite-difference compressible MHD code, formulated in cartesian geometry. PENCIL uses constant values for dissipative coefficients such as viscosity, thermal and magnetic diffusivities, and so operates as a Direct Numerical Simulation (DNS). Pseudo-global simulations of magnetic cycles in a zonally periodic spherical wedge have been presented by Käpylä et al. (2010Käpylä et al. ( , 2012.
The ASH code (Anelastic Spherical Harmonics; see Clune et al. 1999) is a massively parallel spectral code in spherical geometry. Most published solar/stellar simulations (e.g., Brun et al. 2004;Browning et al. 2006;Brown et al. 2010Brown et al. , 2011, and references therein) use a Large-Eddy Simulation (LES) approach, where nonlinear stability is maintained through the use of enhanced values for dissipative coefficients such as viscosity and magnetic diffusivity. Some of the most recent ASH simulations (Nelson et al. 2011(Nelson et al. , 2013a(Nelson et al. , 2013b, including some discussed further below, make use of a dynamical subgrid model which yields significantly lower dissipation levels than prior ASH LES. The EULAG code (EUlerian-LAGrangian; www.mmm.ucar.edu/eulag; see also Prusa et al. 2008) is a publicly available, robust hydrodynamical flow solver formulated in generalized curvilinear coordinates. Its MHD extension (not yet publicly available) is a relatively recent development and is documented in Smolarkiewicz & Charbonneau (2013). EULAG achieves strongly turbulent, low-dissipation flow states through the use of a non-oscillatory forwardin time advection algorithm, and thus belongs to the class of so-called Implicit Large-Eddy Simulations (ILES), where dissipation is delegated to the numerical scheme, rather than being formulated explicitly as in LES. Global simulations of solar-like cycles in the anelastic approximation have been reported in Ghizaru et al. (2010) and Racine et al. (2011). Despite widely varying numerical approaches, presence or absence of an underlying stable fluid layer, treatment of small scales, resolution, and physical parameter regimes, the convective flows produced by these simulations, and the large-scale flows they power via inverse cascade, show remarkable similarities. Convective patterns take the form of broad upflows of hotter fluid, delineated by a network of narrow downflow lanes of colder fluid, typical of thermal convection in a gravitationally stratified background. This pattern is broken only in the equatorial regions, where larger, latitudinally-elongated convective cells are stacked longitudinally, with the flow showing little variations in the direction parallel to the rotation axis, as expected in rotationallyconstrained thermal convection (cf. Fig. 2 in Käpylä et al. 2010, Fig. 1A,B,C in Nelson et al. 2013a; Fig. 1 in Ghizaru et al. 2010;§2.2.3 in Miesch & Toomre 2009). This rotational structuring of convection is also reflected in the large-scale zonal flow, in all cases showing equatorial acceleration but with a strong tendency for cylindrical isocontours throughout the convection zone (cf. Fig. 1A in Käpylä et al. 2010, Fig. 7A in Nelson et al. 2012, andFig. 4 in Racine et al. 2011), unlike the primarily radial isocontours revealed by helioseismology (e.g., Howe 2009, Fig. 19). Moreover, in all three cases considered here the pole-to-equator angular velocity contrast turns out smaller than observed by a factor of about 3; since purely hydrodynamical but otherwise identical parent simulations tend to yield a contrast much closer to observed, this indicates that magnetic forces are playing an important role in zonal dynamics.
Few of these simulations yield a large-scale meridional flow profile (the axisymmetric flow component in [r, θ] planes) that resembles the smooth, steady single flow cell per quadrant commonly used in flux transport dynamos. Instead, the zonally-averaged meridional flow is typically very weak, highly time variable, and characterized by multiple cells in both radius and latitude (cf. Fig. 1C in Käpylä et al. 2012, Fig. 1D in Brown et al. 2010, andFig. 4 in Racine et al. 2011).
Subtracting these axisymmetric large-scale flows extracted from the simulations yields the "small-scale" contribution u ′ , from which the zonally-averaged kinetic helicity h u = u ′ · ∇ × u ′ , measuring the cyclonic character of the turbulence, is readily computed. Once again the kinetic helicity profiles turn out remarkably similar, with h u negative (positive) in the Northern (Southern) hemisphere, peaking at high latitudes with a secondary peak in equatorial regions, often with a sign change near the base of the convection zone (cf. Fig. 1B in Käpylä et al. 2012, Fig. 8A in Brown et al. 2010, and Fig. 15C in Racine et al. 2011. With such striking similarities in the flows powering dynamo action, one might rightfully expect that the large-scale magnetic fields produced by their associated dynamo action should show a comparable level of similarities. Interestingly, this is not the case, although some similarities do exist. Most prominent perhaps is the production of persistent and spatially wellorganized toroidal magnetic flux structures in equatorial regions, undergoing polarity reversals on timescales of a few years, in a manner more or less regular and synchronous across hemispheres, depending on the simulation considered. Figure 1A shows an example taken from a EULAG simulation, in the form of a Mollweide projection of the toroidal magnetic field extracted at r/R = 0.88. The slightly tilted banded structures at low latitudes, maintaining their polarity over a significant longitudinal extent, are also seen in other MHD simulations (cf., e.g., Fig. 2 in Brown et al. 2010;Fig. 4 in Nelson et al 2013a).
Interestingly, zonal averaging of this non-axisymmetric magnetic field still yields a significant axisymmetric component, as shown on Figure   which is what mean-field theory would predict for a kinetic helicity negative in the Northern hemisphere. Using the full α-tensor extracted by Racine et al. (2011) from a similar EULAG simulation, Simard et al. (2013) constructed kinematic α 2 Ω mean-field dynamo models that also exhibit an equatorially-concentrated, poleward-propagating dynamo mode, supportingwithout strictly proving-a dynamo wave interpretation. Cycles on longer timescales also materialize in some of these MHD simulations, when a stably stratified fluid layer is included below the convecting fluid layers. Figure 3 shows an example, taken from a low resolution simulation segment spanning two and a half centuries of simulated time. Time-latitude and time-radius diagrams of the zonally-averaged toroidal magnetic component are shown, in same format as on Fig. 2. Here regular polarity reversals, essentially synchronous in both hemispheres, take place every ∼ 40 yr. Notice how the largescale toroidal field now peaks at mid-latitudes and deep down, slightly below the base of the convecting layers.
Close examination of this simulation reveals a short period oscillation present at low latitudes within the convection zone (see Fig. 4A below). This appears to be the faster dynamo mode already encountered on Fig. 2, superimposed on the deep-seated, dominant mode and, as a consequence, enslaved to the latter's polarity reversals. Observational evidence for such "double-cycle" behavior actually abounds in indicators of solar activity (Benevolenskaya 1995), geomagnetic activity (Mursula et al. 2003), as well as in solar-cycle related shifts in p-mode frequency splittings (Fletcher et al. 2010; and references therein).
The ILES approach underlying these EULAG simulations makes it possible to reach strongly turbulent states on relatively coarse spatial grids 1 , which in turns allows very long temporal 1 A persistent and widely held belief is that low spatial resolution automatically implies high dissipation. It is certainly the case that for a given numerical scheme and/or subgrid scale model, dissipation scales with resolution; but it is also the case that different treatment of subgrid scales can yield different effective dissipation at a given mesh size. The set of ASH simulations recently published by Nelson et al. (2013a) offers a wonderful case in point. By switching from a static, eddy diffusivity-like LES scheme to a dynamical Smagorinsky model, these authors report a ∼ 50-fold decrease in dissipation on a fixed spatial mesh (cf. the quoted values for ν and η for simulations D3b and S3 in their Table 1). This is consistent with decades of community experience in computational fluid dynamics using ILES schemes; see the papers collected in Grinstein et al. (2007) for specific examples, including thermally-driven convection. integrations. Figure 4 shows the simulation of Fig. 3, now extended to 800 yr. The top plot is again a time-latitude diagram of the zonally-averaged toroidal field, now extracted from the inner third of the convecting fluid layer. The well-defined cycle settling in after initialisation, antisymmetric about the equatorial plane, undergoes a switch to a symmetric mode following the failure of the Southern hemisphere to inverse its polarity at around 390 yr. Some 150 years later the cycle seems to vanish, but this turns out to be an artefact of the zonal averaging used to produce this plot. The bottom plot is a time-longitude diagram of the toroidal field at (r/R, θ) = (0.82, 45 • ), and reveals that around 500 yr the dynamo switches now to a nonaxisymmetric, tilted-dipole-like configuration, which continues to undergo polarity reversals, with a slight reduction in both period and amplitude. Ongoing analysis of this simulation suggests that the destablization towards the non-axisymmetric m = 1 mode is mediated by the runaway growth of a MHD instability triggered in the upper portion of the stable fluid layers underlying the convection zone (on this topic see, e.g., Gilman & Fox 1997;Dikpati et al. 2004;Strugarek et al. 2011; and references therein).

Old casks for new wines ? Or is it the other way around ?
It has been well understood for quite some time now that turbulent convective flows of electrically conducting fluid can be powerful dynamos (e.g., Cattaneo 1999). With regards to the solar cycle, the real challenge is to identify the mechanism(s) through which the resulting magnetic fields, generated primarily on temporal and spatial scales commensurate with those of the convective flow, can self-organize into temporally persistent and spatially coherent structures on the scale of the sun as a whole. Traditionally, this self-organization has been viewed as a form of inverse cascade, similar to that powering large-scale flows via Reynolds stresses in purely hydrodynamical turbulence. This process can be formalized through mean-field electrodynamics in a manner such as to allow, under certain (physically restrictive) circumstances, the calculation of source terms for the large-scale magnetic field arising from the inductive action of small-scale flow and field. If a Eclipse on the Coral Sea: Cycle 24 Ascending (GONG 2012, LWS/SDO-5, and SOHO 27) IOP Publishing Journal of Physics: Conference Series 440 (2013) 012014 doi:10.1088/1742-6596/440/1/012014 good separation of scales exists between the large-scale magnetic component and the turbulent flow, such source terms end up embodied in the so-called α-effect, refering to the first term in the tensorial development of the mean turbulent electromotive force ξ ξ ξ ξ in terms of the mean magnetic field B: where α α α α is the α-tensor, whose calculation in terms of the statistical properties of the turbulent flow and field becomes the crux of the matter. Mathematically, representation of the turbulent electromotive force in terms of the large-scale magnetic field via eq. (1) is always possible, although it may not always be physically meaningful. Nonetheless, and with this caveat in mind, mean-field electrodynamics can also be turned around and used as an interpretive tool to analyze turbulent induction in MHD simulations. Accordingly, measurements of the turbulent emf and α-tensor components in turbulent MHD simulations have been carried out, and yielded results often in surprisingly good agreement with the predictions from mean-field theory, even though the simulations operate in physical regimes where many underlying assumptions under which the α-tensor can be computed analytically do not hold (see, e.g., Ossendrijver et al. 2001Ossendrijver et al. , 2002Käpylä et al. 2006Käpylä et al. , 2009Racine et al. 2011;but also Courvoisier et al. 2006). Chief among these restrictive assumptions is the requirement that the turbulence not be affected by the large-scale magnetic field, which is seldom the case in most simulations considered here (on this point see also Hubbard et al. 2009). In fact, in the context of large-scale dynamo action at high magnetic Reynolds number, the constraint of magnetic helicity conservation guarantees that the small and large scales will be coupled through a dual forward/inverse cascade of magnetic helicity of opposite sign (see Brandenburg & Subramanian 2005, and references therein). A different angle on the inevitable coupling of small and large scales in the low dissipation limit can be directly obtained from the MHD induction equation governing the evolution of the large-scale magnetic field: Assume now that because dissipation is very small, the last term on the RHS can be neglected. Then, for the mean field to remain quasi-stationary on the timecale of convection one must evidently have: implying the near-cancellation of the turbulent emf and induction by the large-scale flow.
Remarkably, this is precisely what is found in the EULAG simulation analyzed by Racine et al. (2011;see their Fig. 14): the "long" cycle is driven by the small residual total emf resulting from this near-cancellation of the two inductive contributions. A similar situation was also uncovered in the ASH simulations presented in Brown et al. (2010), which produced steady large-scale magnetic fields within the convecting layers. From the modeling point of view, the implications of this state of affairs are worrisome. Minor changes in the treatment of small scales can alter slightly the turbulent emf (say), but these can translate into large changes in the residual total emf, and therefore large changes in global behavior. Another exciting recent development has been the autonomous production of strongly magnetized (many times equipartition) flux rope-like structures in the ASH simulations of Nelson et al. (2011Nelson et al. ( , 2013aNelson et al. ( , 2013b, a process by all appearance favored by the strong turbulent intermittency characterizing these simulations. Figure 5 shows an example, taken from Nelson et al. (2013a;Fig. 18). Some flux strands building up within toroidal wreaths develop enough of a density deficit in their cores for magnetic buoyancy to contribute to their rise through the convecting layers, after being pulled away from the wreaths by strong convective updrafts. Subsequent detailed analysis by Nelson et al. (2013b) has shown that most of these structures maintain their approximately East-West orientation as they rise, and thus would comply with Hale's polarity rules upon emergence. This offers an intriguing alternative to the prevalent idea that sunspot-forming toroidal magnetic flux ropes must be formed, or at least stored, in the overshoot layer within the tachocline.

Coda: where is the solar dynamo?
Global three-dimensional MHD numerical simulations of solar-like magnetic cycles have landed.
But what have we learned, and what can we still hope to learn, from such simulations? Perhaps the most significant lesson learned is that the buildup of a large-scale magnetic component, retaining spatial coherence over timescales much longer than the convective turnover time, is possible in simulations operating in strongly turbulent regimes. Now, it was not at all obvious a priori that this should be possible. In the turbulent regimes characterizing MHD simulations operating at low dissipation, enhancements of both turbulent induction (stretching and/or constructive folding of magnetic fieldlines) and turbulent dissipation (destructive folding of magnetic fieldlines) are expected. Until recently, whether the former could outpace the latter remained very much an open question. The fact that three distinct classes of MHD simulations, relying on different physical setups and numerical approaches, should all provide the same positive answer is noteworthy indeed. The overall complexity of global dynamo action and magnetic self-organization in the strongly turbulent regime is also well established by these simulations; large-scale flows generated by convection itself are important players in the growth and evolution of the large-scale magnetic field. The turbulent emf is quite complex and attempts to relate it to the mean magnetic field yield an intricate picture where turbulent induction, turbulent pumping and perhaps even time lags all play a significant role. Further complications ensue in simulations including a stably stratified fluid layer underlying the convection zone. In the absence of magnetic fields, such simulations naturally develop a thin, tachocline-like shear layer immediately beneath this interface (e.g. Brun et al. 2011, Fig. 10;Racine et al. 2011, Fig . 3). This shear, although reduced in MHD simulations, can still play a significant role in global dynamo action, and in some cases lead to the buildup of distinct dynamo modes, peaking in different regions of the domain and interfering with one another.
What have these MHD simulations provided in terms of a (tentative) answer to the query raised in the title of this paper? The fact that spatially coherent, persistent magnetic flux Eclipse on the Coral Sea: Cycle 24 Ascending (GONG 2012, LWS/SDO-5, and SOHO 27) IOP Publishing Journal of Physics: Conference Series 440 (2013) 012014 doi:10.1088/1742-6596/440/1/012014 systems can be sustained within the convection zone, and even produce buoyantly rising flux rope-like structures, certainly brings back to the fore the possibility that the dynamo driving the solar cycle may "live" entirely within the convection zone. Käpylä et al. (2012) even produced simulations which, in their nonlinearly saturated regime, exhibit a very solar-like equatorward propagation. Then again, the EULAG simulations producing deep cycles in the tachocline also show a number of solar-like secondary features, such as double cycles (Simard et al. 2013), magnetically-mediated modulation of convective energy transport (Cossette et al., preprint), of the meridional flow (Passos et al. 2012) and rotational torsional oscillation patterns that compare very favorably to observations (Beaudoin et al. 2013). Where is the solar dynamo? We still don't know, but there is a fair chance that the little light we are currently seeing, very far away at the end of the tunnel, for once may not be the proverbial incoming locomotive.