Evidence of active MHD instability in EULAG-MHD simulations of solar convection
EEvidence of active MHD instability in EULAG-MHD simulationsof solar convection
Nicolas Lawson, Antoine Strugarek, and Paul Charbonneau
D´epartement de Physique, Universit´e de Montr´eal, C.P. 6128 Succ. Centre-ville, Montr´eal,Qc H3C 3J7, Canada [email protected]@[email protected]
ABSTRACT
We investigate the possible development of magnetohydrodynamical insta-bilities in the EULAG-MHD “millenium simulation” of Passos & Charbonneau(2014). This simulation sustains a large-scale magnetic cycle characterized bysolar-like polarity reversals taking place on a regular multidecadal cadence, andin which zonally-oriented bands of strong magnetic field accumulate below theconvective layers, in response to turbulent pumping from above in successive mag-netic half-cycles. Key aspects of this simulation include low numerical dissipationand a strongly subadiabatic fluid layer underlying the convectively unstable lay-ers corresponding to the modeled solar convection zone. These properties areconducive to the growth and development of two-dimensional instabilities oth-erwise suppressed by stronger dissipation. We find evidence for the action of anon-axisymmetric magnetoshear instability operating in the upper portions of thestably stratified fluid layers. We also investigate the possibility that the Taylerinstability may be contributing to the destabilization of the large-scale axisym-metric magnetic component at high latitudes. On the basis of our analyses, wepropose a global dynamo scenario whereby the magnetic cycle is driven primar-ily by turbulent dynamo action in the convecting layers, but MHD instabilitiesaccelerate the dissipation of the magnetic field pumped down into the overshootand stable layers, thus perhaps significantly influencing the magnetic cycle pe-riod. Support for this scenario is found in the distinct global dynamo behaviorsobserved in an otherwise identical EULAG-MHD simulations, using a differentdegree of subadiabaticity in the stable fluid layers underlying the convection zone.
Subject headings:
Instabilities — Magnetohydrodynamics (MHD) — Sun: activ-ity — Sun: magnetic field — Stars: activity — Stars: magnetic field a r X i v : . [ a s t r o - ph . S R ] S e p
1. Introduction
The dynamo-based model of the solar cycle put forth by Parker (1955a,b) over half acentury ago has stood the test of time remarkably well. The joint inductive action of dif-ferential rotation and helical turbulence remains at the heart of many contemporary solarcycle models, but helioseismic inversion of the sun’s internal differential rotation has broughtincreased attention to the tachocline, a rotational shear layer straddling the base of the solarconvective envelope (Howe 2009), as the locus of toroidal field amplification and storage,prior to its buoyant destabilization and rise to the photosphere to produce bipolar activeregions (Fan 2009). In these models, turbulent induction within the convection zone and/orgeneration of a surface dipole moment through the decay of active regions provide the regen-erative mechanism required to close the dynamo loop (see Charbonneau 2010, for a surveyof these various models). Moreover, various physical mechanisms have been identified, whichcould power a dynamo contained entirely within the tachocline. For example, Schmitt &Sch¨ussler (1989) (see also Ossendrijver 2000a,b) have argued that helical waves growing alongtoroidal flux tubes stored within the upper stably stratified portion of the tachocline couldprovide an azimuthal electromotive force able to regenerate the poloidal component in situ;likewise, Dikpati & Gilman (2001a) have shown that in the presence of rotation, the jointmagnetohydrodynamical (hereafter MHD) instability investigated by Gilman & Fox (1997),operating in the tachocline, develops a net hemispheric helicity that could provide an analogof the turbulent electromotive force proposed by Parker (1955b) and mathematically for-malized by Steenbeck & Krause (1969) as the “ α -effect” of mean-field electrodynamics (seeDikpati & Gilman 2001b, for an example of flux-transport dynamo models with an α -effectoriginating from such instability).The vast majority of the solar cycle models built using these various regenerative mag-netic field mechanisms operate in the so-called kinematic approximation, whereby the mag-netic back reaction on the inductive flows is altogether neglected or incorporated in themodels through largely ad hoc parameterizations. Global MHD simulations of solar convec-tion do not suffer from this shortcoming, but it is only recently that advances in computingpower and algorithmic design have jointly led to global simulations producing magnetic fieldswell-organized on global scales as well as undergoing (more or less) regular polarity rever-sals (e.g. Ghizaru et al. 2010; Brown et al. 2011; Racine et al. 2011; K¨apyl¨a et al. 2012;Masada et al. 2013; Nelson et al. 2013; Fan & Fang 2014; Passos & Charbonneau 2014;Augustson et al. 2015). However, few of these simulations include a stably stratified fluidlayer underlying the convection zone, and those which do often use strongly enhanced dissi-pative coefficients to ensure numerical stability, which leads to dissipative dynamics in theconvectively stable layers. 3 –The EULAG-MHD simulations reported upon in Ghizaru et al. (2010) (Racine et al.2011; Passos & Charbonneau 2014) offer an interesting exception. In these simulationsnumerical stability is enforced via the advection algorithm itself, which effectively providesan adaptive subgrid model introducing only the minimal amount of dissipation required tomaintain stability in regions of strong shear in the flow or magnetic field, and very littlein smooth regions (see, e.g., Domaradzki et al. 2003). Such simulations thus offer a uniqueopportunity to investigate dynamical effects taking place in the stably stratified layers, inparticular the occurrence of instabilities otherwise suppressed by strong dissipation. Thisis the primary aim of this paper. Working with the EULAG-MHD “millenium simulation”presented in Passos & Charbonneau (2014) and briefly described in §
2, we first investigatein § §
4, following Miesch (2007) we then seek evidence for the development of amagnetoshear instability in the stable layer of the simulation and in § §
2. Simulation characteristics
The foregoing analyses are based on the EULAG-MHD “millenium simulation” pre-sented and analyzed by Passos & Charbonneau (2014). This simulation is based on thenumerical solution of the anelastic magnetohydrodynamical equations in a thick, gravita-tionally stratified shell of electrically conducting fluid, rotating at the solar rate, and sub-jected to thermal forcing driving convection (see Ghizaru et al. 2010; Racine et al. 2011).The simulation is performed using EULAG-MHD (Smolarkiewicz & Charbonneau 2013),a MHD generalization of the robust multi-scale geophysical flow solver EULAG (Prusa &Smolarkiewicz 2003; Prusa et al. 2008). We operate EULAG-MHD in its so-called ImplicitLarge-Eddy Simulation mode, whereby the dissipation required to maintain numerical sta-bility is delegated to the underlying advection scheme, which in this case is analogous to anadaptive subgrid model where the minimal level of dissipation required to maintain stabilityis introduced only where and when it is required. This allows to reach turbulent regimes onrelatively small spatial meshes, in turn allowing temporally extended integrations. This isparticularly important in the solar cycle context, considering the vast disparity of timescalebetween the convective turnover time (hours to days in the outer reaches of our solutiondomain), and the large-scale magnetic cycle, with its multi-decadal period.The millenium simulation used in what follows spans 1600 years, in the course of which39 polarity reversals take place. The cycles are quite regular, with a mean period 40 . ± . . ≤ r/R ≤ .
96, discretized on a modest spatial mesh 128 × ×
48 in longitude, latitude andradius. The background stratification is convectively stable below r/R = 0 . Figure 1 illustrates some characteristics of the magnetic field building up in our sim-ulation. Panel A and B show, respectively, a Mollweide projection of the radial flow andmagnetic field components on a spherical shell at r/R = 0 . ∼
40 yrhalf-period, antisymmetry about the equatorial plane, and good hemispheric synchrony.Here, and in what follows, zonal averages are indicated by angular brackets and com-puted directly from the simulation output as, e.g., (cid:104) B φ (cid:105) ( r, θ, t ) = 12 π (cid:90) π B φ ( r, θ, φ, t )d φ . (1)What we will refer to as non-axisymmetric components, denoted by primes, are then obtainedby subtracting such zonal averages from the total flow ( u ) and magnetic field ( B ) producedby the simulation: u (cid:48) ( r, θ, φ, t ) = u ( r, θ, φ, t ) − (cid:104) u (cid:105) ( r, θ, t ) , B (cid:48) ( r, θ, φ, t ) = B ( r, θ, φ, t ) − (cid:104) B (cid:105) ( r, θ, t ) , (2)Figure 2 shows yet another view of the same simulation data, this time in the form of 5 –Fig. 1.— Snapshots in Mollweide projection of (A) the radial flow component, (B) the radialmagnetic component, and (C) toroidal magnetic component, as produced by the the EULAG-MHD “millenium simulation” used in the foregoing analysis. Panels A and B are constructedon a spherical shell at mid-depth in the convection zone ( r/R = 0 . r/R = 0 . r/R = 0 . (cid:39) . (cid:104) B φ (cid:105) ) φ max ( (cid:104) B θ (cid:105) ) φ + max ( (cid:104) B r (cid:105) ) φ , (3)which reaches values in the range [0 . ,
1] at cycle maxima. A well-organized toroidal com-ponent is also present throughout the bulk of the convecting layers, although with weakeramplitude, reaching ∼ . T in mid-latitudes at cycle maximum.The differential rotation shown on Figure 2B is solar-like, in that it is characterized bysignificant equatorial acceleration and polar deceleration, both vanishing in the stable layeracross a thin shear layer straddling the base of the convective envelope. In this simulation,the pole-to-equator contrast in angular velocity is actually too small by a factor of ∼ r, θ ) plane representations of (A) the zonally-averaged toroidal (colorscale) and poloidal (black and white contours, field lines oriented clockwise and anticlockwiserespectively) magnetic components at a cycle maximum, (B) the rotational angular velocityaveraged zonally and temporally over the full length of the simulation, and (C) a snapshotof the radial flow component at a fixed longitude. The black vertical line is the rotation axis,and the green line indicates the base of the convective zone, at r/R = 0 . In the EULAG-MHD millenium simulation analyzed herein, significant downward pump-ing and accumulation of magnetic fields takes place at the base of the convecting layers. Thisis a robust characteristic of numerical simulations of turbulent convection in density-stratifiedenvironments, and can be traced to the topological asymmetry between strong narrow down-flows of cold fluid, and the gentler broader upflows of warm fluid (see, e.g., Tobias et al. 2001)As can be seen on Figure 2, the magnetic field accumulates and peaks in the outer reachesof the underlying stably stratified fluid layers. This behavior is generally consistent withprevalent views of sunspot formation, which posit the storage of toroidal magnetic flux ropesin the overshoot layer, prior to their buoyant destabilization and rise to the photosphere (Fan(2009); but see Stenflo & Kosovichev (2012); Nelson et al. (2013) for alternative viewpoints).The background stratification used in our EULAG-MHD simulations is characterizedby a very strongly stable stratification in the fluid layers underlying the convection zone.More specifically, we use a layered polytropic model, with index varying from a value 1.5at the base of the convecting layers ( r/R = 0 . r/R = 0 . E Kr associated with the radial component of the flow (solidline) and the root mean square latitudinal deviation ∆Ω of the zonally-averaged plasmaangular velocity (dash-dotted line), both integrated over spherical shells and time-averagedover the whole simulation: E K r ( r ) = ∆ r T (cid:90) T (cid:90) π (cid:90) π ( u (cid:48) r ) ( r, θ, φ, t ) ρ ( r ) r sin θ d θ d φ d t , (4)∆Ω ( r ) = (cid:114)(cid:68)(cid:0) Ω ( r, θ ) − ¯Ω ( r ) (cid:1) (cid:69) θ , (5)where Ω( r, θ ) = 12 πr sin θ (cid:90) π u φ ( r, θ, φ, t )d φ , (6)and ¯Ω( r ) denotes spatial averaging over a spherical shell of radius r , and ∆ r is a shell 9 –Fig. 3.— The solid and dashed lines show the depth variation of the kinetic energy associatedrespectively with radial and horizontal (latitudinal and zonal) fluid motions, integrated onconcentric spherical shells of thickness ∆ r/R = 0 . r/R = 0 .
69 (vertical blueline), followed by slower decrease by another two orders of magnitude down to the base of thestable layer. The finite thickness of the layer across which the kinetic energy drops to zerois due to convective overshoot; the upward-directed buoyancy force below the convectivelyunstable layers does not decelerate instantaneously the strong, convective downflows enteringthe stable layer from above, so that convective mixing persists in a thin layer underlying theconvectively unstable fluid layers. We opted to define our “overshoot layer” as the depthinterval over which the first sharp drop is taking place, i.e., 0 . ≤ r/R ≤ . . ≤ r/R ≤ . E Kr as the top of the simulation domain isapproached from below reflects our (conventional) choice of upper boundary condition onthe flow for such simulations, namely stress-free and impenetrable.The latitudinal differential rotation, as measured by ∆Ω (dash-dotted line on Figure3), is roughly constant in the bottom half of the convecting layers, and also drops in thestable layers, albeit more slowly than the flow kinetic energy. This is due primarily tolarge-scale magnetic torques contributing to radial fluxes of angular momentum down to r/R (cid:39) .
65 in this simulation (see Figure 6 in Beaudoin et al. 2013). Despite a drop bya factor of ∼ r/R = 0 .
04 (Charbonneau et al. 1999a). The slow increase of ∆Ω in the outer half ofthe convecting layers can be traced to the strong equatorial differential rotation building upthere (see Figure 2B).The dashed line on Figure 3 shows the depth variation of the kinetic energy E Ka asso-ciated with the non-axisymmetric horizontal (i.e., latitudinal and zonal) component of theflow, plotted on the same logarithmic scale as E Kr , again integrated over spherical shells andtime-averaged over the whole simulation, similarly to our definition of E Kr in eq. (4): E K a ( r ) = ∆ r T (cid:90) T (cid:90) π (cid:90) π ( u (cid:48) θ ) + ( u (cid:48) φ ) ) ρ ( r ) r sin θ d θ d φ d t . (7)The kinetic energy of the non-axisymmetric horizontal components also stays roughlyconstant through the bulk of the convecting layers, as does E Kr . It also undergoes a sharpdrop in the overshoot layer, but only by two orders of magnitude, as opposed to the 15order-or-magnitude drop observed in E Kr . A further, gradual decrease by another threeorders-or-magnitude takes place between the base of the overshoot layer and the bottom ofthe domain. This indicates the presence of predominantly “horizontal” fluid motions, i.e.fluid motions constrained to constant radius spherical shells, sustained in the stable layers.This represents a first hint of dynamical effects developing in the stable layer, distinct fromthe simple buoyant deceleration of convective downflows. We will revisit this issue in §
4, butas a needed preamble we first examine in greater detail the physical mechanism(s) leadingto magnetic field accumulation and amplification in the stable layers. 11 –
3. Origin of the magnetic fields in the stable layer
Figure 4 offers a more detailed look at the localisation of the large-scale magnetic fieldbands in relation to the overshoot and tachocline layers defined on the basis of Figure 3.This is the same zonally-averaged toroidal magnetic field at cycle peak displayed on Figure2A, but plotted this time in a cartesian radius-latitude plane. The toroidal field bands peakat (cid:39) . ≥ . ∼ S = 1 µ E × B = 1 µ ( u × B ) × B , (8)and integrate its radial component on spherical shells: P ( r, t ) = (cid:90) π (cid:90) π S r ( r, θ, φ, t ) r sin θ d θ d φ . (9)Figure 5 shows the depth variation of this quantity, time-averaged over the whole simu-lation (solid line), as well as equivalent profiles extracted at an epoch of magnetic cycle peak(dashed line) and polarity reversal (dotted line). In all cases the integrated radial Poyntingflux is negative at all depths except in subsurface layers, as a consequence of our imposedimpenetrable upper boundary condition. Magnetic field accumulation is expected whereverthe divergence of the Poynting flux is negative, which here is generally the case in the depthrange r/R ≤ . r/R = 0 . r/R = 0 . r/R = 0 . ◦ latitude. The minima occur shortly after the polarity reversal,consistent with the subsequent buildup of the deep toroidal field at and below the base of theconvecting layers. The downwards slant across the overshoot layer is the direct counterpartof the dashed oblique lines on Figure 6. Note, however, the break of slope at the baseof the overshoot layer, suggestive of an additional —and possibly local— inductive and/ortransport process contributing to the spatiotemporal variations of the Poynting flux.We now compute the total magnetic energy content ( E M ) in the overshoot and stablelayers by direct integration of the simulation output in the depth range 0 . ≤ r/R ≤ . ◦ latitude in the N-hemisphere, over a ∼
200 yr segment of the millenium simulation. The soliddots indicate epochs of minimal shell-integrated Poynting flux at a few depths going fromthe base of the convection zone, through the overshoot layer and into the stable zone. Thegreen and blue dots identify the minima in the correspondingly colored time series on Figure6. 15 –at an epoch of cycle maximum: E M = 12 µ (cid:90) π (cid:90) π (cid:90) . . B r sin θ d r d θ d φ . (10)The timescale for field accumulation beneath the convection zone is then directly obtainedby dividing this quantity by the value of the shell-integrated Poynting flux at r/R = 0 . τ = E M P . (11)The resulting numerical value is τ (cid:39) ∼
40 yr magnetic half-cycleperiod in our simulation. This suggests that sufficient magnetic energy is being providedby downward turbulent pumping to destroy the deep toroidal flux bands of the precedinghalf-cycle and rebuild a band of opposite polarity in each hemisphere. This result is thusconsistent with a minimalistic scenario whereby the buildup and reversal of the toroidal fluxbands below the convection zone are merely a passive side-effect of global dynamo action within the convection zone (on this point, see also Masada et al. 2013).However, τ is but a rough estimate, and other mechanism are potentially at play. Asnoted already in the context of Figure 3, some level of latitudinal differential rotation, bothradial and latitudinal, persists across the overshoot layer and well within the stable zone.With a significant large-scale poloidal magnetic component also present, shearing by dif-ferential rotation can contribute to the induction of a toroidal component. This process iscaptured by the zonal component of the MHD induction equation, which in the ideal limitand for axisymmetric large-scale magnetic field and differential rotation reduces to ∂B φ ∂t = r sin θB r ∂ Ω ∂r + sin θB θ ∂ Ω ∂θ (12)Separately integrating B φ and each term on the RHS of the above expression over the radialand latitudinal extent of the toroidal field at a time of solar maximum allows to computecharacteristic timescales for induction by the radial and latitudinal shear, again by simplydividing the first by the other two. Both of these timescales end up at (cid:39)
30 yr; whilesignificantly larger than the timescale (11) associated with the downward Poynting flux,these are still comparable to the half-cycle period for (cid:39)
40 yr. One can but conclude thatdifferential rotation shear contributes significantly to toroidal field induction in the overshootand stable layers.On the basis of this new estimate, an additional factor is thus added to the minimalscenario outlined above: the buildup and reversal of the deep-seated magnetic field bandsoccurs primarily through turbulent pumping from above, with additional amplification pro-vided by the differential rotation shear. We could stop here, but one key piece of evidencecompels us to push our analysis further, namely the existence of significant horizontal kineticenergy in the stable layer (dashed line on Figure 3). 16 –
4. Searching for the signature of instabilities
As noted earlier, the kinetic energy of the non-axisymmetric horizontal flow components( E Ka shown in Figure 3) remains at significantly high values well into the stable layer, beforefinally dropping abruptly by several orders of magnitudes at the lower boundary, primarilya consequence of the artificial friction terms introduced there. These horizontal flows persistwell below the depth of convective overshoot (blue vertical line on Figure 3), and so cannotbe directly driven by the overlying convection. One possibility is that they arise through thedevelopment of one (or more) local fluid instability.The solar tachocline is the site of significant latitudinal and radial rotational shear, andits upper reaches include the overshoot layer, which is only weakly subadiabatic. Thesecharacteristics are known to be conducive to the development of a wide array of fluid insta-bilities (see, e.g., Tassoul 2000, § ∼ . Numerous analytic, semi-analytic and numerical calculations have identified two-dimensionalMHD instabilities that can become excited in stably stratified regions of the solar interiorin the presence of latitudinal differential rotation and large-scale magnetic field. Particu-larly pertinent to our simulation is the so-called joint MHD instability first investigated byGilman & Fox (1997) (see also Dikpati & Gilman 1999, 2001a; Dikpati et al. 2003; Gilmanet al. 2007). This instability develops in stably stratified environments in the presence of ax-isymmetric latitudinal differential rotation and large-scale toroidal magnetic fields, i.e., the 17 –situation prevailing in the outer reaches of the stable layer in our simulation. The instabilityplanforms are 2D, i.e., they develop on spherical surfaces, and the most unstable modeshave low azimuthal wavenumbers, m = 1 or 2, as magnetic tension provides a restoringforce that strongly suppresses higher wavenumber modes (see Gilman et al. 2007, and refer-ences therein). Depending on the cross-hemispheric phasing of the instability planform, theglobal development of the instability can lead to magnetic reconnection across the equator(“clamshell instability”; see, e.g., Cally et al. 2003), or the two toroidal flux systems can bothtilt while remaining parallel to one another across the hemispheres (“tipping instability”).The magnetoshear instabilities tap into both the kinetic energy of the differential rota-tion, as well as the magnetic energy of the large-scale toroidal magnetic field. It is a “joint”instability, in that both ingredients are required for the instability to grow. The growthrate s of the most unstable mode is typically some fraction of the Alfv´en time based on thetoroidal magnetic field strength: s ∼ Lu A , u A = (cid:115) (cid:104) B φ (cid:105) µ ρ , (13)where L is a typical length scale for the toroidal magnetic field. Using L = R , (cid:104) B φ (cid:105) = 0 . T ,and ρ = 10 − kg m − , appropriate for the upper part of the stable layer in our simulation,yields s (cid:39) t ) = 12 µ (cid:90) . . (cid:90) [ π/ ,π ][0 ,π/ (cid:90) π ( B (cid:48) ) r sin θ d φ d θ d r , (14)TFME( t ) = 12 µ (cid:90) . . (cid:90) [ π/ ,π ][0 ,π/ (cid:90) π (cid:104) B φ (cid:105) r sin θ d φ d θ d r . (15) 18 –Unlike in Miesch (2007)’s model, here the poloidal field is generated autonomously; itwill therefore prove useful to define an analog of (15) for the axisymmetric poloidal field:PFME( t ) = 12 µ (cid:90) . . (cid:90) [ π/ ,π ][0 ,π/ (cid:90) π ( (cid:104) B r (cid:105) + (cid:104) B θ (cid:105) ) r sin θ d φ d θ d r . (16)This proxy turns out to be largely dominated by the contribution from the latitudinalcomponent (cid:104) B θ (cid:105) . Note that in order to ensure cleaner proxy time series, we only integrateover the stable the layer, excluding the overshoot layers. Moreover, as detailed in Passos &Charbonneau (2014), magnetic cycles in the EULAG-MHD millenium simulation can showsmall but significant phase lag between hemispheres; all of the above proxy integrals aretherefore computed separately for the Northern (0 ≤ θ ≤ π/
2) and Southern ( π/ ≤ θ ≤ π )hemispheres.Figure 8 shows a representative 400yr segment of the hemispheric time series for theTFME (solid line) and NAME (dashed line) proxies. Both proxies wax and wane cyclicallywhile maintaining a well-defined phase lag, with NAME peaking in the late descending phaseof TFME, and the onset of the growth phase of NAME almost always occurring near the peakof TFME. This is remarkably similar to the pattern characterizing the forced 2D simulationsof Miesch (2007) (compare his Figure 2 to Figure 8 herein). This suggests —of coursewithout strictly proving— that we are observing in the millenium simulation the same typeof magnetoshear instability investigated by Miesch (2007); however, no external forcing isimposed here, as the latitudinal differential is being maintained by Reynolds stresses withinthe convection zone, and the poloidal field is naturally replenished by turbulent dynamoaction within the convecting layers and subsequent downward pumping of the magneticfield. In particular, the roughly similar amplitudes of variation in TFME and NAME, ∼ × J , are consistent with the cyclic exchange of energy between the axisymmetric andnon-axisymmetric magnetic components characterizing the forced shallow water simulationsof Miesch (2007).Under flux freezing, any global displacement of the toroidal field bands due to themagnetoshear instability —whether developing in its clamshell or tipping variants— must beaccompanied by a corresponding latitudinally-oriented flow constrained to spherical shells.Figure 3 already indicates the predominance of horizontal flows in the stable layer, andthe characterization of such flow can provide further evidence that the instability is indeedoperating. Accordingly, we now define two measures of flow kinetic energy similar to thoseused for magnetic energy: the total kinetic energy of the poloidal (PKE) and toroidal (TKE)components. As before, we integrate separately the Northern and Southern hemisphere:PKE( t ) = 12 (cid:90) . . (cid:90) [ π/ ,π ][0 ,π/ (cid:90) π (cid:0) u θ + u r (cid:1) ρr sin θ d φ d θ d r . (17) 19 –Fig. 8.— Time series of the volume-integrated magnetic energy contained in the axisymmet-ric toroidal field (TFME; solid line) and the total non-axisymmetric field (NAME; dottedline) in the stable zone of the EULAG-MHD millenium simulation, for the Northern (top)and Southern (bottom) hemispheres. Compare to Figure 2 in Miesch (2007).TKE( t ) = 12 (cid:90) . . (cid:90) [ π/ ,π ][0 ,π/ (cid:90) π u φ ρr sin θ d φ d θ d r . (18)Note that we use here the total flow components, but since very little axisymmetricmeridional flow develops in the stable layer, PKE is essentially the flow equivalent of NAME.Moreover, as with the PFME proxy, PKE is dominated by the latitudinal component u θ .Figure 9 shows these two time series over a restricted 200 yr temporal span, together withthe magnetic proxies defined earlier, all in the Northern hemisphere. The dashed verticallines indicate the epoch of peak NAME in the descending phase of each of the five magneticcycles developing over the time period covered. Note that our (axisymmetric) poloidal fieldmagnetic energy PFME also exhibits cyclic behavior, in response to dynamo action in theoverlying convecting fluid layers, unlike in the simulation of Miesch (2007) where the poloidal(latitudinal) field is imposed externally and remains constant.Examination of Figure 9 reveals that both PKE and TKE evolve in phase with NAME,with peak-to-through variations ∼ . . × J, over an order of magnitude lowerthan the corresponding variations in TFME. Indeed, the combined rise of NAME, TKE, andPKE adds up to ∼ . × J , of the same order but still comfortably smaller than theassociated ∼ × J drop in TFME. The in-phase variation of PKE with NAME suggests 20 –Fig. 9.— The top panel shows a 200yr closeup of Figure 8 for the N-hemisphere, and thebottom shows the corresponding time series for the PKE (solid line) and TKE (dotted line)flow energy proxies defined through eqs. (17)–(18). The poloidal magnetic proxy PFME(viz. eq. 16) is also shown. The vertical dashed lines indicate the epochs of peak NAME.Note the in-phase variations of PKE and TKE with NAME. 21 –that the horizontal flow develops simultaneously with the growth of the non-axisymmetricmagnetic component, which is the pattern expected for the clamshell or tipping variants ofthe magnetoshear instability. As plasma is displaced away from or towards the rotation axis,conservation of angular momentum causes zonal deceleration or acceleration, leading to agrowth of TKE in phase with PKE, as is indeed observed on Figure 9. In order to properly identify the clamshell –or tipping variants– of the magnetoshearinstability, we decompose the magnetic field on the spherical harmonics in the stable layer.We use the classical vectorial spherical harmonics basis (see Rieutord 1987; Mathis & Zahn2005; Strugarek et al. 2013) B = (cid:88) l,m α l,m R ml + β l,m S ml + γ ml T ml , (19)where the orthogonal basis ( R ml , S ml , T ml ) is defined by R ml = Y ml e r S ml = ∇ ⊥ Y ml = ∂ θ Y ml e θ + θ ∂ ϕ Y ml e ϕ T ml = ∇ ⊥ × R ml = θ ∂ ϕ Y ml e θ − ∂ θ Y ml e ϕ , with Y ml being the classical spherical harmonics. At a given depth, the magnetic energyspectrum can be decomposed along the spherical harmonic degrees m byME m ( r, t ) = V r µ (cid:88) l ≥ m | α l,m | + l ( l + 1) (cid:0) | β l,m | + | γ l,m | (cid:1) , (20)where V r is the local volume of the spherical shell centered on the cell located at depth r ,included for dimensional consistency.We display the evolution of the magnetic energy for m ∈ [0 ,
5] in Figure 10 at depth r/R = 0 . m = 1 traces the onset of a magnetoshear instabilitywhich is synchronized with TFME (see Figure 9). The m = 1 growth is accompanied byan m = 2 component of the magnetic energy, which saturates later in the cycle, in phasewith the higher m components and with NAME. The non-axisymmetric magnetic energy(NAME) cycle maxima are dominated by the m = 3 or m = 4 modes, depending on theconsidered cycle. Once the tipping instability has taken off, magnetic energy stored in themost unstable m = 1 and m = 2 modes is expected to be non-linearly transfered to thosehigher m ’s due to the simultaneous development of horizontal flows at the same scales. 22 – Time [Year] M E m [ J ] Magnetic energy, r / R = TFME+PFME ( m = ) NAME ( m = ) m = m = m = m = m = Fig. 10.— Magnetic energy evolution at r/R = 0 . m , summed over l . The total, axisymmetric ( m = 0) magnetic energy,which is dominated by the PFME contribution, is indicated in black. The non-axisymmetric(NAME, m (cid:54) = 0) in shown in bold grey. The three vertical lines indicate the time at whichthe mollweide projection in Figure 11 were taken. 23 –In order to geometrically identify the tipping instability that develops on the shellsharboring the strong toroidal magnetic fields, we again follow Miesch (2007) and visualizemagnetic field lines as iso-contours of the magnetic streamfunctions J and ˜ J , defined suchthat B ∼ ˆ z × ∇ J and B − (cid:104) B (cid:105) ϕ ∼ ˆ z × ∇ ˜ J . These fields lines are shown on Figure 11,constructed at a depth r/R = 0 .
680 well below the overshoot layer at epochs of m = 1maximum (left), TFME maximum (middle) and m = 3 maximum (right). J Max m = ˜ J Max m = Max m = Fig. 11.— Magnetic field lines on the spherical surface r/R = 0 . m = 1 maximum, m = 0 maximumand m = 3 maximum. The upper diagrams correspond the streamfunction of the full hori-zontal magnetic field, and the lower diagrams to the streamfunction of the non-axisymmetrichorizontal magnetic field. On these Mollweide projections a purely toroidal magnetic fieldwould have all its field lines oriented horizontally. Compare to Figure 3 in Miesch (2007).The epoch of maxima of m = 1 (left panels) generally occurs concurrently with thepeak of TFME, before the maxima of PFME (middle panels). The field lines are thuspredominantly axisymmetric and composed of longitudinal field at mid and high latitudes,and latitudinal field near the equator. A hint of global m = 1 tilting can be seen in the twoupper left panels, which are confirmed by the predominance of an m = 1 structure in thetwo lower left panels. At the time m = 3 is maximized it dominates the magnetic energyspectrum and no clear longitudinally aligned field lines can be observed on the right panels.Albeit m = 3 dominates, the field lines exhibit a complex pattern of mixed m componentsthat populate the magnetic energy spectrum at the peak of NAME.Further insight is gained by quantifying the spectral energy transfers leading to thesuccessive growth of the m = 1 and m = 3 modes. During its maxima periods, the m = 1 24 –mode completely dominates the non-axisymmetric spectrum. We choose to focus on the m = 3 mode because it provides a typical dominant mode during the NAME maxima periods,and is always one of the few most energetic modes in each cycle (the m = 4 and m = 5 modessometimes dominate the NAME spectrum, but they significantly vary from one cycle to theother). We follow the procedure detailed in Strugarek et al. (2013) to quantify the amount ofenergy transfered to a given m due to the triadic interactions involving two other sphericalharmonics ( l , m ) and ( l , m ). We display in Figure 12 energy transfers maps (summedover l ) towards m = 1 (left column) and m = 3 (right column) during their growth phase fora typical cycle. Hence, the left column was averaged during a different time-window thanthe right column, to account for the time-lag between the growth periods of the two modes(see also Figure 10). Positive (red) and negative (blue) transfers respectively representa source and a sink of magnetic energy for the chosen m . In each panel, the horizontalaxis corresponds to the velocity field spectrum, and the vertical axis to the magnetic fieldspectrum. The contributions are separated in non-axisymmetric ( m (cid:54) = 0) and axisymmetric( m = 0) fields, leading to three possible ways to transfer energy.We show in the upper panels the transfers involving the differential rotation, and in themiddle panels the transfers involving the large scale magnetic field. The axes correspondto the various l couplings, which allows us to identify the dominant transfers from thedifferential rotation ( l = 3, 5, corresponding to the vertical red stripes in the upper panels)and the large-scale magnetic field ( l = 3, corresponding to the horizontal red stripes in themiddle panels). We see that a large range of non-axisymmetric scales l , coupled to theaxisymmetric fields, are involved in the transfer. In the lower panels, we focus on the fullynon-axisymmetric couplings and the axes now represent the m modes of the magnetic andvelocity fields, summed over l . The triadic selection rule naturally leads to diagonal-onlytransfers in the lower panels. The transfers to higher m mode play little to no role in the caseof m = 1, whereas the m = 3 mode transfers a significant amount of energy to the m = 4 to10 modes during its growth, as one would expect from a cascade-like process. The color-scaleis saturated at 25% of the total change of the magnetic energy ˙ E M during the consideredgrowth period in the upper and middle panels, and to 10% in the lower panels. It appearsclearly that the m = 1 and m = 3 modes both gain energy from the differential rotationand the large-scale magnetic field, while the non-linear interactions with higher order modestend to stabilize them. The total magnetic energy transfer from each channel is quantifiedby the percentage in the upper left corner of each panel. The m = 1 mode hence dominantlygrows by receiving energy directly from the large-scale differential rotation, while the m = 3mode preferentially draws energy from the axisymmetric large-scale magnetic field at l = 3.Note that for each column the total is not exactly 100% due to the fact that the transfermaps were averaged during the magnetic energy growth phase. These results are robust for 25 –the different growth phases shown in Figure 10.Using the fully detailed analysis from Strugarek et al. (2013) (not shown here), it isfurther possible to disentangle the origin of the energy transfers with respect to, e.g. , theradial or latitudinal structure of the differential rotation. Our simulation shows that the m = 1 mode dominantly draws energy from the latitudinal differential rotation, while the m = 3 mode grows by extracting energy from the radial structure of the axisymmetric large-scale magnetic field. As a result, these energy transfers are a compelling evidence of theoccurence of an MHD instability akin to the magnetoshear instability in our simulation.The magnetic energy nevertheless largely dominates the energy balance in the stablelayers (see Figure 9). The m = 1 mode grows by receiving energy from both the differentialrotation and the large-scale axisymmetric magnetic field, as a result, it is possible that someother type of MHD instability is instead at play, and we now turn to this possibility.
5. Digging further: the Tayler instability
The magnetoshear instability investigated by Miesch (2007) is far from the only onethat can potentially develop in stably stratified, weakly magnetized differentially rotatingastrophysical environments. For example, the turbulent stresses generated by the magneto-rotational instability (MRI; see Balbus & Hawley 1991) are now believed to dominate theoutward transport of angular momentum in weakly magnetized accretion disks. This is avery powerful, local MHD instability, requiring a significant poloidal magnetic component tooperate, but in the stellar interior context it also requires a significant outwardly decreasingradial differential rotation. The high-latitude regions of the tachocline satisfy this criterion,and could thus be the seat of a localized version of this instability (see, e.g., Parfrey & Menou2007; Masada 2011). However, if present throughout the overshoot layer and tachocline, itwould lead to significant radial mixing. This appears ruled out here, on the basis of Figure3 which indicates that fluid motions in the stable layer are strongly restricted to sphericalsurfaces. The magnetic buoyancy instability (Parker 1955b) is also a potential candidate, butagain it would lead to radial mixing, which we do not observe in the stably stratified layersof our simulation. Moreover, we likely do not have the spatial resolution required to capturethe formation of thin magnetic flux tube-like structures conducive to the development ofthis instability, at least judging from the much higher resolutions simulations of Nelson et al.(2013). To the best of our knowledge, at this writing these remain the only global MHDsimulations of solar convection in which the spontaneous onset of this instability has beenobserved. 26 –
U ( l , m = 0 ) B ( l ) m = 1 U ( l , m = 0 )11.1 % m = 3 U ( l ) B ( l , m = ) l )171.6 % U ( | m | ) B ( | m | ) -7.9 % U ( | m | )-82.6 % − . . . E ne r g y t r an s f e r / ˙ E M − . . . E ne r g y t r an s f e r / ˙ E M Energy transfers, r/R = 0.680
Fig. 12.— Energy transfer maps towards the m = 1 (left column) and m = 3 (rightcolumn) components of the magnetic energy spectrum. The transfers are averaged during atypical growth phase of those components, which are summed over l . The colormap denotespositive (red) and negative (blue) energy transfers and is normalized to the variation of themagnetic energy ˙ E M . The transfers are separated in three panels in each column, allowingto identify transfers from the axisymmetric ( m = 0) and non-axisymmetric ( m (cid:54) = 0) velocity(abscissas) and magnetic (ordinates) fields. The upper panels hence show transfers with theaxisymmetric differential rotation, and the middle panels transfers with the axisymmetricmagnetic field, as a function of l . In the lower panels we focus on the fully non-axisymmetrictransfers for which we chose the energy conversions as a function of m , summed over l .The percentage listed at the upper left in each panel corresponds to the sum over all thecontributions in the panel. 27 –Another class of MHD instabilities, the “Tayler instabilities”, have also been investigatedin detail in the context of various types of large-scale magnetic fields embedded in stellarradiative interiors. Particularly relevant in the present context is the instability of a purelytoroidal magnetic field in the ideal MHD limit Tayler (1973); Spruit (1999). In the absenceof rotation, any such magnetic field B φ ( r, θ ) is unstable, no matter how weak the field is.The instability is most prone to develop close to the magnetic symmetry axis, where B φ = 0,as no restoring force can resist the magnetic force pointing towards the axis (see, e.g., thediscussion in Spruit (1999); also Spruit (2002); Brun & Zahn (2006)). The instability candevelop both axisymmetric and non-axisymmetric planforms, with the lowest zonal mode(azimuthal wavenumber m = 1 in a spherical harmonic expansion) usually most unstable inthe latter case (Zahn et al. 2007), because these develop the least magnetic tension tendingto oppose the instability. As with the magnetoshear instability investigated in the precedingsection, the growth rate for the Tayler instability is of the order of the Alfv´en time. Stabilitycriteria for axisymmetric purely toroidal magnetic fields have been obtained by Goossenset al. (1981) for spherical geometry, and take the form:14 πr sin θ (cid:20) H φ cos θ − sin θ cos θ ∂H φ ∂θ (cid:21) > , m = 0 , (21)14 πr sin θ (cid:20) H φ (cid:0) m − θ (cid:1) − sin θ cos θ ∂H φ ∂θ (cid:21) > , m = 1 , (22)where H φ = b l ( r ) P l (cos θ ). For the purposes of the foregoing analysis we shall simplyassume H φ ≡ (cid:104) B φ (cid:105) ( r, θ ), which leads to:2 cos θ − θ sin θ ∂ ln (cid:104) B φ (cid:105) ∂θ > , m = 0 , (23)1 − θ − θ sin θ ∂ ln (cid:104) B φ (cid:105) ∂θ > , m = 1 . (24)These stability criteria are admittedly obtained in idealized conditions differing signif-icantly from those encountered within our numerical simulation: ideal MHD, no rotation,purely toroidal magnetic field. However, Pitts & Tayler (1985) showed that rotation weak-ens but does not suppress the instability. Likewise, the instability is also weakened, but notsuppressed, in the presence of a large-scale poloidal magnetic field component. In fact, thenumerical simulations of Braithwaite (2006) confirm the generally unstable nature of large-scale magnetic fields in stably stratified, radiative stellar interiors, and also indicate that themost stable large-scale magnetic field configurations have poloidal and toroidal componentsof comparable strengths (Braithwaite & Nordlund 2006). Interestingly, and as noted ear-lier in § m = 0) form of the Tayler instability, in a manner consistent witheq. (21), in her 2D axisymmetric MHD simulations of the solar radiative interior includinga poloidal magnetic field and imposed differential rotation (a forcing setup somewhat as inMiesch 2007). Her simulations show that even in the absence of significant bona fide dy-namo action, as long as the poloidal component persists the differential rotation can inducea toroidal component which, upon becoming unstable to the axisymmetric Tayler instabil-ity, undergoes polarity reversals (see her Figure 6). Even closer to the physical situation ofinterest here, Brun & Zahn (2006) performed 3D MHD simulation of the solar tachoclinein which they also observe the development of what they suggest is the non-axisymmetricform of the Tayler-like instabilities, persisting at all depths and particularly prominent inthe vicinity of the polar axis, as expected of this instability.Figure 13 shows a time-latitude representation of the m = 1 stability criterion (eq. 24),constructed at the upper extent of the stable layer in the simulation ( r/R = 0 . m = 1 Tayler stabilitycriterion, i.e., the LHS of eq. (24). The superimposed thick black lines are the NAME timeseries for the Northern and Southern hemispheres, the latter assigned negative values forthe purpose of clarity and symmetry. The thin dashed line is the (cid:104) B φ (cid:105) = 0 isocontour, andare useful in identifying the spatiotemporal unfolding of magnetic polarity inversion at thisdepth.The large-scale axisymmetric toroidal field does turn out to be stable at most latitudesand phases of the cycle, the stability criterion being violated mostly in the vicinity of thepolarity inversion line, and at high latitudes. As argued earlier, marginal stability is in factexpected across most of the domain in our nonlinearly-saturated, statistically stationarysimulation. One may however also expect the linear stability criteria to be violated whereverand whenever the instability is turning on, and has not yet had time to reach saturation andalter the background magnetic field profile. Taken at face value, Figure 13 indicates that 29 –Fig. 13.— Time-latitude color rendering of the m = 1 instability criterion as given byeq. (24), based on the magnetic field extracted at depth r/R = 0 . ± ◦ latitude, corresponding to the poleward extent of the strong ( (cid:104) B φ (cid:105) ≥ . (cid:104) B φ (cid:105) = 0 isocontour. 30 –the Tayler instability is first triggered at polar latitude, which coincides temporally withthe onset of the growth phase in the NAME proxy (thick black lines), and subsequentlymoves to progressively lower latitudes. Once the “instability front” reaches around ± ◦ latitude, corresponding to the poleward edge of the mid-latitudes magnetic field bands, thefront becomes more slanted, suggesting that the higher magnetic energy available in the fieldbands accelerates the development of the instability, as evidenced by the crossing horizontaland vertical thin black lines on Figure 13. This is also when the NAME proxy reachesits peak before beginning to drop, signaling that the (non-axisymmetric) instability is nowturning itself off via the destruction of the axisymmetric toroidal field bands making up itsenergy reservoir. After the instability shuts off, there follows a “quiet” interval during whichthe axisymmetric toroidal magnetic bands start to rebuild once again, up to the point wherethe instability will once again be triggered at high latitudes.Interestingly, a plot similar to Figure 13 but constructed for the m = 0 stability cri-terion (23) reveals a spatiotemporal pattern closely resembling the m = 1 case, but withreduced amplitude, consistent with the suppression, by rotation, of the axisymmetric Taylerinstability in favor of its non-axisymmetric m = 1 cousin, as suggested by the analysis ofPitts & Tayler (1985).
6. Additional numerical experiments and a plausible scenario
In this paper we have investigated the possible occurrence of non-axisymmetric MHDinstabilities developing in the subadiabatic, stably stratified fluid layers underlying the con-vectively unstable layer in the EULAG-MHD “millenium simulation” described in (Passos &Charbonneau 2014). The cyclic, phase-lagged waxing and waning of the magnetic energiesassociated with the axisymmetric toroidal and non-axisymmetric total magnetic field com-ponents, as extracted from the stable layer in the simulation, bears a striking resemblance tothat characterizing the magnetoshear instability studied by Miesch (2007) using a shallow-water MHD model with forced differential rotation and poloidal magnetic component. Inour simulation this forcing occurs naturally through Reynolds stresses and turbulent dynamoaction taking place in the overlying convecting layers. Our analysis suggests that somethingakin to this magnetoshear instability is operating in our simulation. The analysis of energytransfers unambiguously shows that the magnetic m = 1 mode grows by receiving energyfrom the latitudinal differential rotation, and that the dominant magnetic m = 3 modemainly feeds from the radial structure of the large-scale axisymmetric magnetic field.Motivated by the rather weak latitudinal differential rotation characterizing the stablelayers in our simulation, we also investigated the possibility that the Tayler instability be 31 –operating along with the magnetoshear instability. Further motivation is also found in somerecent numerical studies (Brun & Zahn 2006; Rogers 2011; Parfrey & Menou 2007) whichhave uncovered various elements of evidence for the action of this instability in a similarsolar context. The Tayler instability does not require differential rotation, as it taps onlyinto the magnetic energy of the underlying large-scale magnetic field to power itself, and isparticularly prone to develop in the vicinity of the magnetic symmetry axis. Using linearstability as a guide, we have found some suggestive evidence for the action of the Taylerinstability as an agent contributing to the destruction of the large-scale toroidal magneticfield bands building up in the simulation, especially on their poleward edge.Whatever the exact nature of the instability that may be at play in the stable layers ofour simulation, the crucial question is: does it play a significant —or maybe even important—role in the global dynamo process leading to amplification and polarity reversals of thelarge-scale magnetic field? The fact that the Poynting flux remains downward directedat all phases of the magnetic cycle and at all latitudes in the lower convection zone andovershoot layer certainly suggests that there is no direct upward electromagnetic feedbackfrom the stable layer into the convecting fluid layers in the simulation. The obvious empiricaltest, running a simulation without the stably stratified fluid layer but otherwise identical,is unfortunately inconclusive. Such a simulation does not produce a decadal large-scalecycle, but the differential rotation and cyclonic character of turbulence in the bottom of theconvection zone also turn out markedly different, due to the boundary conditions that mustbe imposed there. It is therefore impossible to tell whether the lack of large-scale magneticcycle is due to a different mode of dynamo action within the convection zone, rather than tothe absence of the stable layer. Dikpati & Gilman (2001a) also proposed that the tachoclineinstabilities could generate a net kinetic helicity leading to a local α -effect contributing insitu to the build-up of the local poloidal field. In our simulation the kinetic helicity patternassociated with the instability is compatible with the Dikpati & Gilman (2001a) scenario,but the estimated amplitude of the associated α -effect, computed following the methodologyintroduced by these authors, remains far too small to have any influence on the magneticcycle timescales.Another option for numerical experimentation is to retain the stable layer, but alter itsdegree of subadiabatic stratification. Figure 14 shows the result of such an experiment. Thissimulation is in all aspects identical to the millenium simulation used herein, except that nowthe polytropic index varies linearly with depth from a value 1.5 at the base of the convectinglayers up to 3.0 at the base of the domain, whereas in the millenium simulation runs thecorresponding linear variation is from 1.5 to 2.5. As a result, this simulation is more stablystratified everywhere in the radiative zone compared to the millenium simulation. The toppanel of Figure 14 shows a time-latitude diagram of the zonally-averaged toroidal magnetic 32 –field at r/R = 0 . (cid:39)
500 yr the Northern hemisphere shuts down,followed by the Southern at (cid:39)
550 yr.Fig. 14.— Time-latitude diagram of the zonally-averaged toroidal magnetic field at r/R =0 .
711 in a EULAG-MHD simulation identical to the millenium simulation, except for a morestrongly subadiabatic polytropic profile in the stable layer. The bottom panel shows a time-longitude diagram extracted at r/R = 0 .
711 and 45 ◦ latitude North, in the same simulation.The apparent disappearance of the cycle on the top panel is a consequence of the dynamoswitching to a non-axisymmetric large-scale mode, still undergoing polarity reversals withrespect to its own symmetry axis (see text).Has the simulation fallen here into a Maunder Minimum-like state of strongly suppressedcyclic activity? The answer is no, as evidenced on the bottom panel of Figure 14. Thisshows a time- longitude diagram of the toroidal magnetic component, not averaged zonally,and extracted at 45 ◦ latitude and the same depth as the top panel. Far from vanishing,the large-scale magnetic cycle persists with similar amplitude and a shorter half-period of (cid:39)
30 yr, but now its symmetry axis has undergone a large tilt with respect to the rotationaxis.At the very least, these numerical experiments indicate that the structure of the stable 33 –layers affects the long term stability of the various types of large-scale cyclic dynamo modesthat can materialize in these EULAG-MHD simulations. But how about the polarity reversaltaking place in the more solar-like axisymmetric mode of large-scale dynamo action arisingin the millenium simulation? The following “global dynamo scenario” is consistent withall analyses presented in this paper, with the global 3D MHD simulations of Masada et al.(2013), as well as with the more geometrical restrictive simulations of Miesch (2007) andRogers (2011):1. Dynamo action is driven primarily within the convection zone, through the differentialrotation and turbulent electromotive force materializing therein. This is consistentwith the analyses of Racine et al. (2011) and Simard et al. (2013), which indicate amode of large-scale dynamo action resembling the so-called α Ω dynamos of mean-fieldtheory;2. Downward pumping of the magnetic field produced in the convecting layers leads tothe buildup of strong zonally-aligned magnetic field bands in the overshoot layer andupper reaches of the underlying stably stratified fluid layer, where further amplifica-tion of the toroidal magnetic component takes place through shearing of the poloidalmagnetic field by differential rotation; downward pumping is observed in virtually allextant MHD simulations of solar convection, and differential rotation shear is a key“ingredient” of most extant solar cycle models;3. Once the toroidal magnetic field bands in the stable layers reach sufficient strength,likely of the order of a few tenths of Tesla, MHD instabilities set in, destabilizing themagnetic field bands and accelerating their dissipation;4. Meanwhile the global magnetic polarity has reversed in the convecting layers, anddownwards pumping of magnetic field of opposite polarity to that having formerlybuilt up in the stable layer begins, eventually leading to polarity reversals therein aswell, closing the dynamo loop towards step 1 above.Additional simulations and analyses are underway to further validate this scenario,which at this juncture remains speculative but plausible. From the point of view of kinematicmean-field and mean-field-like dynamo models of the solar cycle, the instabilities for whichwe have presented evidence in this paper can be considered as contributing to the enhancedturbulent magnetic diffusivity that is essential for such dynamo models to produce solar-likecycles with decadal periods.
Acknowledgements
34 –We wish to thank Sacha Brun useful discussions on the Tayler instability and insightfulcomments upon reading a first draft of this paper. We also thank an anonymous refereefor very useful comments and suggestions that helped strengthen the analysis developed inthis work. The EULAG-MHD “millenium simulation” was originally designed and executedby Mihai Ghizaru. N. Lawson wishes to thanks the physics department of Universit´e deMontr´eal for the award of a graduate fellowship. A. Strugarek is a National PostdoctoralFellow of the Canadian Institute for Theoretical Astrophysics. P. Charbonneau is supportedprimarily by a Discovery Grant from the Natural Sciences and Engineering Research Councilof Canada. All simulations reported upon in this paper were performed on the computinginfrastructures of Calcul Qu´ebec, a member of the Compute Canada consortium.
REFERENCES
Alvan, L., Strugarek, A., Brun, A. S., Mathis, S., & Garcia, R. A. 2015, A&A, 581, A112Augustson, K., Brun, A. S., Miesch, M., & Toomre, J. 2015, ApJ, 809, 149Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Beaudoin, P., Charbonneau, P., Racine, E., & Smolarkiewicz, P. K. 2013, Sol. Phys., 282,335Beaudoin, P., Simard, C., Cossette, J.-F., & Charbonneau, P. 2015, Submitted to ApJBraithwaite, J. 2006, A&A, 449, 451Braithwaite, J., & Nordlund, ˚A. 2006, A&A, 450, 1077Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731,69Brun, A. S., & Zahn, J.-P. 2006, A&A, 457, 665Cally, P. S., Dikpati, M., & Gilman, P. A. 2003, ApJ, 582, 1190Charbonneau, P. 2010, Living Rev. Sol. Phys., 7, 3Charbonneau, P., Christensen-Dalsgaard, J., Henning, R., et al. 1999a, ApJ, 527, 445Charbonneau, P., Dikpati, M., & Gilman, P. A. 1999b, ApJ, 526, 523Cossette, J.-F., Charbonneau, P., & Smolarkiewicz, P. 2013, ApJ, 777, L29 35 –Dikpati, M., & Gilman, P. A. 1999, ApJ, 512, 417—. 2001a, ApJ, 551, 536—. 2001b, ApJ, 559, 428Dikpati, M., Gilman, P. A., & Rempel, M. 2003, ApJ, 596, 680Domaradzki, J. A., Xiao, Z., & Smolarkiewicz, P. 2003, PoF, 15, 3890Fan, Y. 2009, Numerical Modeling of Space Plasma Flows, 416, 489Fan, Y., & Fang, F. 2014, ApJ, 789, 35Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. 2010, ApJ, 715, L133Gilman, P., & Dikpati, M. 2014, ApJ, 787, 60Gilman, P. A. 2015, ApJ, 801, 22Gilman, P. A., Dikpati, M., & Miesch, M. S. 2007, ApJS, 170, 203Gilman, P. A., & Fox, P. A. 1997, ApJ, 484, 439Goossens, M., Biront, D., & Tayler, R. J. 1981, Ap&SS, 75, 521Hathaway, D. H., Upton, L., & Colegrove, O. 2013, Science, 342, 1217Howe, R. 2009, Living Rev. Sol. Phys., 6, 1K¨apyl¨a, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22Masada, Y. 2011, MNRAS, 411, L26Masada, Y., Yamada, K., & Kageyama, A. 2013, ApJ, 778, 11Mathis, S., & Zahn, J.-P. 2005, A&A, 440, 653McIntosh, S. W., Wang, X., Leamon, R. J., & Scherrer, P. H. 2014, ApJ, 784, L32Miesch, M. S. 2007, ApJ, 658, L131Miesch, M. S., & Toomre, J. 2009, Annu. Rev. Fluid Mech., 41, 317Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2013, ApJ, 762, 73Ossendrijver, M. A. J. H. 2000a, A&A, 359, 364 36 –—. 2000b, A&A, 359, 1205Parfrey, K. P., & Menou, K. 2007, ApJ, 667, L207Parker, E. N. 1955a, ApJ, 122, 293—. 1955b, ApJ, 121, 491Passos, D., & Charbonneau, P. 2014, A&A, 568, 113Pitts, E., & Tayler, R. J. 1985, MNRAS, 216, 139Prusa, J. M., & Smolarkiewicz, P. 2003, J. Comp. Phys., 190, 601Prusa, J. M., Smolarkiewicz, P., & Wyszogrodzki, A. A. 2008, Computers & Fluids, 37, 1193Racine, ´E., Charbonneau, P., Ghizaru, M., Bouchat, A., & Smolarkiewicz, P. 2011, ApJ,735, 46Rieutord, M. 1987, GApFD, 39, 163Rogers, T. M. 2011, ApJ, 735, 100Schmitt, D., & Sch¨ussler, M. 1989, A&A, 223, 343Simard, C., Charbonneau, P., & Bouchat, A. 2013, ApJ, 768, 16Smolarkiewicz, P., & Charbonneau, P. 2013, J. Comp. Phys., 236, 608Spruit, H. C. 1999, A&A, 349, 189—. 2002, A&A, 381, 923Steenbeck, M., & Krause, F. 1969, Astro. Nach., 291, 49Stenflo, J. O., & Kosovichev, A. G. 2012, ApJ, 745, 129Strugarek, A., Brun, A. S., Mathis, S., & Sarazin, Y. 2013, ApJ, 764, 189Tassoul, J.-L. 2000, Stellar Rotation (Cambridge University Press)Tayler, R. J. 1973, MNRAS, 161, 365Tobias, S. M., Brummell, N. H., Clune, T. L., & Toomre, J. 2001, ApJ, 549, 1183Watson, M. 1981, GApFD, 16, 285