Issue 
A&A
Volume 589, May 2016



Article Number  A56  
Number of page(s)  24  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201527002  
Published online  13 April 2016 
Multiple dynamo modes as a mechanism for longterm solar activity variations
^{1}
ReSoLVE Centre of Excellence, Department of Computer ScienceAalto
University,
PO Box 15400
00076
Aalto
Finland
email: maarit.kapyla@aalto.fi
^{2} Department of Physics, Gustaf Hällströmin katu 2a (PO Box
64), 00014 University of Helsinki, Finland
^{3}
Nordita, KTH Royal Institute of Technology and Stockholm
University, Roslagstullsbacken
23, 10691
Stockholm,
Sweden
^{4}
Department of Astronomy, AlbaNova University Center, Stockholm
University, 10691
Stockholm,
Sweden
^{5}
JILA and Department of Astrophysical and Planetary Sciences, Box
440, University of Colorado, Boulder, CO
80303,
USA
^{6}
Laboratory for Atmospheric and Space Physics,
3665 Discovery Drive,
Boulder, CO
80303,
USA
^{7}
MaxPlanckInstitut für Sonnensystemforschung,
JustusvonLiebigWeg 3,
37077
Göttingen,
Germany
^{8}
Tartu Observatory, 61602
Tõravere,
Estonia
Received: 20 July 2015
Accepted: 19 January 2016
Context. Solar magnetic activity shows both smooth secular changes, such as the modern Grand Maximum, and quite abrupt drops that are denoted as grand minima, such as the Maunder Minimum. Direct numerical simulations (DNS) of convectiondriven dynamos offer one way of examining the mechanisms behind these events.
Aims. In this work, we analyze a solution of a solarlike DNS that was evolved for roughly 80 magnetic cycles of 4.9 years and where epochs of irregular behavior are detected. The emphasis of our analysis is to find physical causes for such behavior.
Methods. The DNS employed is a semiglobal (wedgeshaped) magnetoconvection model. For the data analysis tasks we use Ensemble Empirical Mode Decomposition and phase dispersion methods, as they are well suited for analyzing cyclic (nonperiodic) signals.
Results. A special property of the DNS is the existence of multiple dynamo modes at different depths and latitudes. The dominant mode is solarlike (equatorward migration at low latitudes and poleward at high latitudes). This mode is accompanied by a higher frequency mode near the surface and at low latitudes, showing poleward migration, and a lowfrequency mode at the bottom of the convection zone. The lowfrequency mode is almost purely antisymmetric with respect to the equator, while the dominant mode has strongly fluctuating mixed parity. The overall behavior of the dynamo solution is extremely complex, exhibiting variable cycle lengths, epochs of disturbed and even ceased surface activity, and strong shortterm hemispherical asymmetries. Surprisingly, the most prominent suppressed surface activity epoch is actually a global magnetic energy maximum; during this epoch the bottom toroidal magnetic field obtains a maximum, demonstrating that the interpretation of grand minimatype events is nontrivial. The hemispherical asymmetries are seen only in the magnetic field, while the velocity field exhibits considerably weaker asymmetry.
Conclusions. We interpret the overall irregular behavior as being due to the interplay of the different dynamo modes showing different equatorial symmetries, especially the smoother part of the irregular variations being related to the variations of the mode strengths, evolving with different and variable cycle lengths. The abrupt lowactivity epoch in the dominant dynamo mode near the surface is related to a strong maximum of the bottom toroidal field strength, which causes abrupt disturbances especially in the differential rotation profile via the suppression of the Reynolds stresses.
Key words: convection / turbulence / dynamo / Sun: magnetic fields / Sun: activity / stars: activity
© ESO, 2016
1. Introduction
Solar activity manifests itself through the wellknown 11year sunspot cycle, where sunspots are formed progressively closer to the equator. This is best seen in the timelatitude domain resulting in the socalled butterfly diagram. The cycle is not strictly periodic: both cycle length and amplitude are known to be variable over time. While there are signs of a longterm cyclic modulation on a centennial time scale referred to as the Gleissberg cycle, also frequent grand minima with very low activity indicators are known in solar history. The Maunder Minimum (1645−1715) and the Dalton Minimum (1790−1830) are two prime examples of such minima in existing historical sunspot records. Several such events have also been deduced indirectly, from cosmogenic radionucleid data over several millennia (e.g., Usoskin et al. 2007). Especially the actual duration and level of activity of the Maunder Minimum (hereafter MM) continues to raise questions. For example, Zolotova & Ponyavin (2015) claimed to have found historical evidence showing that some observers did not mark down all the sunspot observations on purpose because of the influence of religious or philosophical dogmas, resulting in the underestimation of spottedness during the MM, which − according to their interpretation − was actually a rather typical cyclic activity during a regular minimum of the centennial Gleissberg cycle. In the light of all the other available activity indicators (auroral sightings, cosmogenic radionuclide data, solar eclipse observations) analyzed by Usoskin et al. (2015), this interpretation seems unlikely. Moreover, during the MM, the latitude range where sunspots appeared (i.e., the width of the butterfly wings) was narrower (e.g., Ribes & NesmeRibes 1993; Ivanov & Miletsky 2011; Usoskin et al. 2015) and strong hemispherical asymmetry was present (e.g., Ribes & NesmeRibes 1993; Sokoloff & NesmeRibes 1994). Analysis of sunspot proper motion seems to indicate slower rotation and stronger latitudinal surface differential rotation during the MM than for more prominent activity states (Eddy et al. 1976; Ribes & NesmeRibes 1993).
The solar magnetic field is thought to arise as an interplay of rotation, nonuniformities related to it, and the collective inductive effects of smallscale convective turbulent motions that amplify and sustain the magnetic field against intense turbulent mixing (see, e.g., Ossendrijver 2003; Charbonneau 2014, and references therein). The classical hydromagnetic dynamo picture relies on significant turbulence effects throughout the convection zone described by tensorial turbulent transport coefficients accounting for, for example, the α effect, turbulent pumping, and turbulent diffusion (e.g., Moffatt 1978; Krause & Rädler 1980). Meanfield models of these socalled distributed dynamos usually employ subsets of turbulent transport coefficients that are either analytically derived (e.g., Pipin & Seehafer 2009) or extracted from local convection simulations (Käpylä et al. 2006). The other dynamo paradigm is the socalled fluxtransport scenario which relies on the existence of highly localized field generation regions. One such region is at the bottom of the convection zone in the tachocline and the other is at the top. In the top layer the twist of the buoyantly rising flux tubes due to the Coriolis force generates a poloidal field from the underlying toroidal field (the BabcockLeighton mechanism). The two regions are connected through turbulent diffusion and meridional flow acting as a conveyor belt (e.g., Choudhuri et al. 1995; Dikpati & Charbonneau 1999; Karak et al. 2014). Both of these scenarios are capable of explaining the regular part of the solar cycle, manifested by the equatorial symmetry properties and the migration direction of surface magnetic fields, provided the turbulence effects are suitably parameterized (e.g., Charbonneau 2010).
To explain the grand minimatype events with meanfield dynamo models, one usually invokes fluctuations in the main dynamo drivers, i.e., differential rotation, meridional circulation, and the smallscale turbulence effects. By including such fluctuations, dynamo models are usually found to exhibit irregular behavior reminiscent of grand minimatype events (e.g., Ossendrijver et al. 1996; Moss et al. 2008; Choudhuri & Karak 2012). The direct observational knowledge of the change of these quantities during a grand minimumtype event, however, is very limited, and therefore the meanfield modeling approach is problematic. Another possibility is to seek answers from direct numerical simulations of turbulent convection either in local (e.g., Käpylä et al. 2013a; Masada & Sano 2014) or global domains (e.g., Ghizaru et al. 2010; Käpylä et al. 2012; Augustson et al. 2015; Fan & Fang 2014; Mabuchi et al. 2015; Simitev et al. 2015). This is particularly promising in the latter case, where it is possible to directly track the change of all relevant dynamo drivers, provided that a desired type of dynamo solution is found. Important exceptions are the turbulent quantities, namely the inductive and diffusive parts of the mean turbulent electromotive force that require special techniques to be properly separated. One such method is the socalled testfield method (Schrinner et al. 2005, 2007, 2012), the proper application of which to solarlike global, spherical magnetoconvection solutions has recently been performed (Warnecke et al. 2016).
Oscillatory dynamo solutions from global magnetoconvection models have been known for a long time (Gilman 1983; Glatzmaier 1985). The first solarlike solutions, however, were obtained only recently (Schrinner et al. 2011; Käpylä et al. 2012; Augustson et al. 2015), and only a couple of them have been run up to time scales of interest for detecting the irregular variations: the EULAG code Millennium simulation at solar rotation (Passos & Charbonneau 2014; Norton et al. 2014, hereafter EULAGMHD) covers roughly 20 magnetic cycles of 80 years cycle length, while the ASH code simulation (Augustson et al. 2015, hereafter ABMT), rotating at three times the solar rate, covers roughly 24 cycles with a cycle length of 6.2 years. The former does not exhibit significant irregular behavior and only very weak latitudinal migration of the magnetic field, while the latter produces a clearer equatorward migrating branch at lower latitudes and a poleward migrating branch at higher latitudes, especially pronounced in the radial field. In addition, ABMT report a particular grand minimumtype event, where magnetic activity is suppressed in the equatorial surface region and polarity reversals are not seen in the polar surface regions for roughly five half cycles. The interpretation of the origin of an oscillatory magnetic field and equatorward migration varies considerably: from a turbulent dynamo picture producing an αΩ oscillatory solution exhibiting a latitudinal dynamo wave (Käpylä et al. 2012; Warnecke et al. 2014) to the magnetic field feeding back to differential rotation via the Lorentzforce, the field migration being related to the variations in the differential rotation and possibly to a nonlinear dynamo wave (Augustson et al. 2015). However, the latter interpretation might suffer from an inaccurate calculation of the α effect; see Warnecke et al. (2016).
In this paper we analyze a semiglobal (wedgeshaped) magnetoconvection simulation similar to those reported earlier (Käpylä et al. 2012; Warnecke et al. 2014) with slightly varied parameters, but integrated over a much longer time span. The obtained simulation shows solarlike migration patterns of the toroidal magnetic field, and exhibits a dominant cyclic dynamo mode with an average cycle length of roughly 5 years. The simulation covers roughly 80 such cycles. In addition to this basic mode, two other significant modes are detected and characterized with suitable time series analysis tools designed for nonperiodic signals: Ensemble Empirical Mode Decomposition (Wu & Huang 2009) and phase dispersion statistics (Pelt 1983). In the following, we refer to these methods as EEMD and D^{2} statistics, respectively. One of the main goals of this paper is to investigate the significance of these multiple dynamo modes to the global dynamics of the system. The solution also exhibits several epochs of abrupt irregular behavior (disappearance of the surface activity, sudden switches of parity, sudden changes in cycle length and migration direction of the toroidal field), and a physical explanation for these events is sought by computing proxies for the different dynamo drivers, namely differential rotation, meridional circulation and the α effect, during different activity states.
2. Model
Our magnetohydrodynamic (MHD) model has been described in many earlier studies, in particular in Käpylä et al. (2013b), and the details will not be repeated here. We perform computations in a spherical wedge using spherical polar coordinates (r,θ,φ) corresponding to radius, colatitude, and longitude, respectively. The computational domain spans from r_{0} = 0.7 R_{⊙} to r_{1} = R_{⊙} in the radial direction, where R_{⊙} = 7 × 10^{8} m is the solar radius; from θ_{0} = π/ 12 to θ_{1} = 11π/ 12 in colatitude (±75° latitude); and 90° in longitude, i.e., a quarter of a full sphere. Our setup is semiglobal in the sense that we exclude the polar regions. Including the poles would require prohibitively short timesteps. Within this domain, we solve the standard compressible magnetohydrodynamic equations for logarithmic density lnρ, specific entropy s, velocity u, and magnetic vector potential A, which gives the magnetic field as B = ∇ × A. To close the system of equations, we assume that the fluid obeys the ideal gas law with p = (γ − 1)ρe, where γ = c_{P}/c_{V} = 5 / 3 is the ratio of specific heats at constant pressure and volume, respectively, and e = c_{V}T is the specific internal energy, where T is temperature. The fluid is subject to gravitational acceleration, g = −GM_{⊙}r/r^{3}, where G = 6.67 × 10^{11} m^{3} kg^{1} s^{2} is the gravitational constant and M_{⊙} = 2.0 × 10^{30} kg is the mass of the Sun, and to rotation, the rotation vector being Ω_{0} = (cosθ, − sinθ,0)Ω_{0}. We neglect selfgravity of the gas in the convection zone.
As explained in detail in Käpylä et al. (2013b), since we cannot get anywhere near the high Rayleigh numbers of real stars, we must use higher diffusivities. In the present model this implies a roughly 10^{6} times higher luminosity in the model in comparison to the Sun. This allows us to reach the KelvinHelmholtz time scale in our simulations, implying that our runs are thermally relaxed. As the convective energy flux scales as F_{conv} ~ ρu^{3}, the convective velocity u is roughly 100 times greater in the simulations than in the Sun. To obtain the same rotational influence on the flow as in the Sun, we must therefore increase Ω by the same factor. In general, the scaling of velocity and rotation rate can be written as (1)where L_{ratio} = L_{0}/L_{⊙}, with L_{0} and L_{⊙} ≈ 3.84 × 10^{26} W being the luminosities of the model and the Sun, respectively. In the current model, as in Warnecke et al. (2014), we then renormalize our simulation to the rotation rate Ω_{0} ≡ 5Ω_{⊙}, where Ω_{⊙} ≈ 2.7 × 10^{6} s^{1} is the mean solar rotation rate, corresponding to Ω_{⊙}/ 2π = 430 nHz.
In what follows we express our results in solar units so that for the velocity, for example, we quote . The scaling used here is based on dimensional arguments. It is supported by mixing length theory (Vitense 1953) and simulations (Brandenburg et al. 2005; Miesch et al. 2012), and should be applicable as long as the energy transport is not yet affected by rotation (see, e.g., Yadav et al. 2013). Furthermore, we assume that the density at the base of the convection zone (r = 0.7R_{⊙}) has the solar value ρ = 200 kg m^{3}.
The simulations are performed with the Pencil Code^{1}, which uses a highorder finite difference method for solving the compressible equations of magnetohydrodynamics.
2.1. Initial and boundary conditions
The initial hydrostatic state is described by a polytrope with index n_{ad} = 1.5, i.e., an isentropic stratification. Our model does not include a stable overshoot region in the bottom of the convection zone. The density stratification is roughly 30 in the beginning, while it decreases to about 20 in the course of the simulation. In the final state, the number of density scale heights covered by the model is roughly 3. To speed up the thermal relaxation, the initial condition is chosen not to be in thermostatic equilibrium, but closer to the final convecting state. We choose the heat conduction profile K(r) such that it decreases toward the surface as r^{15}. This means that a very small portion of the energy is carried by radiative diffusion near the surface (Brandenburg et al. 2005; Käpylä et al. 2014). We introduce weak smallscale Gaussian distributed noise as velocity and magnetic field perturbations in the initial state.
Our simulations are defined through the energy flux imposed at the bottom boundary, F_{b} = −(K∂T/∂r)  _{r = r0}; the values of the rotation rate Ω_{0}; the kinematic viscosity ν; the magnetic diffusivity η; and the subgrid scale diffusivity for the entropy in the middle of the convection zone, . The radial profile of χ_{SGS} is similar to that shown in Fig. 1 of Käpylä et al. (2011) with the surface value being , corresponding to 6.2 × 10^{8} m^{2} s^{1} in physical units. The control parameters of the simulation are quantified by the thermal and magnetic Prandtl numbers, and Pr_{M} = ν/η, respectively, and the Rayleigh number , where the subscript “hs” indicates a hydrostatic nonconvecting state and Δr = r_{1} − r_{0} = 0.3 R_{⊙} is the depth of the convection zone.
The radial and latitudinal boundaries are assumed to be impenetrable and stress free for the flow; see Eqs. (8), (9) of Käpylä et al. (2013b). The magnetic field is assumed radial on the outer boundary at r_{1}, while the latitudinal and inner radial boundaries are perfectly conducting; see Eqs. (10)−(12) of Käpylä et al. (2013b). Density and specific entropy have vanishing first derivatives on the latitudinal boundaries. On the outer radial boundary we apply a black body condition σT^{4} = −K∇_{r}T − χ_{SGS}ρT∇_{r}s, where σ is a modified Stefan−Boltzmann constant, which is chosen such that the flux at the surface carries the total luminosity through the boundary in the initial nonconvecting state.
Input and diagnostic parameters together with the kinetic and magnetic energies realized in the simulation in units of 10^{5} J m^{3}.
2.2. Diagnostic quantities
The most important nondimensional diagnostic quantities of our model are the fluid and magnetic Reynolds numbers, Re = u_{rms}/νk_{f} and Re_{M} = u_{rms}/ηk_{f}, where is an estimate of the full threedimensional rms velocity based on the values in the meridional plane, and the subscripts indicate averaging over r, θ, φ, and a time interval where the run is thermally relaxed, and is an estimate of the wavenumber of the largest eddies. The azimuthal velocity component is not included in the computation because it is dominated by the differential rotation, which would then not characterize the level of turbulence, which is what we are interested in. We also use the meridional distribution of turbulent velocities , where the fluctuating velocity is defined via , where the overbar denotes an azimuthal average (see Sect. 3 for a more specific definition of the averages). Furthermore, the rotational influence on the flow is measured by the Coriolis number Co = 2Ω_{0}/u_{rms}k_{f}.
Our simulations are characterized by Re = Re_{M} ≈ 30, Pr_{SGS} = Pr_{M} = 1, Co ≈ 9.5, and Ra ≈ 1.0 × 10^{7}. The corresponding numbers computed for the total velocity field, i.e., including the azimuthal velocity due to differential rotation, read Re_{tot} = Re_{M}_{tot} ≈ 55 and Co_{tot} ≈ 5.1. The most important input parameters and diagnostic quantities are also listed in Table 1.
This run has therefore a smaller Pr_{SGS} than the runs with equatorward migrating magnetic fields of Käpylä et al. (2012) and Warnecke et al. (2014, 2015), but twice the values of Pr_{SGS} and Pr_{M} than the run with poleward migrating field of Warnecke et al. (2014) by leaving the other parameters the same. The simulation of Augustson et al. (2015) has 2−4 times smaller Prandtl numbers, but otherwise comparable parameters. The EULAGMillennium simulation setup is similar to that described in Ghizaru et al. (2010) and Racine et al. (2011), but with some explicit diffusion added for stability reasons. The estimated values of their Re and Re_{M} are in the range 30−60, and the magnetic Prandtl number is of the order of unity.
3. Data analysis
Because MHD processes are nonlinear and the resulting data nonstationary, it is necessary to choose analysis tools that accurately describe its cyclic components locally and adaptively. For example, the solar cycle is known to be of varying length and amplitude. While Fourier analysis is the most common data analysis technique used to extract periodicities from harmonic signals, it requires constant amplitudes and phases and is not well suited to the problem (Barnhart 2011). In this work we have chosen to utilize the EEMD and D^{2} statistics, both of which are suited to the analysis of nonstationary signals. They are presented in Appendices A and B, respectively.
We define mean quantities as azimuthal averages (over the φcoordinate) and denote them by overbars. Fluctuations about the mean are denoted by a prime, e.g., . We also often average the data in time over the period of the simulations where thermal energy and differential rotation have reached statistically saturated states. To study the hemispherical asymmetries, we analyze azimuthally averaged quantities separately for the northern and southern hemispheres. We also often employ the toroidal vs. poloidal decomposition for the azimuthally averaged, and therefore axisymmetric, magnetic field, i.e., (2)with , where ê_{φ} is the unit vector in the azimuthal direction, and where the radial and latitudinal components form the poloidal component . A similar decomposition is used for the velocity field, , where . The three main depths at which we plot/analyze the time series are near the surface, R_{s} = 0.98 R_{⊙}; at the middle, R_{m} = 0.85 R_{⊙}; and close to the bottom, R_{b} = 0.72 R_{⊙}, of the convection zone.
Fig. 1 Time evolution of the mean toroidal magnetic field in the convection zone for three depths (from top to bottom, R_{s} = 0.98 R_{⊙}, R_{m} = 0.85 R_{⊙}, and R_{b} = 0.72 R_{⊙}). The extrema are [−11.0,10.3 ] kG at R_{s}, [−23.3,20.7 ] kG at R_{m}, and [−29.9,31.7 ] kG at R_{b}. 
Fig. 2 Evolution of the mean radial magnetic field in the convection zone for three depths (from top to bottom, R_{s}, R_{m}, and R_{b}). The color ranges have been clipped at half of the maxima/minima, i.e., [−11.7,10.4 ] kG at R_{s}, [−6.1,5.7 ] kG at R_{m}, and [−3.3,3.8 ] kG at R_{b}, to emphasize the weaker structures at low latitudes. 
Fig. 3 a) Total (, purple dashed), total mean (, black solid), toroidal (, blue), and poloidal (, red) rms values as functions of time as averages over the whole volume. b) in the north (blue) and south (red). The black solid line shows the ratio of the northern to southern energy as a function of time, and the orange line the mean of the ratio computed over the whole time span. c) The same for the poloidal field. d) The bottom (at R_{b}, blue line) and surface (at R_{s}, red line) toroidal magnetic field strengths as functions of time. The thin orange horizontal lines indicate the mean of the quantities over the whole simulation time. In a) and d) thick orange lines show the smoothed curve of the total mean magnetic field a) and the bottom/surface toroidal magnetic field strength d); the yellow horizontal lines show 1.1 and 0.9 times the mean levels, respectively. 
Fig. 4 a) Rms value of the mean azimuthal velocity as a function of time together with the average over the whole time series plotted with a black horizontal line. The dashed black line shows the volumeaveraged rms value of the fluctuating velocity field, u_{rms}. The times considered “high” are indicated in red, “low” in blue, and “nominal” activity states according to the strength of the total magnetic field in orange (D1, see Sect. 4.4). b) As in a), but for the rms value of the total magnetic field as a function of time. The black horizontal line indicates the mean over the whole time series. c) Rms value of the meridional velocity as a function of time (black solid line) with the mean over the whole time series (orange horizontal line) and a 4.9year moving average (thick red solid line). d) Zoomin to 50−65 years of evolution of the mean azimuthal velocity (black solid), mean meridional velocity (blue line, multiplied by a factor of 10 to fit the figure), and the total mean magnetic field strength (red line, scaled to fit the figure). 
4. Results
We ran our simulation for nearly 430 years in physical time units. The solution we obtain is oscillatory, but it shows very complex behavior. The definition of a mean cycle period is not unique, but − based on the nearsurface toroidal magnetic field oscillation − the mean cycle length is, in physical units, 4.9 years and therefore, on average, our simulation covers 80 cycles. In comparison to the Sun, the obtained cycle length is roughly a factor of 4 shorter. If, however, this was expressed in solar units, the length of our simulation would roughly correspond to an evolution of two millennia. From a similar simulation with equatorward migrating fields (Warnecke et al. 2015), we find that the cycle length increases with slower rotation, showing that with more realistic parameters it is possible to obtain cycle lengths comparable to the solar value. Evolving a simulation with a longer cycle over the time span presented here, however, is computationally prohibitively expensive. In comparison to EULAGMHD, the number of cycles covered here is roughly two times larger, and four times larger in comparison to ABMT. Furthermore, more pronounced solarlike latitudinal migration and irregular behavior are observed in our simulation than in EULAGMHD. In the ABMT model, a clear solarlike equatorward migration is seen with one MMtype event. Some fundamental differences also relate to the data analysis, especially concerning the EULAGMHD: in Passos & Charbonneau (2014) the analogue to the sunspot cycle is defined using the toroidal magnetic field strength at the bottom of the convection zone at high latitudes, and a poloidal field proxy is built using the radial field in the nearsurface polar regions, while our analysis covers the whole convection zone. Moreover, we restrain from making a detailed comparison to the observed statistical properties of the solar cycle; we feel that such attempts are pointless owing to the vast difference between simulated and real parameter regimes. Instead, we concentrate on the general effects that may cause the irregularities that are seen to occur in our simulation.
4.1. Overall evolution of the magnetic field
In Fig. 1 we show longitudinal averages of the toroidal magnetic field as a function of time and latitude (butterfly diagram) for three depths in the convection zone (R_{s}, R_{m}, and R_{b}), and in Fig. 2 the corresponding plot for the radial field. The dynamo solution is, in general, very similar to those reported by Käpylä et al. (2012) and Warnecke et al. (2014, 2015) with cyclic behavior consisting of an equatorwardmigrating dynamo wave at low latitudes and a polar branch at high latitudes near the surface. The dominating component near the surface is the radial one with roughly twice the strength of the azimuthal component, while the azimuthal field becomes dominant at greater depths, being roughly twice (four) times stronger than the radial field in the middle (at the bottom) of the convection zone, see also Fig. 5e in Warnecke et al. (2014) for a plot of a similar run. The lack of strong toroidal magnetic field in the top layer is also related to radial field boundary condition at the surface as investigated in detail in Warnecke et al. (2015). The marked difference to our earlier results is the irregularity and hemispheric asymmetry seen in the field evolution. During some epochs, the surface activity ceases and/or becomes disturbed (i.e., the cycle length changes significantly), which can happen asynchronously in different hemispheres. The radial field evolution is more regular than that of the toroidal field, even though the ceased surface activity is perhaps even better seen in the radial field evolution, see Figs. 1 and 2. During the most severe event (20−40 years in the south; 35−45 years in the north), both the radial and toroidal fields near the surface practically vanish at midlatitudes. Simultaneously, regions of strong toroidal magnetic field are seen in the nearequator regions (at R_{s}) and also in the bottom layers (at R_{b}). Even visually, the toroidal field at R_{b} exhibits another toroidal dynamo mode with considerably longer period than the basic cycle of 4.9 yr seen throughout the simulation. In addition to these dynamo modes, there appears to be a highfrequency (shortperiod) poleward dynamo mode in the upper part of the convection zone, confined to the nearequator region, similar to the one reported in Käpylä et al. (2012) and Warnecke et al. (2014, 2015). This dynamo mode might be related to a locally operating α^{2} dynamo in the top layers (Käpylä et al. 2013b; Warnecke et al. 2014, 2016). The quantitative characterization of these modes is presented in Sect. 4.3.
The volume and timeaveraged energies of different flow and magnetic field components are listed in Table 1. In Fig. 3a, we show the volumeaveraged rms values of the total and mean magnetic field together with its mean toroidal and poloidal components as functions of time. The dominant component is the fluctuating part containing roughly twice as much energy as the mean field. In a wedge covering only one quarter of the sphere in longitude, the largescale magnetic field is almost purely axisymmetric and contains only a negligible contribution from nonaxisymmetric modes. Of the mean fields, the dominant component is the toroidal one, containing roughly twice the energy of the poloidal component. This is consistent with the conclusion of Warnecke et al. (2014) that the dynamo is of αΩ type. The poloidal field evolution shows clear cyclicity with modest longterm amplitude variations, whereas the longterm variations of the toroidal field are much stronger. The irregularities are stronger in the early phases of the simulation; the lowered and disturbed surface activity in the butterfly diagrams of Figs. 1 and 2 are clearly seen as a global maximum of the volumeintegrated magnetic field energy, especially in the toroidal field (see Fig. 3). The early phases (about 3/4 of the simulation length) show markedly irregular behavior, while the later phases (the last 1/4) shows more regular behavior with slightly decreased overall magnetic activity level.
In panels (b) and (c) of Fig. 3, we plot the rms values of toroidal (b) and poloidal (c) components for both hemispheres (north − blue; south − red). The solid black lines in these panels show the ratio of northern to southern energy (expected to attain the value of unity if the hemispheres are totally symmetric). Although averaged over the full time series, the asymmetry measure is close to unity for both the toroidal (2.7% dominance of the northern hemisphere) and poloidal components (1.7% for the north), the variance around the mean is significant (around 0.1) for both field components, although the volumeintegrated quantities indicate no clear correlation between the asymmetry and extrema in the total energy. We will investigate the northsouth asymmetry in more detail in Sect. 4.6.
In Fig. 3d we plot the toroidal magnetic field strength near the surface (at R_{s}, red line) and the bottom (at R_{b}, blue line). The toroidal fields located at the bottom of the convection zone are roughly three times stronger than those at the surface. The surface field shows some variability, but significantly less than those associated with the bottom toroidal magnetic field. While the surface toroidal field strength shows no systematic trend when the early and late parts of the simulation are compared, the bottom field is clearly decaying in strength. Therefore, the overall decrease and change in the magnitude of the irregularities in the total mean field evolution can be attributed to the changes seen in the bottom toroidal field.
4.2. Overall evolution of the velocity field
In Fig. 4 we show the rms values of the azimuthal, fluctuating (a), and meridional (c) velocities. The system energetics is dominated by differential rotation, the energy of which comprises roughly 57 percent of the total kinetic energy, which is roughly 2.6 times stronger than the total magnetic energy. The energy contained in the meridional flow is negligibly small; in terms of the rms velocities, the meridional flow is roughly 15 times weaker than the azimuthal value. The energy of the fluctuations is roughly 1.3 times weaker than the energy of the axisymmetric part. Evidently, the azimuthal velocity is strongly affected by the dynamically significant magnetic field. A strong global magnetic field quenches the azimuthal velocity, while during weaker magnetic field epochs, the average azimuthal velocity is higher, as also indicated by the colorcoding applied in Figs. 4a and b for high, low, and nominal states defined based on the global magnetic energy level (D1; see Sect. 4.4 for a detailed definition). The average meridional flow is an order of magnitude weaker, but shows signs of weak systematic variation when a running average with the mean magnetic cycle length is applied. In Fig. 4d we show a zoomin over a time epoch (50−65 years) where no significant extrema are seen either in the total magnetic field strength or the azimuthal velocity. This figure reveals that both the azimuthal and meridional velocities show a variation over the magnetic cycle, stronger in the former, weaker (yet noticeable) in the latter. The azimuthal velocity and magnetic field are phase shifted in such a way that the azimuthal velocity grows when the magnetic field is weak, and decays when it is strong, in broad agreement with Warnecke et al. (2012) and Karak et al. (2015). The meridional velocity shows a similar trend, even though it is much more difficult to detect as the fluctuations dominate the weak signal. Therefore, in terms of the magnetic cycle, the angular velocity variations, often called torsional oscillations, occur with twice the frequency of the magnetic cycle itself (see Fig. 7 and Sect. 4.3 for a more thorough analysis). Our results are in fair agreement with those obtained for the EULAGMHD and ABMT models during the regular periods of evolution.
In contrast to the magnetic field, the mean velocity field components show only very little hemispherical asymmetry; the evolution of the quantities are virtually identical for both, and are therefore not separated for north and south in Fig. 4. It is evident that the strong, frequently occurring, shortterm asymmetries are a property of the magnetic field alone and any asymmetry in the flow is probably caused by the asymmetry in the magnetic field. In the top panel of Fig. 7 we plot the timelatitude diagram of the angular velocity variations at R_{s}. From this plot it can be seen that in addition to the rather weak, regular oscillation with twice the magnetic cycle frequency, which is evident at all depths, there is a rather strong longterm variation of the angular velocity near the equator in the top layers. The origin of this variation is analyzed in detail in Sect. 4.5.1.
4.3. Multiple dynamo modes and their significance
The results presented in this section were obtained by applying EEMD analysis (see Appendix A) on the azimuthally averaged quantities using an ensemble size of 100 and Gaussian white noise with standard deviation equal to that of the time series being analyzed. The time series were chosen by sampling the simulation domain with 12 points in the radial direction and with 18 points in the latitudinal direction. Our analysis is focused on the components of mean magnetic field, mean velocity field, and kinetic helicity; we aim to extract the most significant modes for the given time series globally and at three depths (R_{b}, R_{m}, and R_{s}). We also use the D^{2} statistic (see Appendix B) on some selected data sets to validate the EEMD results, but also to investigate the coherence of the cycles over time.
Summary of the mode quantities.
On average we detected a total number of 13 modes. As expected, it turned out that most of the energy was contained in a limited number of modes. In the following we use a mode numbering scheme such that the modes are ordered by their average frequency from higher to lower values (the highest frequency has an index of 1). To express the energy contained in each mode quantitatively we use the notion defined as (3)where A is the physical quantity analyzed, i is the mode index, is the mean square amplitude of the mode, and summation is performed over all modes and integration either over latitude at selected radii or over the full meridional plane.
Fig. 5 Distribution of mean amplitudes of the most significant modes of mean magnetic and mean velocity fields found from EEMD. The figures are labeled with the physical variable followed by mode index (e.g. 7 indicates mode 7 of the mean latitudinal magnetic field). Colors reflect the mean amplitude (blue − high, yellow − low) of the mode at the given location; we note that the scales of separate figures are different. Contours on the plot represent the lines of constant amplitude with values given in the legend. 
Fig. 6 Phase dispersion analysis results for the components of . The calculations were done for at latitude ± 22° and radius 0.94 R_{⊙}, for at latitude ± 66° and radius 0.94 R_{⊙}, and for at latitude ± 49° and radius 0.82 R_{⊙}. Left (right) refers to the northern (southern) hemisphere. 
However, calculating only the energy content does not allow us to detect weak, but more regular (harmonic), modes. For that purpose, for each mode we calculate the spectral entropy H = ^{∫}p(ν)log p(ν)dν, where p(ν) is the power spectral density at frequency ν. The average value of the given quantity per each mode compared to the expected value in the case of white noise serves as a regularity measure of the mode (i.e., how much more Gaussian or nonuniform the spectrum is compared to the spectrum of the white noise). An alternative, but more straightforwardly calculated measure is the ratio of energies in each mode to the expected values of white noise. A large value of this ratio for a given mode is another indication of the presence of a more regular signal. This idea is captured by the formula (4)where E_{i} is a mode energy fraction defined in Eq. (3) and is the expected energy fraction of white noise for the given mode. As is explained later, the results for both regularity measures are in good agreement with each other and so in the following we only show quantitative values calculated based on the second definition.
The results of the EEMD analysis are gathered in Table 2, where the mode period represents a descriptive period calculated based on the zerocrossings of the most dominant instance of the mode of given order (i.e., the intrinsic mode function, hereafter IMF, with the highest average mean energy over the full meridional plane). In the table we also show fractions E_{i} of the energy contained in the ith mode calculated over the full meridional plane and over the latitude at radii R_{b}, R_{m}, and R_{s}. The latitude values are approximately equal to the reconstruction error of the initial signal given that the ith mode would be omitted from the sum. For clarity, only the values above 10% are shown and lower values have been marked with “−”. We can see that in the case of the magnetic field components most of the energy is contained in mode 7 (hereafter M7), which can be associated with the dominant cyclicity of the 4.9year cycle length deduced from D^{2} statistic using surface toroidal field. Only for in the bottom layer is the energy more spread between modes with a longer period. In Fig. 5 the two upper rows show the spatial distribution of the mean amplitudes of the most significant modes found in the magnetic field; it is clear that for the toroidal component of the magnetic field, M7 is mainly confined to the region where the equatorward migration of the field is observed (compare Fig. 5, second panel from left in the upper row, and Fig. 12, top panel). For the poloidal field, however, the M7 amplitude is the highest at higher latitudes: for the radial field near the surface close to the latitudinal boundaries, and for the latitudinal field near the bottom and high latitudes. The modes with shorter cycles (mode 1, hereafter M1, contains most of the energy) have the highest amplitudes near the equatorial region and near the latitudinal boundaries, both regions being concentrated near the surface. The long modes (the most notable is mode 11, hereafter M11, which also has an increased regularity measure) reside at the very bottom of the convection zone, confined mainly to lowlatitude regions.
In the case of , the energy spectrum is flatter with roughly equal distribution of energy between modes 6−10; see Table 2. From the third row of Fig. 5 it is evident that the highest amplitudes of these modes all occur in the equatorial region, the region extending to higher latitudes for modes 6−7, and getting narrower for the ones with a longer cycle. For the other velocity components − as well as for kinetic helicity defined as , where ω′ = ∇ × u′− energy is primarily contained only in modes 1 and 2. From the last row of Fig. 5 we can see that these short cycles are most prominent near the surface close to the latitudinal boundaries for the radial component of the velocity, and near the equator at the surface for the latitudinal velocity component. The modes in the kinetic helicity show high amplitudes near the equatorial region extending deeper to the convection zone, but also in a narrow band near the surface over the whole latitudinal extent.
In the penultimate column of Table 2, we show a measure of mode regularity as defined in Eq. (4). These values are calculated at those radii and latitudes where the energy of the given mode has its maximum. Higher values of this number indicate higher energy content in the given mode compared to the white noise level. According to the values shown, in addition to M7 of the magnetic field, we have one additional mode where the regularity clearly has a high value. This is M11, which has, however, a significant energy fraction only in the case of in the bottom layer. It is also quite obvious that the leading modes of , , and H_{kin} represent primarily noise, because the value of regularity in these cases is very low. It is noteworthy, however, that for all these quantities the regularity measure of M11 is high, indicating the presence of a longterm regular cycle. Here we also note that the mean spectral entropy calculations lead to similar results to those above: the spectrum of M7 of the magnetic field is the most regular one, while modes 10 and 11 are less regular, although significant deviations from the entropy of the white noise spectrum are still detected. It can be concluded, therefore, that the dynamo mode at the bottom of the convection zone has an influence on the overall dynamics of the system although it is insignificant in terms of the total energetics of the system.
The last column of Table 2 represents the equatorial symmetry of the modes, calculated as a time averaged parity defined by Eq. (13). The EEMD analysis confirms our earlier conclusion that the mean azimuthal velocity field is mostly symmetric between the hemispheres because the parity is high. For all the magnetic field modes, however, values significantly deviating from perfect antisymmetry (parity values of − 1) or perfect symmetry (parity values of + 1) are found. The higher frequency modes show parity values fluctuating around zero, while the lower frequency modes show an increasing tendency for antisymmetric parity as the period of the cycle increases. The equatorial symmetry will be discussed in more detail in Sect. 4.6.
As an independent check for the EEMD results we choose the D^{2} phase dispersion statistic (see Appendix B). Based on the spatial distribution of the most significant modes, we calculate the D^{2} spectra at selected depths for the components of . Ideally, we would expect to see a match between the average mode periods found from EEMD and cycle lengths found from D^{2}.
To focus on the most pronounced activity regions we chose latitude ±22° in the layer at R_{s} for , latitude ± 66° in the layer at R_{s} for , and latitude ± 49° in the layer at R_{m} for to detect the short cycles such as M7 (see Table 2). Owing to the limited length of the data set, detecting the cycle is not feasible with the D^{2} analysis. However, by using lowpass filtering on the raw data, the existence of M11 can be easily confirmed. The results of the D^{2} analysis are depicted in Fig. 6. Evidence of M7 with a cycle length of around 5 years can be seen for all components. It is interesting, however, that the cycle is modulated differently for the northern and southern hemispheres. For the southern hemisphere there is an indication of amplitude modulation with a long period, while for the northern hemisphere the modulation seems to be more complex. If we assume a simple amplitude modulation for the southern hemisphere, and this is justified because we see a relatively symmetric split of the spectrum when moving from short coherence times to longer ones, then by measuring the difference of these peaks from the spectrum taken at the high end of the coherence time, we can estimate the period of this modulation. Given the low precision that we can achieve using this simple procedure we obtained the value for the modulating period to be around 100 ± 10 yr.
We also analyzed the higher frequency region (0.1 to 2 yr) with the D^{2} statistic, but did not detect any strong and/or singular minima in this range. Given that the regularity measures obtained from EEMD also indicate very low values for these modes, we have to conclude that this cyclic mode is coherent only over very short time intervals; the variation of the length of this cycle is comparable to the length of the cycle itself. Therefore, the signal of this cycle can be spread between multiple EEMD modes, and the spectrum broadened over a large frequency band.
4.4. Definition of the activity levels
As was evident from the results presented in Sect. 4.1, a decreased surface magnetic activity does not directly imply a global magnetic energy minimum, especially during those epochs when the longperiod dynamo mode at the bottom of the convection zone is strong (roughly speaking the first 3/4 of the simulation). Therefore, the definition of a low/high state is not obvious. Here we use three different definitions, all based on computing running averages with a window coinciding with the main periodicity present at different depths and latitudes (4.9 years):

D1.
If the smoothed volumeaveraged magnetic energy is higher (lower) than 1.1 (0.9) times its average over the whole simulation, the state is termed high (low). All other epochs are termed nominal states of activity; see Fig. 3.

D2.
Otherwise identical, but the surface magnetic field energy at R_{s} in comparison to its smoothed value is used as the criterion.

D3.
Otherwise identical, but the energy contained in the magnetic field mode of the bottom layer at R_{b} is used as the criterion.
In Fig. 3(a; D1) and (d; D2 and D3) we demonstrate these criteria by overplotting the smoothed curves (orange thick lines) and the corresponding upper/lower limits (yellow horizontal lines) on top of the original data. In the remaining parts of this section, we use these definitions to analyze some key quantities during the different activity states. In practice this is done by computing the mean of a quantity over the whole simulation run, dividing the data points into different states according to each criterion, and computing the mean profiles for the three activity states separately. If the profiles obtained in this way differ significantly from the global mean, i.e., the difference is larger than the standard deviation of the quantity, we take note of this dependence.
Fig. 7 From top to bottom: timelatitude diagram of the angular velocity variations Ω^{v}, variations of the largescale Lorentz force , and the Reynolds stress components R_{rφ} and R_{θφ} near the surface (at R_{s}). For the Reynolds stress plots, the contours have been clipped at half of the maxima/minima. 
4.5. Variation of the dynamo numbers
In this section, we investigate the variability of the three key factors involved in the dynamo process, namely differential rotation and meridional circulation (Sect. 4.5.1), as well as the inductive effects due to convective turbulence (α effect; through kinetic and magnetic helicities as its proxy, Sect. 4.5.2), using the different definitions of the activity level.
4.5.1. Rotation and its nonuniformities
As stated in Sect. 2, the input angular rotation velocity is five times the solar rotation rate. The differential rotation profile generated in the simulation is solarlike, although the equator is rotating only approximately 7% faster than the poles when the magnetic fields have grown and saturated at dynamically significant strengths. Therefore, the obtained latitudinal differential rotation is approximately three times weaker than in the Sun. We have also performed a hydrodynamical counterpart simulation, in which the latitudinal differential rotation is found to be two times stronger (poleequator difference of 14%). The distribution of angular velocity is fairly similar to the MHD state, discussed later in this section, and therefore the HD results are not shown separately. Furthermore, the hydrodynamic state is steady, and shows no oscillatory behavior similar to the dynamo run. Therefore, the variations seen in the MHD state (Fig. 4) arise as a consequence of the dynamically significant magnetic field to the flow field. Such a backreaction can occur through two pathways: a largescale Lorentzforce (the socalled MalkusProctor effect; see Malkus & Proctor 1975) or through turbulence effects either directly on the convective motions and thereby the generators of differential rotation (the Λeffect: see Rüdiger 1989) as recently explored in Karak et al. (2015), or through the turbulent Maxwell stress (e.g., Käpylä et al. 2004). At the same time as it quenches the flow field, the magnetic field suppresses its own generators. One argument to explain both the regular and irregular parts of the solar cycle is related to this mechanism (see, e.g., ABMT). The classical theory of turbulent hydromagnetic dynamos, however, allows for the excitation of oscillatory solutions independent of the preexistence of such modulation in the flow field (see, e.g., Parker 1955). In this section we investigate the role of the changes seen in the mean flow to the longterm evolution of the magnetic field, especially to its disturbed states.
We begin by investigating the angular velocity variations Ω^{v} over time in more detail. Here Ω^{v} = Ω− ⟨ Ω ⟩ _{t}, where and ⟨ Ω ⟩ _{t} is timeaveraged over the saturated stage. We plot a timelatitude diagram of this quantity at R_{s} (Fig. 7, top panel), at which depth the most prominent variations occur in the equatorial regions. In addition to the strongest, seemingly irregular changes symmetrically distributed around the equator, a weaker modulation corresponding to the poleward migrating magnetic cycle can be seen at high latitudes throughout the simulation. In Fig. 7, second panel, we investigate the role of the largescale Lorentz force in causing the changes seen in the angular velocity. It is evident that the highlatitude variations in the angular velocity, following the mean magnetic cycle, are mediated through this effect. The sudden drops of the surface magnetic field strength are very clearly seen as similar behavior in the largescale Lorentz force and occur simultaneously with the surface field disappearances, especially pronounced during the first 50 years of the simulation. Only very mild enhancements of differential rotation at higher latitudes are accompanied by the sudden drops in the Lorentz force. The strong variations seen in the angular velocity near the equator, on the other hand, cannot be explained by the largescale Lorentz force.
Their symmetric character with respect to the equator hints that they are related to the Reynolds stress component . Instead of inspecting all the Reynolds stress components, we concentrate here on examining the variability of those known to be the most significant for the generation of differential rotation, namely R_{rφ} for vertical and for horizontal differential rotation, and postpone the full analysis of the turbulent quantities to a forthcoming paper. The timelatitude plot of the aforementioned Reynolds stress components are shown in the two lowermost panels of Fig. 7. As is evident from these figures, both stress components undergo variations on the time scale of the shortest dynamo cycle (M1), while showing no modulation by the dominant magnetic cycle (M7). The stresses have nonzero values only near the equatorial regions, where the strongest angular velocity variations are seen. The minima/maxima of the stresses coincide very well with the minima/maxima of the angular velocity and therefore indicate a strong connection between these quantities.
Fig. 8 Top row, from left to right: normalized timeaveraged rotation profile ⟨ Ω ⟩ _{t}/ Ω_{0}, and difference to the rotation profile ΔΩ_{state} during the high, low, and nominal state of global activity (D1). Middle (bottom) row, the same for mean radial (latitudinal) differential rotation. 
Fig. 9 Meridional flow, indicated as arrows, and timeaveraged dynamo parameter C_{u}, indicated with color contours. On the left, we show the meridional circulation profile averaged over the whole simulation time span, and the next three panels show the difference ΔC_{u, ⟨ state ⟩} = C_{u} − C_{u, ⟨ state ⟩} to the mean during a high, low, and nominal global activitystate (D1). 
Next we apply the three different definitions of high and low magnetic activity states (D1−D3) introduced above, and look for significant deviations, indicative of the sensitivity of the mean flow field on the quantity chosen. The mean flow is almost completely insensitive to the variations of the surface magnetic activity level (D2), and not significantly different when the bottom activity level is used (D3). Some significant differences, however, can be detected using the global magnetic activity level (D1) as an indicator.
Results with D1 are depicted in Fig. 8, where we compare the timeaveraged rotation profile to the different activity states (low, nominal, and high) denoted and defined as (5)where ⟨ Ω ⟩ _{t} is the timeaveraged rotation profile, computed over the whole simulation excluding only the initial kinematic stage when the magnetic field is still growing, and the quantities including angular brackets denote averages over a certain state. Using this definition, enhancement and quenching with respect to the average value during a certain state will show up as negative and positive values in the plot, respectively.
The timeaveraged rotation profile is very similar to the ones obtained and reported in various earlier works (Käpylä et al. 2012, 2013b), being solarlike, but with a somewhat more cylindrical profile, showing a midlatitude region of slower rotation leading to a region of reversed sign of radial and latitudinal differential rotation. This region was found to be instrumental in causing the equatorward migration of the magnetic field in Warnecke et al. (2014). The variations in the rotation profile averaged over different global magnetic activity states are weak (roughly 0.4 percent), while the instantaneous variations near the surface (see Fig. 7) could be as large as 5 percent during the early stages of the simulation (consistent with the slowly rotating convection of Karak et al. 2015). During a low/high state, slightly faster/slower equatorial rotation is obtained. This is consistent with the picture that the stronger the magnetic field, the stronger the suppression of the velocity field. The region with a reversed gradient of Ω, however, is observed to persist during all activity states, with virtually no change in its magnitude.
In Fig. 8 (middle row, radial differential rotation; lower row, latitudinal differential rotation) we investigate how the magnitude and distribution of differential rotation change depending on the global magnetic activity level. Here we define the radial shear as Δ_{r} = ∂Ω /∂r and the latitudinal shear as Δ_{θ} = ∂Ω /r∂θ, and the profiles computed over different activity states as (6)where the quantities with angular brackets are averages over the whole time series except for the kinematic stage (index t) or over a certain state (index i). It is evident that during high/low states of magnetic activity, the radial differential rotation is also quenched/enhanced in a similar way to the rotation itself, the magnitude of the change being roughly 1.5 percent in the equatorial region. No significant change is seen for the latitudinal differential rotation, except at the high latitude regions, which show hemispheric asymmetry especially during the low state. Finally, we define a dynamo number describing the magnitude of the radial differential rotation as (7)where the turbulent diffusivity is approximated by , and where is the local turbulent correlation time, α_{MLT} the mixing length parameter, and H_{p} = −(∂lnp(r,θ) /∂r)^{1} the pressure scale height. Following our previous studies, we use α_{MLT} = 5 / 3. The spatial distribution of C_{Ω} follows closely that of the radial gradient of the angular velocity (the leftmost panel in the middle row of Fig. 8), and therefore we do not plot it separately, but we merely report its spatial extrema [−30,110 ] averaged over the whole time span of the simulation, which is similar to those obtained in Käpylä et al. (2013b) and Warnecke et al. (2014).
Next we examine the meridional flow profiles and compute a dynamo number (8)where and other quantities are defined in the same way as in Eq. (7). The largest deviation from the timeaveraged profile is, similarly to the differential rotation, obtained when the activity states are defined based on the global magnetic field energy (D1), the magnitude of the change being maximally 6 percent in a thin layer concentrated in the equatorial region near the surface. The results are presented in Fig. 9, zoomed into the equatorial region. However, there are also some enhanced meridional flows generated near the latitudinal boundaries, but as these are most likely artifacts arising from the latitudinal boundaries, they are not shown here. The meridional flow pattern (Fig. 9, leftmost panel) is very similar to those obtained earlier by Käpylä et al. (2012, 2013b) and Warnecke et al. (2013, 2015), i.e., consisting of three cells in radius, aligned with the inner tangent cylinder. Additionally, there is one counterclockwise cell confined near the surface and equator, the nearsurface poleward flow being the strongest at this location. The strongest variations in magnitude over time are related to this region. The dynamo number varies spatially between 0 and 3 averaged over the whole simulation time span, the dynamo number of the meridional flow being more than thirty times weaker than of the Ω effect (radial differential rotation). The high (second panel in Fig. 9) and low (third panel in Fig. 9) states show no marked difference in the spatial distribution nor strength of the main cells; the small nearequator cell undergoes the strongest variations. The most interesting observation is the marked hemispheric asymmetry pronounced especially during the low state: the meridional flow in the southern hemisphere is more strongly quenched than the northern one. A similar effect, but more localized at the surface regions, is seen during the high state when the northern surface flow gets reduced, while the southern one gets stronger.
In conclusion, if the MM was an epoch of a global magnetic energy minimum, our results would predict faster surface rotation, stronger differential rotation, and hemispherically asymmetric meridional flow pattern. If a global magnetic energy maximum was actually occurring in the deeper layers of the convection zone during the MM, then our prediction would be reversed to slower surface rotation and weaker differential rotation, while the meridional flow would retain its asymmetric character. Neither of these scenarios is consistent with the actual observations of sunspot proper motion during the MM (Eddy et al. 1976; Ribes & NesmeRibes 1993). Overall, however, our analysis suggests that the surface activity disturbances have no direct and/or straightforward relation to changes in the rotational speed nor the strength of differential rotation. Moreover, the meridional circulation is so weak that it is most likely incapable of influencing the global dynamics significantly.
During the first 100 years of the simulation there are strong dips in the angular velocity, accompanied by simultaneous maxima in the global magnetic field strength. However, only one of these events (during t = 20−45 yrs) leads to the disappearance of the surface activity. Even then, the dip in the angular velocity does not coincide with the beginning and the duration of the disturbed surface activity period. Therefore, we can conclude that the surface irregularities do not originate from the variations seen in the angular velocity or meridional circulation, at the same time noting that a disturbed surface activity may not be indicative of a decrease in the global magnetic energy level.
4.5.2. Proxy of the α effect
We build a proxy for the isotropic α effect following the definition of Pouquet et al. (1976), (9)where the current helicity consists of the fluctuating current J′ = ∇ × B′/μ_{0} and the fluctuating magnetic field . When looking at the kinetic and magnetic contributions to the α effect proxy, it can be observed that the kinetic helicity is dominating the signal, the current helicity being roughly an order of magnitude weaker. Therefore, in the following we mostly concentrate on examining the properties of the kinetic helicity, while also computing a suitable dynamo number.
First we examine the time evolution of the azimuthally averaged kinetic helicity . This quantity shows a clear mean signal as a function of latitude, so that negative (positive) values are obtained in the north (south), increasing towards the poles where maximum values are obtained, with a strong localized enhancement seen in both hemispheres in the lowlatitude regions, extending to at least half of the depth of the convection zone. Close to the bottom of the convection zone at highlatitude regions, the sign of this quantity is reversed. These results are in agreement with earlier ones both in local and global domains (e.g., Käpylä et al. 2009, 2013b). In Fig. 10a, we show this quantity with the mean profile over the whole time series subtracted, showing that most of the variation occurs near the equatorial region. The signal is dominated by the highfrequency modes (M1 and M2), while at higher latitudes the dominant mode M7 in the magnetic field is detectable, but not significant in comparison to the highfrequency signal. The zoomedin plots of Figs. 10b and c reveal the presence of the weak basic cycle modulation in the kinetic helicity, but with half of the length of the magnetic cycle. Especially during times of clear antisymmetry (Fig. 10d) a modulation pattern is seen near the surface both at low latitudes (location of the strongest signal in kinetic helicity, solid lines in the figure), and at higher latitudes (location of the toroidal magnetic field maximum, dashed lines). The modulation pattern is such that the extrema of helicity roughly coincide with those of the radial magnetic field at midlatitudes (dotted lines). During times when M7 shows stronger symmetry with respect to the equator (Fig. 10e) all the observed weak dependencies in the earlier antisymmetric phase of the simulation break down, although signs of immediate restoration of the modulation pattern can be detected when the symmetry abruptly switches back (at around t = 234 yrs) to antisymmetric. Interestingly, signs of a longerterm modulation are seen in Fig. 10d, as during that epoch the kinetic helicity signal is systematically more negative (positive) in the north (south) in comparison to the mean value over the whole time span. This is also manifested in Fig. 10b as solid blue (yellow) stripes at low northern (southern) latitudes. Detected by EEMD only through the high value of the regularity measure, however, this modulation is of much lower energy than the highfrequency modes.
Fig. 10 a) Timelatitude plot of the kinetic helicity variations . Zoomedin plots of over t = 50−65 yrsb) and t = 225−240 yrsc). d) H_{kin} at ±10° (solid) and ± 25° (dashed) latitude for north (blue) and south (red) with the means over the whole time series indicated with an orange horizontal line. The data is averaged with a running mean over a 1year window, to smooth out the dominant highfrequency component. The dotted lines show the radial magnetic field at ± 25° latitude, scaled to fit the plot. All quantities are plotted at R_{s}. 
Fig. 11 From left to right: C_{α} timeaveraged over saturated stage. Difference to the C_{α} profile during the high, low, and nominal state of surface magnetic activity (D2). 
Next, we define and compute a dynamo number (10)This quantity is plotted in Fig. 11, left panel, as the time average over the whole length of the simulation. The profile is similar to those obtained earlier, e.g., in Runs C1 and C2 of Käpylä et al. (2013b).
This dynamo number is almost an order of magnitude larger than the corresponding value from meridional circulation, but roughly five times smaller than the value estimated from radial differential rotation. The strongest deviation from the mean profile is obtained when high and low states are determined according to the surface activity level (D2), while selecting the time points based on the global (D1) and bottom energy (D3) levels produce no significant deviation from the mean. During the low state, the α effect is strongly enhanced i.e., significantly higher positive and negative values are obtained in these regions in north and south, respectively. During the high state, there is also a mild enhancement, while in the nominal state the α effect is strongly reduced. The magnitude of the enhancement/quenching is the largest of all the dynamo drivers, roughly 30 percent, but acts in the opposite direction from that expected; i.e., with a locally and temporally larger C_{α} the naive expectation is to obtain a more efficient dynamo near the surface, but weaker surface activity is seen.
We also compute the migration direction of the magnetic field predicted by the ParkerYoshimura sign rule (Yoshimura 1976) (11)We plot this quantity together with the rms value of the mean toroidal magnetic field in the convection zone in Fig. 12, top panel. In the entire region of negative radial shear, manifested by the region of slower rotation in the differential rotation profile, equatorward migration is predicted to occur in both hemispheres. This region coincides with a toroidal field belt having a maximum roughly at ± 28° latitude and 0.85 R_{⊙} depth. There is another equatorward migration region at the very bottom of the convection zone, coinciding with the location of the strongest toroidal field belt roughly in a similar latitude range. In the upper third of the convection zone, however, the predicted migration direction is poleward. This matches with a toroidal magnetic field belt at ± 22° latitude and 0.9 R_{⊙} depth. No strong radial migration is predicted for the lower parts of the convection zone, where both upward and downward migration regions are seen, but without a clear systematic pattern. The good agreement of the predicted and actual migration direction as well the migration pattern itself fit well with the work of Warnecke et al. (2014, 2015). The migration pattern does not change significantly using any of the activity level definitions (D1−D3), which is likely because the gradients of the angular velocity and the kinetic helicity, and therefore the α effect, are sensitive to different activity indicators (the former to the global, the latter to the surface) and the signal is cancelled out. Therefore, we note that the α effect shows the most interesting behavior during the extrema seen in the surface activity, but our approach to investigating its effects via simple scalar proxies of kinetic and magnetic helicities is inadequate; we will return to this problem in a forthcoming paper.
Finally, to investigate whether the different cycle lengths of the dynamo modes and their distinct spatial distribution is due to a systematic, strong dependence of the magnetic diffusivity or magnetic Reynolds number (12)where η = 10^{8} m^{2} s^{1} is the molecular diffusivity profile either on radius or latitude, we plot their radial profiles at two different latitudes in Fig. 12, lower panel. To explain the coexistence of dynamo solutions with more than two orders of magnitude varying periods between the bottom and the top with the dependence being due to a weaker diffusivity in the bottom parts, a two orders of magnitude increase in the magnetic Reynolds number as a function of radius should be seen. Even though the trend seen is indeed increasing as a function of depth, the magnitude of the increase is at most four, and is too weak to provide an explanation.
Fig. 12 Top panel: color contours of in kGauss averaged over the whole simulation run, zoomedin to show the near equatorial region. The overplotted arrows show the predicted (Eq. (11)) timeaveraged mean migration direction ξ_{migr} normalized to unity. The dashed red lines show the radius at ± 25 degrees of latitude. Lower panel: radial dependence of the magnetic Reynolds number at two latitudes (north − solid; south − dashed). 
4.6. Equatorial symmetry
In Fig. 13 we plot the equatorial symmetry of the magnetic field, defined as (13)where E_{even} is the energy of the quadrupolar (symmetric) and E_{odd} the energy of the dipolar (antisymmetric) mode of the magnetic field. In Fig. 13a we show the global parity as a function of time. This quantity is obtained by computing, for each latitude pair for certain depth, the energies contained in the even and odd modes, and taking the mean of the obtained data. Evidently, the parity is mixed throughout the whole simulation, the dominant mode being the dipolar (solarlike) mode (average parity − 0.15 over the whole time series). A difference can again be observed between the first 3/4 of the simulation (average parity of − 0.26) in comparison to the last quarter, when there are larger fluctuations around the clearly positive mean parity of 0.1.
The best correlation between the different activity level definitions can be obtained when D3 is used. The parity not only shows a clear modulation with a similar longterm cycle to the toroidal field in the bottom of the convection zone, but also a clear pattern is revealed for the first 3/4 of the simulation with the D3 colorcoding: during a bottom mode minimum, the parity also obtains a minimum (− 0.42), while during a high state, the parity is the least solarlike (− 0.14). During the last quarter of the simulation, the bottom mode is so weak, that it falls below the low activity limit for the whole epoch. Nevertheless, the longterm cycle persists both in the bottom mode and in the parity, suggesting a relation between these quantities.
As can be seen from Table 2, the lowfrequency bottom mode (M11) is predominantly of dipolar symmetry (mean parity − 0.36). From Fig. 13b, it is evident that rather strong parity fluctuations are related to this mode, the excursions to positive parity being more short lived than the periods spent in the negative parity region. The majority of the parity variation, however, is related to the dominant basic cycle (M7), which on average has a parity of roughly − 0.12, but its symmetry strongly fluctuates as a function of time, as is evident from Fig. 13b. The parity fluctuates even more strongly than that of M11 between nearly purely dipolar and quadrupolar states. Most of the time M7 and M11 are in antiphase, i.e., when M11 attains the largest positive value, M7 is the most negative. There is a tendency for irregular behavior when the parities of these modes are in phase, such as during t = 20−45 yrs and t = 250−300 yrs. In Figs. 13c and d we show two zoomins of the parity evolution, panel (c) showing a clearly regular, antisymmetric behavior of the system and the surface toroidal field (t = 50−65 yrs), while panel (d) shows a predominantly symmetric state that during only a half a cycle changes into an antisymmetric one (t = 225−235 yrs). During the former epoch, the bottom toroidal magnetic field strength attains a minimum, while during the latter, a maximum is observed.
Fig. 13 a) Parity P as function of time. The colorcoding is defined using Definition D3. b) Parity of M7 (black) and M11 (red) as a function of time; their means are indicated with the horizontal lines with the same colorcoding. c) Parity (black) and toroidal magnetic field near the surface (blue: north, red: south) at ± 25° latitude, during the low state of the bottom toroidal mode, panel d) shows the same quantity, but for a high state of the bottom toroidal mode. 
Fig. 14 Top panel: magnetic cycle length determined from the toroidal magnetic field, , at ± 25° latitude and r = R_{s}, blue/red indicating positive/negative parity, north with triangle and south with plus symbols. Middle panel: time evolution of butterfly wing inclination. Bottom panel: butterfly wing widths, red/orange/blue indicating high/nominal/low bottom magnetic field energy (according to D3), north with triangle and south with plus symbols. 
4.7. Variation of general cycle properties
Evidently, the overall properties of the magnetic cycles show quite a large variation over time. In Fig. 14 we plot three indicators, namely the cycle length of the toroidal field near the surface at ± 25° latitude, the inclination of the butterfly wing, and the butterfly width in degrees for north (triangles) and south (plus signs) separately.
As for the toroidal field cycle length (Fig. 14, top panel), two distinct types of events can be seen. During suppressed magnetic activity near the surface layers (t = 20−45, 275−285, and 310−320 yrs), the field ceases to change its sign; therefore, some cycles are missed, and then cycles that are twice as long are detected as a result. The other type of variation are the less abrupt fluctuations around the mean cycle length, the shortest cycle lengths being slightly less than 4 years and the longest one around 8 years. The cycle length is not sensitive to any of the activity level definitions (D1−D3), but there is a tendency for the cycles to be shorter if the parity is more symmetric (blue symbols), while greater cycle length variations occur for more solarlike (antisymmetric; red symbols) parities. Also, all the “missing” cycles occur during the more solarlike parity.
To further analyze the general cycle properties we filtered the surface magnetic field data so that only points with were retained. This way the concentrations of strong magnetic field are clearly grouped into separate halfcycles or wings. The relative properties of these structures, such as the latitudinal extent (width) and duration can then be directly measured. Based on these properties we define the inclination of the butterfly wing as the half cycle width divided by its duration. For north/south, a solarlike cycle would have negative/positive inclinations, while for cycles where migration is weak or nonexistent, values close to zero are expected. As is evident from the middle panel of Fig. 14, there is a tendency for the missing cycles to appear with weak migration only. Antisolar cycles are also observed a number of times both for north and south, but not at simultaneous epochs. It appears that the shorter the cycle, the higher the probability for antisolar migration pattern.
The corresponding butterfly width (Fig. 14 bottom panel) is most often large for the missing cycles, which is in striking contrast to the observational evidence for the MM (see Usoskin et al. 2015, and references therein), if these two types of events in reality represent one and the same phenomenon. There is a tendency for narrow widths to occur when the bottom magnetic field is weak.
Fig. 15 From top to bottom: zoomins of the timelatitude diagram of the mean toroidal magnetic field near the surface (at R_{s}) and bottom (at R_{b}), the angular velocity variations, and the Reynolds stress components R_{rφ} and R_{θφ} near the surface (at R_{s}) during t = 20−65 yrs, showing that the magnetic cycle and surface activity are strongly disturbed (t = 20−45 yrs), after which the regular state is quickly recovered (t = 45−65 yrs). 
4.8. Irregular event during t = 20−45 yrs
Finally, we present a closeup of the most clearly pronounced irregular activity epoch during t = 20−45 yrs. In Fig. 15, we reproduce the plot shown in Fig. 7, but restrict the plotting range to 20−65 years, which shows both the irregular evolution, followed by the resumption of the regular behavior, indicating that the system changes itself very abruptly after the distortion. The time resolution in the figure is high enough to clearly show the highfrequency (M1) mode of the magnetic field (top), not clearly discernible in any other longer time span figure; this cyclicity is seen to persist regardless of the irregularities seen in the modes M7 and M11.
The irregular behavior is first seen in the southern hemisphere, see Fig. 15a, and the event appears to be preceded by a field enhancement in the bottom mode (M11; Fig. 15b). Even though the toroidal field of the mode M7 gets very weak, the surface activity does not totally cease, but the field evolves with significantly longer period without a polarity change. The southern cycle, therefore, appears truly as a missed one; meanwhile, the northern hemisphere continues as usual. After twice the time of a normal cycle has elapsed, the bottom magnetic field attains another strong maximum. This is reflected by a strong minimum in the Reynolds stresses, Figs. 15d and e, and surface differential rotation, Fig. 15c. Immediately after this, the toroidal magnetic field practically vanishes from the northern hemisphere, while the south seems totally unaffected by the northern disruption. Again, after the south has produced two normal magnetic cycles, normal activity is also rapidly restored in the north. After that, the bottom magnetic field exhibits no further strong extrema, and the mode M7 continues its evolution in a regular manner. It is also noteworthy that the bottom magnetic field changes its polarity during the disturbed epoch.
During this epoch, it is evident even by eye that there is especially strong equatorial (northsouth) asymmetry related to the surface toroidal magnetic field. This is reminiscent of the observations of the Maunder and Dalton minima (Ribes & NesmeRibes 1993; Usoskin et al. 2009), which indicate strong hemispheric asymmetry and a relatively strong quadrupolar component of the magnetic field. The global parity does not strongly reflect this local surface anomaly as a sudden, persistent change towards a quadrupolar state (+1), although the quantity fluctuates strongly between the values [− 0.9, 0.7]. This is most likely due to the fact that the symmetry properties of the very strong bottom toroidal magnetic field dominate the signal, at all times remaining close to a dipolar configuration.
Although similarly dramatic events cannot be isolated from the rest of the time series, it is clear that both extrema and polarity reversals of the bottom mode M11, frequent during the first 3/4 of the simulation time, make the system more irregular than during the last 1/4 when especially the energy contained in the bottom mode M11 gets weaker.
5. Conclusions
In this paper we have analyzed a semiglobal (wedgeshaped) DNS that produces a solarlike oscillatory dynamo with surface migration properties of the magnetic field resembling those from observations. The mean cycle length is 4.9 years, and the simulation is evolved over 80 such magnetic cycles. If scaled to solar units using the fact that the Sun has a magnetic cycle of roughly 22 years, our simulation would correspond to two millennia of solar evolution. The chosen parameters, most notably the modest stratification and the values of the Reynolds numbers, are still far from the real Sun; nevertheless, this combination of parameters produces a solarlike dynamo, and it is computationally affordable to integrate it over a large number of cycles. The rotation rate is five times solar, resulting in a differential rotation that is roughly three times weaker than in the Sun. The rotation profile is solarlike, but somewhat more cylindrical, and the meridional circulation consists of multiple cells in radius. Nevertheless, we expect that this simulation can easily be used to study the longterm evolution of solartype dynamos.
A general property of the dynamo solution is its cyclic nature: even though there is a clear magnetic cycle, the various nonlinearities in the system drive it away from an exactly periodic (harmonic) state. Therefore, special time series analysis techniques are needed that can cope with nonperiodic signals. In this study, we perform a statistical analysis over the whole convection zone, investigating the properties of all relevant quantities at different depths and latitudes. The methods we choose to employ are the EEMD and D^{2} statistics.
The behavior of the dynamo solution is extremely complex. In addition to changing cycle length, we observe epochs of disturbed and even ceased surface activity, and strong shortterm hemispherical asymmetries. The major general findings related to the overall dynamics of the system include the following:

The hemispheric asymmetries are related to the magnetic field alone, while the velocity field remains almost perfectly symmetric at all times.

The epochs when the surface activity has decreased or is practically nonexistent are not global magnetic energy minima but maxima, as strong magnetic fields are stored in the deeper parts of the convection zone.
The main goal of this study is to find the causes for the irregular behavior rather than to make direct comparisons with the observed characteristics of the solar magnetic cycle. Therefore, we have concentrated our efforts toward analyzing the dynamo solution itself to investigate how the most important properties (cycle length, migration properties, energetics) change during the irregular epochs. We also investigated all the key factors in the dynamo process, namely the rotation and its nonuniformities, meridional circulation, the inductive action arising from turbulent convection (α effect), and the changes of these as functions of a relevant activity level measures. The major findings from this part of the analysis can be summarized as follows:

The specific dynamo solution analyzed here contains not only one, but three separate modes, having distinct cycle lengths and symmetry properties, and are located in different parts of the convection zone. The dominant mode is the nearsurface 4.9 year cycle (denoted with M7) showing equatorward migration at lower latitudes and poleward migration at higher latitudes. This mode is accompanied by a weaker highfrequency poleward migrating mode (M1, 0.11 years) in the equatorial region and a lowfrequency mode at the bottom of the convection zone (M11, roughly 50 years).

The crucial property of the different dynamo modes is their different symmetries. While the dominating surface mode has a mixed equatorial symmetry that undergoes strong fluctuations over time, the bottom mode exhibits nearly pure antisymmetry and is also more regular than the surface mode.

There is a close relationship between the global magnetic and kinetic energies such that strong magnetic fields quench the flow field in a manner that the angular velocity is significantly reduced, differential rotation gets somewhat weaker, and the meridional circulation attains an asymmetric character with modifications of the circulation magnitude especially in nearequatorial surface regions. We expect that this is what would have been seen during the MM, if it was an event corresponding to the abrupt disappearance of surface activity with a simultaneous global energy maximum in the bottom of the convection zone (see, e.g., Küker et al. 1999; Karak 2010, for similar arguments).

A more common way of interpreting MM is an overall drop in the magnetic activity level. Our simulation produces no such events, which is naturally not proof that they do not occur in the real Sun. The only quantity markedly linked to the surface magnetic activity level in our model is the magnitude of the α effect.

Two kinds of irregularities are separable from the dynamo solution itself: smooth variations in the cycle properties and abrupt changes that lead to missing cycles and ceased surface activity. All these can be satisfactorily explained as the interplay of the different dynamo modes and their influence on the flow field.
Even though this simulation was run for a considerable length of time, we still observe a secular trend − especially for the longterm cycle in the bottom of the convection zone. To fully capture and understand this trend, and to find out whether it is a transient or a real cycle, the simulation would still need to be continued further. Although we present a rather elaborate and time consuming data analysis of the results, we have merely touched upon the turbulent quantities. In particular, the α effect was treated in a very simplified manner, only using a proxy from the kinetic helicity, which is far from adequate. We recently undertook an effort to compute the full tensorial representation of both the diffusive and antidiffusive contributions to the electromotive force using the spherical testfield method (Warnecke et al. 2016) from a similar, but shorter, solarlike dynamo simulation. It would therefore be useful to compute and compare the turbulent transport coefficients from the different epochs of the long run presented here. The simplified analysis presented here, however, reaffirms our earlier conclusion that the general dynamo solution and its migration can be explained in terms of a turbulent αΩ dynamo (Käpylä et al. 2012; Warnecke et al. 2014).
Acknowledgments
We thank the anonymous referee for useful comments that helped us to improve the contents and presentation of the manuscript. The simulations were performed using the supercomputers hosted by the CSC − IT Center for Science Ltd. in Espoo, Finland, which is administered by the Finnish Ministry of Education and in the HLRS supercomputing centre in Stuttgart, Germany, through the PRACE allocation “SOLDYN”. Financial support from the Academy of Finland ReSoLVE Center of Excellence (grant No. 272157; MJK, NO) and grants No. 136189, 140970, 272786 (PJK), the Swedish Research Council grants 62120115076 and 20125797 (AB), and from the Estonian Research Council (Grant IUT401; JP) are acknowledged, as are the MaxPlanck/Princeton Center for Plasma Physics and funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/20072013) under REA grant agreement No. 623609 (JW). BBK is supported by the NASA Living With a Star Jack Eddy Postdoctoral Fellowship Program, administered by the University Corporation for Atmospheric Research.
References
 Augustson, K., Brun, A. S., Miesch, M., & Toomre, J. 2015, ApJ, 809, 149 [Google Scholar]
 Barnhart, B. L. 2011, Ph.D. Thesis, University of Iowa [Google Scholar]
 Brandenburg, A., Chan, K. L., Nordlund, Å., & Stein, R. F. 2005, Astron. Nachr., 326, 681 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P. 2010, Liv. Rev. Sol. Phys., 7, 3 [Google Scholar]
 Charbonneau, P. 2014, ARA&A, 52, 251 [Google Scholar]
 Choudhuri, A. R., & Karak, B. B. 2012, Phys. Rev. Lett., 109, 171103 [NASA ADS] [CrossRef] [Google Scholar]
 Choudhuri, A. R., Schüssler, M., & Dikpati, M. 1995, A&A, 303, L29 [NASA ADS] [Google Scholar]
 Cohen, L. 1995, Timefrequency Analysis, Electrical engineering signal processing (Prentice Hall PTR) [Google Scholar]
 Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508 [NASA ADS] [CrossRef] [Google Scholar]
 Eddy, J. A., Gilman, P. A., & Trotter, D. E. 1976, Sol. Phys., 46, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Fan, Y., & Fang, F. 2014, ApJ, 789, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Flandrin, P., Rilling, G., & Goncalves, P. 2004, Signal Processing Letters, IEEE, 11, 112 [Google Scholar]
 Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Gilman, P. A. 1983, ApJS, 53, 243 [NASA ADS] [CrossRef] [Google Scholar]
 Glatzmaier, G. A. 1985, ApJ, 291, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, N. E., Shen, Z., Long, S. R., et al. 1998, Proc. Roy. Soc. London A, 454, 903 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, N. E., & Wu, Z. 2008, Rev. Geophys., 46 [Google Scholar]
 Ivanov, V. G., & Miletsky, E. V. 2011, Sol. Phys., 268, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2004, A&A, 422, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Tuominen, I. 2006, Astron. Nachr., 327, 884 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Korpi, M. J., & Brandenburg, A. 2009, A&A, 500, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2011, Astron. Nachr., 332, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, ApJ, 755, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2013a, Geophys. Astrophys. Fluid Dynam., 107, 244 [Google Scholar]
 Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013b, ApJ, 778, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karak, B. B. 2010, ApJ, 724, 1021 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B. B., Jiang, J., Miesch, M. S., Charbonneau, P., & Choudhuri, A. R. 2014, Space Sci. Rev., 186, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B. B., Käpylä, P. J., Käpylä, M. J., et al. 2015, A&A, 576, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krause, F., & Rädler, K.H. 1980, Meanfield Magnetohydrodynamics and Dynamo Theory (Oxford: Pergamon Press) [Google Scholar]
 Küker, M., Arlt, R., & Rüdiger, G. 1999, A&A, 343, 977 [NASA ADS] [Google Scholar]
 Lindborg, M., Mantere, M. J., Olspert, N., et al. 2013, A&A, 559, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mabuchi, J., Masada, Y., & Kageyama, A. 2015, ApJ, 806, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Malkus, W. V. R., & Proctor, M. R. E. 1975, J. Fluid Mech., 67, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Masada, Y., & Sano, T. 2014, ApJ, 794, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Miesch, M. S., Featherstone, N. A., Rempel, M., & Trampedach, R. 2012, ApJ, 757, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Moffatt, H. K. 1978, Magnetic Field Generation in Electrically Conducting Fluids (Cambridge: Cambridge University Press) [Google Scholar]
 Moss, D., Sokoloff, D., Usoskin, I., & Tutubalin, V. 2008, Sol. Phys., 250, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Norton, A. A., Charbonneau, P., & Passos, D. 2014, Space Sci. Rev., 186, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Olspert, N., Käpylä, M. J., Pelt, J., et al. 2015, A&A, 577, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ossendrijver, M. 2003, A&ARv, 11, 287 [Google Scholar]
 Ossendrijver, A. J. H., Hoyng, P., & Schmitt, D. 1996, A&A, 313, 938 [NASA ADS] [Google Scholar]
 Parker, E. N. 1955, ApJ, 122, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Passos, D., & Charbonneau, P. 2014, A&A, 568, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pelt, J. 1983, in Statistical Methods in Astronomy, ed. E. J. Rolfe, ESA SP, 201, 37 [Google Scholar]
 Pipin, V. V., & Seehafer, N. 2009, A&A, 493, 819 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pouquet, A., Frisch, U., & Léorat, J. 1976, J. Fluid Mech., 77, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Qian, S. 2002, Introduction to Timefrequency and Wavelet Transforms (Prentice Hall PTR) [Google Scholar]
 Racine, É., Charbonneau, P., Ghizaru, M., Bouchat, A., & Smolarkiewicz, P. K. 2011, ApJ, 735, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Ribes, J. C., & NesmeRibes, E. 1993, A&A, 276, 549 [NASA ADS] [Google Scholar]
 Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solartype Stars (Berlin: Akademie Verlag) [Google Scholar]
 Schrinner, M., Rädler, K.H., Schmitt, D., Rheinhardt, M., & Christensen, U. 2005, Astron. Nachr., 326, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Schrinner, M., Rädler, K.H., Schmitt, D., Rheinhardt, M., & Christensen, U. R. 2007, Geophys. Astrophys. Fluid Dyn., 101, 81 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Schrinner, M., Petitdemange, L., & Dormy, E. 2011, A&A, 530, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schrinner, M., Petitdemange, L., & Dormy, E. 2012, ApJ, 752, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Simitev, R. D., Kosovichev, A. G., & Busse, F. H. 2015, ApJ, 810, 80 [Google Scholar]
 Sokoloff, D., & NesmeRibes, E. 1994, A&A, 288, 293 [NASA ADS] [Google Scholar]
 Usoskin, I. G., Solanki, S. K., & Kovaltsov, G. A. 2007, A&A, 471, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin, I. G., Mursula, K., Arlt, R., & Kovaltsov, G. A. 2009, ApJ, 700, L154 [Google Scholar]
 Usoskin, I. G., Arlt, R., Asvestari, E., et al. 2015, A&A, 581, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vitense, E. 1953, Z. Astrophys., 32, 135 [NASA ADS] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2012, Sol. Phys., 280, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Mantere, M. J., & Brandenburg, A. 2013, ApJ, 778, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, ApJ, 796, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke, J., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2015, A&A, submitted [arXiv:1503.05251] [Google Scholar]
 Warnecke, J., Rheinhardt, M., Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2016, A&A, submitted [arXiv:1601.03730] [Google Scholar]
 Wu, Z., & Huang, N. E. 2009, Advances in Adaptive Data Analysis, 01, 1 [Google Scholar]
 Yadav, R. K., Gastine, T., Christensen, U. R., & Duarte, L. D. V. 2013, ApJ, 774, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Yoshimura, H. 1976, Sol. Phys., 50, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Zolotova, N. V., & Ponyavin, D. I. 2015, ApJ, 800, 42 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Empirical mode decomposition
Fourier spectral analysis provides a general method for examining the global energy−frequency distributions. For its validity, some crucial conditions must hold: the system must be linear and the data must be strictly periodic or stationary. If these restrictions are not satisfied, the nonstationary nature of the data causes the energy to be spread over a wide frequency range. Additionally, deformed waveprofiles, which are a direct consequence of nonlinear effects, need additional harmonic components to be fitted. As a result, the energy−frequency distribution can be misleading and hard to interpret (Huang et al. 1998).
Different analysis methods have been developed for the purpose of describing nonstationary signals. Wellknown examples are shorttime Fourier transform and wavelet transform (e.g., Qian 2002; Cohen 1995). An alternative to the aforementioned methods is the HilbertHuang transform (HHT; Huang et al. 1998). In contrast to almost all the other methods, HHT works directly in the temporal domain rather than in the corresponding frequency space, and the basis functions, also known as intrinsic mode functions (IMFs), are derived from the data not selected a priori. The decomposition implicitly makes the simple assumption that, at any given time, the data may have many coexisting simple oscillatory modes of significantly different frequencies, one superimposed on the other (Huang & Wu 2008).
The HilbertHuang transform is performed in two steps: first, using an algorithm called Empirical Mode Decomposition (EMD), by which the signal is separated into a set of IMFs; second, for each extracted mode, Hilbert spectral analysis is applied which allows each mode to be described as an analytical signal having the form A(t)exp(iϕ(t)), where A(t) is an instantaneous amplitude and ϕ(t) an instantaneous phase. For signals that result in such a form, the low frequency content is in the amplitude term and the high frequency content in the exponential term (Cohen 1995). Having obtained a decomposition into IMFs that satisfy the analytic signal conditions, we can localize any event on the time and on the frequency axis. If local timedependent aspects of the IMFs are not of interest the second step in the above algorithm can be replaced by a crude approach that estimates the average mode period and amplitude by counting zerocrossings and determining points of extrema.
One of the major drawbacks of the original EMD was the problem known as mode mixing, which is a consequence of signal intermittency. To overcome this problem, a noiseassisted data analysis method called Ensemble EMD (EEMD) was proposed (Wu & Huang 2009). In Flandrin et al. (2004) it was shown that when applied to Gaussian white noise, EMD acts as a dyadic filter bank. Utilizing this property allows robust and statistically significant IMFs to be extracted.
The Ensemble EMD method was previously applied to time series of different solar activity proxies. An example of the application to total solar irradiance and sunspot data can be found in Barnhart (2011).
As an illustration of how EEMD works we selected a time series of taken at R_{b} and a latitude of 22°. This time series was decomposed into 12 IMFs, most of the energy being distributed over modes 7−11, each of which individually contribute more than 10% to the full energy. The given modes, with the corresponding instantaneous amplitudes are shown in the top 5 panels in Fig. A.1. In the bottom panel the original time series of and the sum of these five modes are shown. The small difference between the curves comes from the fact that the less significant modes (1−6, 12 and 13) have been excluded from the sum. We also note that the form of the IMFs does not precisely satisfy the criteria of an analytic signal on the one hand due to the finite precision limit in the IMF extraction algorithm and on the other hand due to the relatively small ensemble size (the noise does not fully cancel out).
Fig. A.1 Modes 7−11 of at R_{b} and a latitude of 22°. 
Appendix B: D^{2} phase dispersion statistic
The D^{2} phase dispersion statistic was first introduced in Pelt (1983). Recent applications of the method can be found in Lindborg et al. (2013), Karak et al. (2015), and Olspert et al. (2015). Similarly to EEMD, the D^{2} statistic is well suited for nonstationary data. However, while EEMD aims at decomposing the initial time series into a set of quasiharmonic signals, the purpose of D^{2} is to find a set of average cycle lengths, P_{m}, that are consistent with the full data set. The D^{2} statistic is defined as (B.1)where f(t_{i}),i = 1,...,N is the input time series, σ^{2} is its variance, g(t_{i},t_{j},P,Δt) is the selection function, which is significantly greater than zero only when where P is the trial period and Δt is the socalled coherence time, which is the measure of the width of the sliding time window wherein the data points are taken into account by the statistic. More precisely, in the given study the function g was chosen as the product of cosine and Gaussian functions. Here the former depends only on the phase difference of given points with respect to the given trial period and the latter only on the time difference between the corresponding points:
where frac((t_{j} − t_{i}) /P) removes the integer part of (t_{j} − t_{i}) /P.
As Δt is made shorter, we match nearby cycles in a progressively narrower region, and consequently estimate a certain mean cycle length P_{m}, not necessarily coherent for the full time span. At large values of Δt, i.e., when the time window width approaches the time span of the data, the D^{2} spectrum approaches the result that would be obtained using Fourier transform, for example. More precisely, P_{m} is determined from the dispersion spectrum calculated for the longest coherence time for which the shape of the minimum is still symmetric and singular (with a single local minimum within a reasonably narrow search region to exclude possible minima from higher harmonics). Denoting this optimal coherence time by Δt_{opt} we define the mean cycle length as (B.7)When visualizing the spectrum, we plot the cycle length P on vertical axis, but instead of directly plotting the coherence time on horizontal axis, we use Δt/P instead. This is just a count of cycles with a corresponding period that fit into the time interval equal to given coherence time. We call this entity coherence length.
All Tables
Input and diagnostic parameters together with the kinetic and magnetic energies realized in the simulation in units of 10^{5} J m^{3}.
All Figures
Fig. 1 Time evolution of the mean toroidal magnetic field in the convection zone for three depths (from top to bottom, R_{s} = 0.98 R_{⊙}, R_{m} = 0.85 R_{⊙}, and R_{b} = 0.72 R_{⊙}). The extrema are [−11.0,10.3 ] kG at R_{s}, [−23.3,20.7 ] kG at R_{m}, and [−29.9,31.7 ] kG at R_{b}. 

In the text 
Fig. 2 Evolution of the mean radial magnetic field in the convection zone for three depths (from top to bottom, R_{s}, R_{m}, and R_{b}). The color ranges have been clipped at half of the maxima/minima, i.e., [−11.7,10.4 ] kG at R_{s}, [−6.1,5.7 ] kG at R_{m}, and [−3.3,3.8 ] kG at R_{b}, to emphasize the weaker structures at low latitudes. 

In the text 
Fig. 3 a) Total (, purple dashed), total mean (, black solid), toroidal (, blue), and poloidal (, red) rms values as functions of time as averages over the whole volume. b) in the north (blue) and south (red). The black solid line shows the ratio of the northern to southern energy as a function of time, and the orange line the mean of the ratio computed over the whole time span. c) The same for the poloidal field. d) The bottom (at R_{b}, blue line) and surface (at R_{s}, red line) toroidal magnetic field strengths as functions of time. The thin orange horizontal lines indicate the mean of the quantities over the whole simulation time. In a) and d) thick orange lines show the smoothed curve of the total mean magnetic field a) and the bottom/surface toroidal magnetic field strength d); the yellow horizontal lines show 1.1 and 0.9 times the mean levels, respectively. 

In the text 
Fig. 4 a) Rms value of the mean azimuthal velocity as a function of time together with the average over the whole time series plotted with a black horizontal line. The dashed black line shows the volumeaveraged rms value of the fluctuating velocity field, u_{rms}. The times considered “high” are indicated in red, “low” in blue, and “nominal” activity states according to the strength of the total magnetic field in orange (D1, see Sect. 4.4). b) As in a), but for the rms value of the total magnetic field as a function of time. The black horizontal line indicates the mean over the whole time series. c) Rms value of the meridional velocity as a function of time (black solid line) with the mean over the whole time series (orange horizontal line) and a 4.9year moving average (thick red solid line). d) Zoomin to 50−65 years of evolution of the mean azimuthal velocity (black solid), mean meridional velocity (blue line, multiplied by a factor of 10 to fit the figure), and the total mean magnetic field strength (red line, scaled to fit the figure). 

In the text 
Fig. 5 Distribution of mean amplitudes of the most significant modes of mean magnetic and mean velocity fields found from EEMD. The figures are labeled with the physical variable followed by mode index (e.g. 7 indicates mode 7 of the mean latitudinal magnetic field). Colors reflect the mean amplitude (blue − high, yellow − low) of the mode at the given location; we note that the scales of separate figures are different. Contours on the plot represent the lines of constant amplitude with values given in the legend. 

In the text 
Fig. 6 Phase dispersion analysis results for the components of . The calculations were done for at latitude ± 22° and radius 0.94 R_{⊙}, for at latitude ± 66° and radius 0.94 R_{⊙}, and for at latitude ± 49° and radius 0.82 R_{⊙}. Left (right) refers to the northern (southern) hemisphere. 

In the text 
Fig. 7 From top to bottom: timelatitude diagram of the angular velocity variations Ω^{v}, variations of the largescale Lorentz force , and the Reynolds stress components R_{rφ} and R_{θφ} near the surface (at R_{s}). For the Reynolds stress plots, the contours have been clipped at half of the maxima/minima. 

In the text 
Fig. 8 Top row, from left to right: normalized timeaveraged rotation profile ⟨ Ω ⟩ _{t}/ Ω_{0}, and difference to the rotation profile ΔΩ_{state} during the high, low, and nominal state of global activity (D1). Middle (bottom) row, the same for mean radial (latitudinal) differential rotation. 

In the text 
Fig. 9 Meridional flow, indicated as arrows, and timeaveraged dynamo parameter C_{u}, indicated with color contours. On the left, we show the meridional circulation profile averaged over the whole simulation time span, and the next three panels show the difference ΔC_{u, ⟨ state ⟩} = C_{u} − C_{u, ⟨ state ⟩} to the mean during a high, low, and nominal global activitystate (D1). 

In the text 
Fig. 10 a) Timelatitude plot of the kinetic helicity variations . Zoomedin plots of over t = 50−65 yrsb) and t = 225−240 yrsc). d) H_{kin} at ±10° (solid) and ± 25° (dashed) latitude for north (blue) and south (red) with the means over the whole time series indicated with an orange horizontal line. The data is averaged with a running mean over a 1year window, to smooth out the dominant highfrequency component. The dotted lines show the radial magnetic field at ± 25° latitude, scaled to fit the plot. All quantities are plotted at R_{s}. 

In the text 
Fig. 11 From left to right: C_{α} timeaveraged over saturated stage. Difference to the C_{α} profile during the high, low, and nominal state of surface magnetic activity (D2). 

In the text 
Fig. 12 Top panel: color contours of in kGauss averaged over the whole simulation run, zoomedin to show the near equatorial region. The overplotted arrows show the predicted (Eq. (11)) timeaveraged mean migration direction ξ_{migr} normalized to unity. The dashed red lines show the radius at ± 25 degrees of latitude. Lower panel: radial dependence of the magnetic Reynolds number at two latitudes (north − solid; south − dashed). 

In the text 
Fig. 13 a) Parity P as function of time. The colorcoding is defined using Definition D3. b) Parity of M7 (black) and M11 (red) as a function of time; their means are indicated with the horizontal lines with the same colorcoding. c) Parity (black) and toroidal magnetic field near the surface (blue: north, red: south) at ± 25° latitude, during the low state of the bottom toroidal mode, panel d) shows the same quantity, but for a high state of the bottom toroidal mode. 

In the text 
Fig. 14 Top panel: magnetic cycle length determined from the toroidal magnetic field, , at ± 25° latitude and r = R_{s}, blue/red indicating positive/negative parity, north with triangle and south with plus symbols. Middle panel: time evolution of butterfly wing inclination. Bottom panel: butterfly wing widths, red/orange/blue indicating high/nominal/low bottom magnetic field energy (according to D3), north with triangle and south with plus symbols. 

In the text 
Fig. 15 From top to bottom: zoomins of the timelatitude diagram of the mean toroidal magnetic field near the surface (at R_{s}) and bottom (at R_{b}), the angular velocity variations, and the Reynolds stress components R_{rφ} and R_{θφ} near the surface (at R_{s}) during t = 20−65 yrs, showing that the magnetic cycle and surface activity are strongly disturbed (t = 20−45 yrs), after which the regular state is quickly recovered (t = 45−65 yrs). 

In the text 
Fig. A.1 Modes 7−11 of at R_{b} and a latitude of 22°. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.