Issue 
A&A
Volume 568, August 2014



Article Number  A113  
Number of page(s)  16  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201423700  
Published online  02 September 2014 
Characteristics of magnetic solarlike cycles in a 3D MHD simulation of solar convection^{⋆}
^{1} CENTRA, Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049001 Lisboa, Portugal
email: dariopassos@ist.utl.pt
^{2} GRPS, Départment de Physique, Université de Montréal, CP 6128, Centreville, Montréal, Qc, H3C3J7 Canada
^{3} Departamento de Física, Universidade do Algarve, Campus de Gambelas, 8005139 Faro, Portugal
Received: 23 February 2014
Accepted: 30 May 2014
We analyse the statistical properties of the stable magnetic cycle unfolding in an extended 3D magnetohydrodynamic simulation of solar convection produced with the EULAGMHD code. The millennium simulation spans over 1650 years, in the course of which forty polarity reversals take place on a regular ~40 yr cadence, remaining wellsynchronized across solar hemispheres. In order to characterize this cycle and facilitate its comparison with measures typically used to represent solar activity, we build two proxies for the magnetic field in the simulation mimicking the solar toroidal field and the polar radial field. Several quantities that characterize the cycle are measured (period, amplitudes, etc.) and correlations between them are computed. These are then compared with their observational analogs. From the typical GnevyshevOhl pattern, to hints of Gleissberg modulation, the simulated cycles share many of the characteristics of their observational analogs even though the simulation lacks poloidal field regeneration through active region decay, a mechanism nowadays often considered an essential component of the solar dynamo. Some significant discrepancies are also identified, most notably the inphase variation of the simulated poloidal and toroidal largescale magnetic components, and the low degree of hemispheric coupling at the level of hemispheric cycle amplitudes. Possible causes underlying these discrepancies are discussed.
Key words: dynamo / Sun: interior / Sun: activity / Sun: magnetic fields / magnetohydrodynamics (MHD) / sunspots
Appendix is available in electronic form at http://www.aanda.org
© ESO, 2014
1. Introduction
The recent, extended phase of very low activity having delineated sunspot cycle 23 from the current cycle 24 brought home once again the fact that we do not understand what sets the duration and amplitude of successive solar magnetic activity cycles. Sunspot cycles, welldocumented for four centuries, are characterized by a mean period of about eleven years, half the period of the underlying magnetic cycle; sunspots form in the same way irrespective of the polarity of the solar internal magnetic field. However, successive cycles do not unfold as a steady harmonic oscillation, but show instead considerable variations in amplitude and duration (see Hathaway 2010, and references therein). These fluctuations are not a consequence of the vagaries of sunspot formation, as they appear equally clearly in other quantitative indicators of solar activity with long temporal records, such as the F10.7 radio flux. At their most extreme, cycle fluctuations include the socalled grand minima, extended periods of strongly suppressed activity spanning many sunspot cycles. The best known is the 1645–1715 Maunder Minimum (Eddy 1976), but indirect indicators of solar activity have revealed many more such grand minima in the pretelescopic era, irregularly interspersed within epochs of normal activity (Usoskin 2013).
It is now generally agreed that the solar magnetic cycle is powered by a dynamo mechanism associated with the inductive action of fluid flows within the solar interior. Considering that the outer 30% of the sun’s radius is in a state of strongly turbulent thermallydriven convection, it is perhaps not surprising to observe cycletocycle fluctuations, simply as a result of stochastic perturbation of the dynamo. Moreover, the nonlinear dynamical backreaction of the induced magnetic field on these same flows can also generate deterministic modulations of the cycle on a variety of timescales. These stochastic and deterministic effects have been modelled extensively in the context of (relatively) simple meanfield and meanfieldlike dynamo models (see Charbonneau 2010, Sect. 5 for some illustrative examples). At least at a qualitative level, many observed solar cycle fluctuation patterns can be reproduced by proper adjustment of model parameters and free functionals. Yet, at this point in time, comparing such modelling results with observations still does not allow us to unambiguously establish the relative importance of stochasticity from nonlinear deterministic effects (Weiss 2010).
In principle, global magnetohydrodynamic (MHD) simulations of solar convection embody selfconsistently both stochastic forcing and nonlinear magnetic backreaction on all spatial and temporal scales resolved by the simulations. In the past few years, many such simulations, although distinct in their computational design, have achieved the production of axisymmetric largescale magnetic fields undergoing more or less regular polarity reversals, some even producing sustained solarlike cycles (Augustson et al. 2014; Brown et al. 2010, 2011; Charbonneau & Smolarkiewicz 2013; Ghizaru et al. 2010; Käpylä et al. 2010, 2012, 2013; Racine et al. 2011). Although they still lack some potentially important field regeneration mechanisms – most notably perhaps the surface decay of active regions and associated buildup of a surface dipole moment, – they represent unique virtual laboratory experiments allowing us to investigate the dynamics of convectivelydriven magnetic cycles.
This paper presents an observational analysis of a very long MHD simulation of solar convection, akin to those presented in Ghizaru et al. (2010) and Racine et al. (2011), spanning 20 full magnetic cycles, equivalent to 40 sunspot cycles. An overview of the simulation’s design and overall characteristics is first presented in Sect. 2, followed in Sect. 3 by the definition of the various simulationbased proxies built to allow comparison with sunspot data. In Sect. 4 we examine the statistical properties of our pseudosunspot cycles, and the degree to which they exhibit known patterns in sunspot cycle data, namely the socalled Waldmeier and GnevyshevOhl Rules. In Sect. 5 we repeat some of these analyses treating each solar hemisphere separately, which allows us to ascertain the degree of hemispheric coupling characterizing this specific simulation. We present in Sect. 6 a speculative discussion of the possible physical origin of the patterns extracted from the simulation, and close in Sect. 7 by summarizing our results, and their implications for the prediction of solar cycle characteristics.
2. The simulation
The simulation used in the foregoing analysis solves the anelastic MHD equations in a thick, stratified spherical shell of electrically conducting fluid, spanning from r/R_{⊙} = 0.63 to 0.96 in fractional solar radius and rotating at the solar rate. The stratification in the outer two thirds of the domain (0.713 ≤ r/R_{⊙} ≤ 0.96) is weakly convectively unstable, and the portion below is stably stratified. Thermal forcing within the convectively unstable layers sustains highly supercritical, turbulent convective fluid motions. The overall simulation design is described in Ghizaru et al. (2010), Racine et al. (2011), and Guerrero et al. (2013). The governing equations are solved using the EULAGMHD code, the MHD generalization described in Smolarkiewicz & Charbonneau (2013) of the robust multiscale flow solver EULAG (Prussa et al. 2008).
EULAGMHD belongs to the class of socalled Implicit LargeEddy Simulations (ILES), in which nonlinear stability is maintained by implicit dissipation introduced at the level of the numerical advection scheme, rather than through explicit enhanced physical dissipation and/or subgrid model. EULAG experience garnered over a wide range of hydrodynamical multiscale flow simulations indicate that its underlying advective scheme operates as adaptive subgrid model, turning on wherever and whenever dissipation is required to maintain stability, and remaining minimally dissipative otherwise (Smolarkiewicz & Margolin 2007; Piotrowski et al. 2009). Moreover, when it activates, dissipation is concentrated at the smallest scales resolved by the spatial mesh. This allows the simulation of turbulent flows characterized by a spectra with a good inertial range, even on relatively small spatial meshes. This property, in turn, permits long time integration in a reasonable amount of computing time, as needed when simulating fluid systems whose evolution is characterized by widely differing timescales. In this respect solar MHD convection is an extreme case in point, with convective turnover times of hours to days for subphotospheric convection, and a magnetic cycle with a halfperiod of ≃11 yr.
The simulation’s set up used to produce the new data set analyed in what follows (reference name millennium simulation) is similar to that described in Ghizaru et al. (2010); Racine et al. (2011) but with the addition of a small amount of explicit diffusivity (which seems to help stabilizing the solution). It is computed on a modest spatial mesh, 128 × 64 × 47 cells in longitude, latitude and fractional radius, and is characterized by magnetic and viscous Reynolds number in the range 30–60 and a magnetic Prandlt number of order unity. In such ILES, lacking explicit dissipative terms, these quantities can only be estimated indirectly. The range quoted above results from four distinct estimates, namely (1) ratio of injection to dissipation scale; (2) extend of inertial range in the turbulent spectrum; (3) comparison with other LES simulations; (4) numerical experiments whereby increasing explicit dissipation is added to the simulation until the global behavior is altered. The simulation is run over 2.9 × 10^{7} 30 min time steps, adding up to a full simulated time span of 1650 years, in the course of which 41 polarity reversals of the largescale, axisymmetric magnetic component building up in the simulation have taken place at a regular cadence. This amounts to 40 sunspotlike cycles, which present variations in amplitude and periodicity but remained generally wellsynchronized across both hemispheres. We thus define the largescale magnetic field to be associated with our simulated solar cycles as the zonal average of the total magnetic field.
Figure 1A shows timelatitude diagrams of the toroidal (zonallyoriented) largescale magnetic component at the base of the convecting layers (r/R_{⊙} = 0.718). This represents our simulation’s equivalent of the sunspot butterfly diagram. The toroidal component peaks at midlatitudes immediately beneath the base of the convection, reaching 0.6T there. It only shows a hint of equatorward propagation, and remains confined to latitudes higher than ± 30°, in contrast to the sunspot butterfly diagram which builds up at lower latitudes and reaches all the way to the equator. The largescale magnetic field is antisymmetric about the equator, in agreement with Hale’s polarity Laws, but its cycle period is here almost four times larger than in the sun. Figure 1B shows a timelatitude diagram of the near surface radial magnetic component (here at 0.94 r/R_{⊙}). It is characterized by a welldefined dipole moment strongly concentrated in the polar caps, as observed on the sun (but stronger by a factor of ~20) and its polarity reversals occur in phase with the internal toroidal component, unlike in the sun where a π/ 2 phase lag is observed.
Fig. 1 Panel A) Time latitude diagram of the zonallyaveraged toroidal field taken at r/R_{⊙} = 0.718, which marks the interface between the convecting layers and underlying stably stratified region. Panel B) Zonallyaveraged radial field extracted from the simulation’s subsurface layer, at r/R_{⊙} = 0.938. The color scale encodes the magnetic field strength, values in Tesla. Panel C) Evolution of the total magnetic energy integrated over the whole computational domain. The simulation uses a hydrostatic initial state, subjected to small perturbations in fluid velocity and magnetic field. Convection reaches a statistically stationary state after about a year, but the buildup of a stable, cyclic largescale axisymmetric magnetic component sets in after about fifty years, here first in the southern hemisphere. Dipolar parity is maintained here through 20 full magnetic cycles. 
Fig. 2 Meridional plot of the large scale, zonallyaveraged toroidal field (panel A)) and zonallyaveraged radial field (panel B)), extracted at t = 534 yr, a cycle maximum indicated by the vertical dotted black line in Fig. 1. The color scale on the right is in Tesla. The bounding boxes represent the integration regions for the proxies computation, individually for the northern and southern hemispheres. The dashed circular arc corresponds to r/R_{⊙} = 0.718, the interface between the convecting layers and underlying convectively stable fluid. 
Figure 1C shows a time series of total magnetic energy, integrated over the full simulation domain. The cycle of the largescale axisymmetric magnetic field shows up quite prominently as a regular oscillation with welldefined period and amplitude, both showing relatively modest cycletocycle fluctuations.
This simulation presents other features (not shown here) that are reasonably solarlike and similar to those found in previous data sets produced with EULAGMHD, Ghizaru et al. (2010) and Racine et al. (2011). This includes the internal differential rotation profile, and associated torsional oscillations driven by the magnetic cycle, which show the same amplitude and phasing as inferred helioseismically (Beaudoin et al. 2012). The meridional flow also shows cyclic modulation, with the buildup of a secondary flow cell at high latitudes in the descending phase of the cycles, in qualitative agreement with solar surface Doppler measurements (Passos et al. 2012). Moreover, the simulation shows a weak cyclic modulation of the convective energy flux, varying in phase with the magnetic cycle, as does the total solar irradiance (Cossete et al. 2013). While operating in a physical parameter regime still far removed from solar interior conditions, the cyclic dynamics exhibited by these simulations may hold important clues to the operation of the true solar cycle.
Fig. 3 Time series of the fulldisk toroidal field proxy B_{F}. The vertical dashed line represents the beginning of the sequence used for the foregoing analysis, which spans 39 sunspotlike cycles. Cycles are numbered from one minimum to the next, in keeping with the numbering convention of sunspot cycles. 
3. Building proxies for the solar cycle
Our first task is to build, from the simulation output, proxies for observational quantities that are used to characterized the solar cycle, most notably the sunspot number. Our first proxy, hereafter denoted B, is computed by integrating the largescale, (zonally averaged) axisymmetric toroidal field ⟨ B_{φ} ⟩ over a thin layer immediately beneath the base of the convective layers, where the largescale toroidal component reaches peak amplitude (see Fig. 2A). This is physically consistent with current understanding of sunspot formation, which places the formation and storage of the participating flux ropes in the weakly subadiabatic upper reaches of the tachocline (see Fan 2009, Sect. 3, and references therein). The second proxy (poloidal), hereafter denoted A, is aimed at representing the Sun’s dipole moment. We follow the same procedure as before but this time we integrate the zonally averaged radial magnetic field component, ⟨ B_{r} ⟩, in subsurface layers of the polar caps (see Fig. 2B).
Mathematically, the proxies are built using the following expressions
where r_{1} ≃ 0.67R_{⊙}, r_{2} ≃ 0.72R_{⊙}, θ_{1} ≃ 33°, θ_{2} ≃ 73°, and r_{3} ≃ 0.91R_{⊙}, r_{4} ≃ 0.95R_{⊙}, θ_{3} ≃ 62°, θ_{4} ≃ 87° (see masks in Fig. 2).
This procedure is carried out separately for the northern and southern hemispheres, thus yielding four time series B0_{N}, A0_{N} and B0_{S}, A0_{S} correspondingly. We combine the northern and southern toroidal time series in order to build a full disk proxy, B0_{F} = (B0_{N} + B0_{S})/2, since our comparison target, the standard international sunspot number reconstruction, currently does not distinguish solar hemispheres. The toroidal time series are then normalized at the maximum value of B0_{F} yielding the final toroidal proxies: B_{F}, B_{N} and B_{S}. A similar procedure is used to build the poloidal proxies, A_{N} and A_{S}. In order to facilitate the identification of cycles features (maxima, minima, etc.), these time series are smoothed by convolving them with a Gaussian filter (width = 24 months for the toroidal proxies and 36 months for poloidal proxies), following a procedure analogous to that described by MuñozJaramillo et al. (2013a).
4. Full disk analysis
We first investigate the statistical properties of cycles as defined through our full disk proxy. The associated time series for the toroidal proxy B_{F} is shown in Fig. 3. Significant variations are seen in the cycle amplitude, as well as from one minimum to the next. For all the cycles, maxima and minima along with their times of occurrence were identified, and then used to define other quantities. This is done semiautomatically using an analysis pipeline designed in Mathematica for this purpose. At first the analysis pipeline smooths the data slightly and computes its first derivative. The zeros of the derivative are used to identify local maxima that are a posteriori classified into cycle maxima and cycle minima. In the places where the algorithm fails to recognized these points, a manual correction is introduced. Once these points are identified all other quantities are straightforwardly computed.
For comparison with our full disk toroidal proxy, B_{F}, we chose the international sunspot number, SSN. We opted to use the monthly averaged SSN (from 1749 to Dec. 2012) data because it remains the most commonly used proxy of solar magnetic activity and it has the same onemonth temporal cadence as our computational proxies. We retrieved the SSN data from the SIDC^{1} database and ran it through the same analysis pipeline used to characterize our proxies. Although the sunspot number time series has been analysed exhaustively in countless works, we decided to reanalyse it again using our pipeline to avoid biases introduced by different analysis methodologies or temporal intervals, so as to maximize statistical homogeneity when comparing the results obtained for our proxies with those obtained for observational data (Dikpati et al. 2008; Ogurtsov & Lindholm 2011).
Fig. 4 Histograms of frequency counts for the cycle periods for our fulldisk toroidal proxy B_{F} (panel A)) and the smoothed monthly SSN (panel B)). The vertical dashed lines represent the distributions means at 40.5 and 11 years correspondingly. Although the means differ by a factor of almost four, the two distributions show similar standard deviations and skewness (see text). 
We adopt the usual cycle numbering convection, from one minimum to the next, with the minimum of cycle N defined as the minimum that follows cycle N maximum. The other quantities defined are: maxima are peak amplitudes of the cycles; minima are the cycles’ amplitudes at minimum; period is the cycle’s duration defined as the time between two consecutive minima; rising time is the time between minimum N − 1 and maximum of cycle N. Sometimes around solar maxima the cycle shows 2 or 3 distinctive local peaks and the highest one is not always the first, conditioning the rise time; mean growth rate is the mean value of 1st derivative of B averaged over the rising time. This is our way of gauging the steepness of the cycle and is similar to the method presented in Cameron & Schüssler (2008). This yields a better measure of how intense the rising phase of the cycle is. While we still use the times of minimum and maximum to define the averaging period for this quantity, the averaging procedure actually compensates for cycles with double or triple peaks. Between these peaks the derivative’s sign changes and averages out to small very values after the first peak. This avoids problems related to the correct definition the rising phase if the highest local peak is not the first. For more details consult Cameron & Schüssler (2008) and Karak & Choudhuri (2011).
Our first analysis step is to seek correlations between these various quantities. Table 1 compiles some of the results of this exercise, for our toroidal proxy as well as for the smoothed monthly SSN (rightmost column). See the Appendix for a complete list of the correlations calculated. The results obtained for the latter using our analysis pipeline are compatible to the recent analysis by Ogurtsov & Lindholm (2011) for the same data and time interval.
Comparison between the Pearson’s linear correlation coefficients, r, obtained between several quantities measured at cycles N and N + 1 for our simulated B_{F} (39 cycles) and the observed SSN (23 cycles).
Examination of Table 1 reveals that simulation and observations show some common correlations, but also some noticeable differences. Some of the correlations listed in this table are commonly known as Waldmeier rules and are often used to characterize the solar cycle. For example the relationship between cycle duration (period) and peak amplitude (maxima) within the same cycle presents similar values in both data sets. In contrast, the period of cycle N presents a high (anti)correlation with N + 1 maxima (i.e. the amplitude of the following cycle) for the SSN, but a much weaker correlation for our B_{F} proxy. The same happens for rise time vs. maxima. Nevertheless for most of the relationships the sign of the correlations agree.
Besides the obvious fact that our simulation is not the real Sun, and even taking in consideration the absence of surface dynamical phenomena, there are additional facts to keep in mind when comparing the results obtained for both data sets. Part of the differences in the values of the correlations can be attributed to the different number of cycles (low in both cases) used to compute these statistical properties. Moreover, the deep toroidal field which we integrate to define our B_{F} proxy shows very little latitudinal propagation, and therefore the associated magnetic flux systems of successive cycles exhibits no temporal overlap around the times of cycle minima. For the SSN this is not the case, as there exists an overlap of the decaying phase of cycle N − 1 with the rising phase of cycle N (Cameron & Schüssler 2007). This overlap affects the timing of SSN minima, and in some cases also the onset of the rising phase of the new cycle. Such effects may also introduce some differences in the computed correlations.
The cycles measured in our toroidal B_{F} proxy have an average period of 40.5 yr, almost four times the 11.0 yr average period of the SSN. The frequency distribution of proxy cycle periods, shown in Fig. 4A, has a mildly asymmetrical shape similar to that constructed from the SSN time series, plotted on 4B. More specifically, both have similar skewness, S = 0.57 for B_{F}, vs. S = 0.64 for the SSN. Both also have similar standard deviations, σ = 1.5 yr for B_{F} and σ = 1.2 yr for SSN. However, with the B_{F} mean period four times longer, comparable standard deviations indicate that the period shows better relative stability in our simulation than in the sun.
Another feature present in both our proxy and the SSN is the presence of longterm modulation of the cycle amplitude. In the case of the SSN, the bestknown such long periodicity is the socalled Gleissberg cycle with period ~90 yr, with even longer modulation periods identified in the cosmogenic isotope record (Usoskin 2013). Figure 5 shows power spectra of our fulldisk toroidal proxy B_{F} (A) and SSN (B).
Fig. 5 Panels A) and B) Power spectrum for B_{F} and for SSN respectively. The dashed vertical lines mark important frequencies in the spectra and are defined, from left to right as ω_{1}, ω_{2}, ω_{3}, ω_{4} and ω_{5} (the latter is not highlighted for the SSN). 
Fig. 6 Histograms of frequency counts for the cycle peak amplitude for A) our fulldisk toroidal proxy B_{F} and B) the SSN. The vertical dashed lines represent the distributions means. 
The primary peak, labelled ω_{3}, is of course the main cycle (~40 yr for B_{F} and 11 yr for the Schwabe SSN cycle). Besides this primary spectral peak in both distributions, we have highlighted on their right two other frequencies, ω_{4} and ω_{5}, that correspond to the first and second harmonics of the main peak of B_{F}. For the SSN, ω_{5} is below noise level and its not identified by eye. These higher harmonics are simply reflecting the fact that in both cases the cycles are not pure sinusoids. On the lower frequency side of the power spectra, we have also highlighted two peaks as ω_{1} and ω_{2} which correspond to 279 yr and 170 yr, for B_{F} and 89 yr and 53 yr for SSN, respectively. These peaks are not welldefined spectrally due to the limited length of the time series, but are statistically relevant as they amount to the same power as the first harmonic of the main peak. In the SSN record, the longer period is usually associated with the Gleissberg cycle, longer than the Schwabe cycle by a factor ~10. This periodicity in the SSN is usually considered wellestablished, as it shows up in other records of solar activity with long temporal records. Interestingly, the longer period found in our B_{F} proxy is also a factor of ~10 longer than the primary cycle. Whether this is truly significant will require, at the very least, pushing the simulation much farther in time to better spectrally resolve the longer periodicities apparent in Fig. 5A.
Figure 6 shows the frequency count histograms for cycle maxima in the B_{F} proxy in panel A, and for the SSN in panel B. Both distribution are reasonably symmetrical about their mean value, but for the SSN the relative spread is larger, with the standard deviation amounting to 37% of the mean, vs. 16% for our B_{F} proxy.
Another classical SSN pattern is the socalled GnevyshevOhl rule, whereby evennumbered sunspot cycles tend to have lower amplitudes than their oddnumbered counterparts, which leads to an alternance between higherthanaverage and lowerthanaverage cycle amplitudes. This pattern is best sought in our proxy time series by first establishing a temporal sequence of cycle maxima (B_{F})_{k} which is then detrended by subtracting its own 121 running mean: (3)Figure 7A illustrates this procedure, with the black solid and dashed lines showing the raw and 121 smoothed sequences of maxima, and panel B the result of the above detrending procedure. The GnevyshevOhl pattern shows up here as a regular alternance of positive and negative detrended values, a long such sequence materializing here between simulated cycles 16 and 30. Panel C shows the result of applying the same procedure to the smoothed monthly SSN, with the GnevyshevOhl pattern sustained here from cycle 9 to 21.
Fig. 7 The GnevyshevOhl pattern in our fulldisk toroidal proxy. Panel A) shows the sequence of cycle maxima (black), and the corresponding 121 running mean thereof (dashed); subtracting the latter yields the detrended sequence plotted in panel B). The GO pattern shows up here as a sustained regular alternance between positive and negative detrended values, e.g. between simulated cycles 16 and 30 (thick segments). Panel C) shows a similarly detrended cycle amplitude sequence for the monthly smoothed SSN; here the GO pattern holds uninterrupted from sunspot cycle 9 to 21 (thick segments). 
5. Hemispheric analysis
The existence of continuous, systematic observation campaigns, some beginning over a century ago, has made it possible to characterize solar activity individually for each hemisphere, a distinction that has not yet been made for the complete SSN data. Hemispheric asymmetries are found in the sunspot number and areas (Temmer et al. 2006; Norton & Gallagher 2010; MuñozJaramillo et al. 2013a), in polar fields (Svalgaard et al. 2005; Schatten 2005; MuñozJaramillo et al. 2012) and in many other solar activity indexes. Besides looking at the individual characteristics of each hemispheric proxy, we also intend to establish the level of crosshemispheric coupling in our simulation. Figure 8 shows the time series of northern (black) and southern (red) hemispheric toroidal (panel A) and poloidal (panel B) proxies. Southern hemisphere proxies are arbitrarily assigned a negative value, for plotting purposes. For comparison with our toroidal hemispheric proxies, B_{N} and B_{S} (Fig. 8A), we chose the hemispheric sunspot area (SSA) data recently compiled and calibrated by MuñozJaramillo et al. (2013a) and encompassing 15 complete sunspot cycles (from 1835 to 2010). Similarly, for comparison with our poloidal proxies, A_{N} and A_{S} (Fig. 8B), we use the absolute value of the signed polar flux (PF) time series from MuñozJaramillo et al. (2012), spanning 10 activity cycles from 1907 to 2010. The SSA time series is smoothed in the same way as our toroidal proxy (24 month Gaussian filtering) while the PF time series is not subjected to smoothing because of the low resolution available (1 point per year). These observational time series were also processed by our analysis pipeline.
Fig. 8 Solar cycle proxies built from the millennium simulation. Black and red lines represent magnetic activity on the northern and southern hemispheres correspondingly. Panel A) shows the toroidal proxies, B_{N} and −B_{S}. In panel B) we present the poloidal (dipolar) proxies, A_{N} and −A_{S}. The southern proxies are plotted with a negative sign to better illustrate the coupling differences between both hemispheres. 
For each hemisphere, we segment the proxy time series into cycles, following the same strategy and numbering convention introduced earlier for the fulldisk proxy. Just as in the previous section, we first examine the correlations found between various cycle characteristics obtained for our simulated hemispheric toroidal proxies and the SSA time series. The main results are listed in Table 2, with complete correlation tables again presented in the Appendix.
Comparison between the correlation coefficients obtained between several quantities for the millennium simulation toroidal proxies B_{N} and B_{S}, and observational SSA data for the north and south.
As in the case of the analysis carried out for the fulldisk proxies (viz. Table 1), we observe both similarities and differences in the levels of correlation obtained for the simulationbased proxies, as opposed to the SSA, although once again in most cases all correlations have the same sign. At this point we would like to note that the SSN and SSA time series have different statistical properties (Norton & Gallagher 2010). This should be taken into consideration before attempting to make a direct comparison between the SSA correlations in this section and SSN correlations on the previous one.
Figure 9 shows correlation plots for four of the most noteworthy correlations found in the proxies B_{N} and B_{S}, colorcoded according to hemisphere. The period vs. maxima anticorrelation (Fig. 9A) is again present at a somewhat stronger level than in SSA data, as was the case with the fulldisk proxy B_{F} as compared to the SSN. In the simulation as well as in observed data, the strongest correlation found is between mean growth rate vs. maxima (Fig. 9B), an interesting feature for the shortterm prediction of cycle peak amplitude, also present in our fulldisk toroidal proxy as well as in the SSN (see Table 1). A reasonably good correlation is also found between period vs. rise time (Fig. 9C), stronger in the simulation than in SSA data. Rise time anticorrelates moderately well with maxima in the northern hemisphere but not in the south (Fig. 9D); no significant such anticorrelation is found in SSA data.
Fig. 9 Waldmeiertype relationships for the magnetic proxies built from the millennium simulation. A) Period vs. Maxima; B) Mean Growth Rate vs. Maxima; C) Period vs. Rise Time; and D) Rise Time vs. Maxima. The similar anticorrelation observed in Period vs. Maxima (in A)) and Rise Time vs. Maxima (in D)) is a consequence of the good correlation between Period and Rise Time (in C)). 
5.1. More on toroidal proxies: B_{N} vs. B_{S}
Examination of Fig. 8A reveals that the northern and southern hemispheric toroidal proxies behave similarly, except very early on as for this specific simulation the largescale magnetic field builds up sooner in the southern hemisphere (see also Fig. 1). Cycles in both hemispheres remain very well synchronized over more than 1600yr, and show nearly identical mean periods over that time span. Figure 10 shows the frequency distributions of maxima for the northern and southern toroidal proxies. Both distributions have similar means, while the southern hemisphere shows a slightly larger standard deviation (σ_{BN} = 0.163 and σ_{BS} = 0.210). Interestingly, both hemispheres exhibit bimodality, in the form of a secondary peak of lowamplitude cycles well separated from the group of highamplitude cycles – to the extent that such a separation can be confidently ascertained based on statistics from 39 cycles in each hemisphere. Interestingly, recent reconstructions on longterm solar activity from the abundances of cosmogenic radioisotopes also show a similar bimodal distribution, characterized by a secondary population of lowactivity values (Usoskin et al. 2014).
Fig. 10 Histograms of the cycle peak amplitudes, based on the toroidal proxies, measured in the A) north and B) south hemispheres. Although the distributions are built from a relatively small sample of cycle maxima, both hemispheres show a clear hint of bimodality, in the form of a subgroup of lowamplitude cycles. Vertical dashed lines indicate mean values 
Despite the overall correlation between the complete time series of B_{N} and B_{S} being a reasonable r = 0.66 with significance level P = 99.9%, and their similar bimodal statistical distributions of amplitudes, the cycletocycle variations of peak amplitudes seem to unfold almost independently in the northern and southern hemispheres. In particular, one may verify, upon examination of Fig. 8A, that the lowamplitude cycles (with amplitudes below the mean value minus 1σ) belonging to the second mode of the amplitude distributions never occur simultaneously in both hemisphere. In fact lowamplitude cycles in one hemisphere rather tend to be accompanied by a high amplitude cycle in the other, as also shown by a correlation coefficient r = −0.25 when correlating the maxima of one hemisphere against the other. This is further illustrated in Fig. 11. The dashed lines indicate the (normalized) amplitude values B_{N} = B_{S} = 0.6, which is approximately the mean minus one standard deviation for either hemispheric distribution. We adopt this value as delineating the two modes in the distributions of Fig. 10. Please note how the lower left quadrant segmented by these two lines is empty; now, with B_{N} ≤ 0.6 for 8/39 cycles, and B_{S} ≤ 0.6 for 7/39 cycles, the probability of none of 39 cycles having B_{N} ≤ 0.6 and B_{S} ≤ 0.6 is 0.23; this is certainly not minuscule, but hints at a possible crosshemispheric influence. It also explains why the bimodality so clear in Fig. 10 – again, smallN statistics notwithstanding, – does not show up in Fig. 6A.
Fig. 11 Correlation between northern and southern maxima A) and period B) in the time series of toroidal proxies B_{N} and B_{S}. In the two panels the thick blue cross indicates the mean values ±one standard deviation. In A) the horizontal and vertical dashed lines, drawn at B_{N} = 0.6 and B_{S} = 0.6, delineate the cycles belonging to the two modes of the distributions in Fig. 10. The diagonal light gray line indicates a perfect correlation (B_{N} = B_{S}). 
Turning to proxy amplitude at cycle minima, no crosshemispheric correlation is found (r = −0.07 with P = 32%). For the SSA time series the crosshemispheric amplitude correlations are much higher, with r = 0.74 and P = 99% for maxima and r = 0.5 with P = 96% for minima.
The mean cycle period in the northern and southern hemispheres, 40.65 yr and 40.57 yr, fall well within one another’s standard deviations (σ_{N} = 1.81 yr and σ_{S} = 1.65 yr, respectively). Both hemispheres thus show the same primary periodicity, but their frequency distributions of periods, shown in Fig. 12, differs more than their amplitude counterpart. Even though the southern hemisphere distribution is characterized by a narrower span than in the north, it shows again bimodality, with a group of 5 longperiod cycles separated from a welldefined main group by a gap around 42.5 yr (see Fig. 11B). We note that this subgroup of longduration cycles does not map onto the lowamplitude group of cycles in the peak amplitude distribution; in fact, examination of Fig. 9A reveals that the lowamplitude cycles (B_{S} ≤ 0.6) in the southern hemisphere span the whole range of measured cycle periods, with the two lowest amplitude cycles falling at opposite ends of the period distribution. Unlike with cycle amplitude, the scatter plot of southern vs. northern cycle period (Fig. 11B) shows no obvious deficit in any quadrant defined through the mean values and standard deviation.
Fig. 12 Histograms for the periods measured in the A) north and B) south hemispheres. Signs of bimodality are apparent again here, although not as welldefined in the northern hemisphere as in the south. However, these subgroups of longperiod cycles do not map onto the subgroups of lowamplitude cycles in Fig. 10 (see text). 
Fig. 13 The power spectrum for A) B_{N} and for B) B_{S}. The dashed vertical lines mark important frequencies in the spectra and are defined, from left to right as ω_{1}, ω_{2}, ω_{3}, ω_{4} and ω_{5}. 
Fig. 14 Time lag measured between the two hemispheres for minima (blue, dashed line) and maxima (black, thick line). Negative values indicate that the northern hemisphere cycle leads its southern counterpart. The numbers inside the squares represent the simulated cycle number. The mean lag values are −2.6 yr for minima and −2.3 yr for maxima, amounting to about 6% of the mean cycle duration. The largest individual lag values approach 10 years for cycles 3, 20 and 33, i.e. 25% of the mean cycle duration. 
As in the previous section, we test our hemispheric toroidal proxies for longer periodicities. Again the power spectrum of both proxies show the same characteristics as the power spectrum for B_{F} (Fig. 13). While the main peak and its first and second harmonics (ω_{3}, ω_{4} and ω_{5}), are practically indistinguishable between both hemispheres, the lower frequencies show some noticeable differences. For B_{N}, ω_{1} and ω_{2} correspond to periods of 339 yr and 210 yr while for B_{S} they correspond to 279 yr and 170 yr. Since these frequencies in the southern hemisphere present a higher power than their counterpart in the north, they dominate the B_{F} power spectrum in this frequency range. Once again, longer simulated time series are required to better resolve these tantalizingly solarlike long periodicities.
5.2. Lag analysis
The northern and southern hemisphere cycles in Fig. 8 are clearly wellsynchronized, yet careful examination reveals that minima and maxima do not always occur at exactly the same times in both hemispheres, with a sustained lag sometimes extending over many successive cycles, e.g. between simulation years 800 and 1000. In order to quantify these lag effects we define the following two additional quantities: maximum lag, as the difference in time between maxima in the north vs. the south, i.e. time of maxima in the north minus the time of maxima in the south. A negative value means that the maximum was reached first in the north then in the south; minimum lag is defined the same way but using minima.
In order to estimate the average lag between north and south we conduct a minima and maxima time lag analysis on the hemispheric toroidal proxies. The time lags computed for the 39 cycles are presented in Fig. 14. On average, the northern hemisphere leads over the south by 2.3 years for maxima and 2.6 years for minima. Coincidentally for the SSA data, the northern hemispheres also leads over the south by 0.5 years for maxima and 0.08 years for minima. The relative lags, measured as a fraction of the cycle period, are thus similar for maxima, but larger in the simulation for minima. The correlation between these minima and maxima lag times for our proxies is r = 0.8 with P = 99.9% while for the SSA data is r = −0.17 with P = 49%. The low correlation between these lags in the SSA data is mainly attributed to the fact that around cycle maxima, the data show several peaks. Since the maxima is defined as the highest peak amplitude, the fact that often the first peak is not the strongest, heavily influences the statistics of all parameters dependent on maxima and its timing. This also explains the main differences between correlations obtained for SSN and SSA. For a more detailed discussion of this subject consult Dikpati et al. (2008); Cameron & Schüssler (2007, 2008).
In the simulation run analysed here, the southern hemisphere develops a largescale cycle first, before the northern hemisphere joins in (see cycle 0 in Fig. 8). Nevertheless, once cycles have stabilized in both hemispheres (starting in cycle 1), the activity in the northern hemisphere ends up leading throughout most of the simulation. The slow, erratic evolution of the lag with time, apparent in Fig. 14, is reminiscent of a random walk. Yet, considering the low degree of crosshemispheric correlation in the amplitude of the cycles, it remains noteworthy here that after 40 cycles, the accumulated hemispheric lag remains below the 10% level. We performed a Monte Carlo simulation where two independent 39member sequences of cycle periods are extracted from Gaussian distributions with the mean and standard deviations of the hemispheric distributions of Fig. 12, and calculated the probability distribution of final lags at the end of 10^{5} such pairs of sequences. The probability of this lag falling within ± 10% of one mean period is ≃0.2, which is small though again certainly not minuscule.
Correlations between the two simulation proxies, B vs. A for both hemispheres and for the observational data, sunspot areas vs. polar flux timeseries.
We repeated the lag analysis using this time the poloidal proxies plotted in Fig. 8B. The analysis is more arduous because that proxy is noisier, with a greater tendency to exhibit multiple peaks within one cycle, and multiple minima in between two cycles which forces specific extrema to be chosen for the lag analysis. The minima lag ends up again at ~10% of the mean period after 39 cycles, but the lag based on maxima wanders about zero lag to within ± one year for 20 of the 30 last cycles of the simulation, with a final lag of about 9 months; the same Monte Carlo exercise as described above now yields a probability of ≃0.04 for such a low final lag to materialize after 39 cycles.
The lag analysis thus supports the notion that despite the low degree of hemispheric correlation characterizing cycle amplitude, some selfregulating process is keeping the two hemisphere in steps with regards to polarity reversals. We shall return to this intriguing issue in Sect. 6 below.
5.3. GnevyshevOhl pattern
The GnevyshevOhl pattern of alternating higher and lowerthanaverage peak cycle amplitude, characterizing our full disk proxy B_{F} also materialize in the individual hemispheric proxies, but the patterns are not necessarily coincident across hemispheres, nor in phase; for example, from cycles 1 to 8 the southern hemisphere shows a continuous GO pattern, while a shorter, outofphase GO pattern is sustained from cycle 5 to 9 in the north, with the consequence that the fulldisk proxy only shows a brief GO episode from cycle 1 to 4 (viz. Fig. 7). A similar situation occurs at the end of the simulation, leaving only a 5cyclelong GO pattern in the full disk proxy. On the other hand, the 15cyclelong GO pattern observed between cycles 16 and 30 (thicker line segments in Fig. 7) is largely due again to the southern hemisphere, but now with the northern hemisphere making a shorter inphase contribution from cycles 22 to 27. The GO pattern is clearly a property of the dynamo operating in this simulation, but shows little synchrony or phase coherence across hemispheres. This points again to a significant hemispheric autonomy in whichever physical process sets the cycle’s amplitudes.
5.4. Toroidal vs. poloidal
Recent years have witnessed increased interest in the relationships between various proxies of the solar poloidal magnetic field, and proxies of the internal toroidal field, because of the potential use of the former as a predictive precursor of the latter. As we now turn in this section to the relationships existing between our toroidal and poloidal proxies, two caveats are in order. Often in the literature, we find references to a significant positive correlation between the intensity of the polar flux of cycle N and the amplitude of sunspot cycle N+1 (Schatten 2005; Svalgaard et al. 2005; MuñozJaramillo et al. 2013a). In these works, this PF intensity (sometimes called average amplitude of the polar flux during minimum) is an integrated measure of the PF over an interval of a few years around minimum. This is different from what we define as amplitude of PF here. In our case we are simply referring to the cycle’s maximum amplitude of the PF and not some kind of integrated value. We opted for this definition because our A proxies are synchronized (and in phase) with our B proxies, so we cannot compute any meaningful integrated averaging during minimum. Second, since only 10 cycles are available in the PF data, the simple statistical quantities derived in this work should be taken as purely indicative. For a more thorough study of the relationships between these observational poloidal and toroidal proxies see MuñozJaramillo et al. (2013a,b).
In Table 3 we present some of the correlations found between the toroidal and poloidal proxies in this study.
Despite (or perhaps because of) their noisier character, cycles defined on the basis of the hemispheric poloidal proxies are in some ways better behaved than those defined through the toroidal proxies. The amplitude distributions are fairly symmetrical about their respective means and show no bimodality. Examination of Fig. 8B also reveals that the two poloidal proxies also show somewhat better correlation across hemispheres, with lowamplitude cycles now often occurring simultaneously, e.g, cycles 3, 8, 31 , 32 (but not always, cf. cycles 7, 22 and 33). The two time series have a higher correlation coefficient, r = 0.81, than do their toroidal counterpart (r = 0.66). The period distributions are more compact than in Fig. 12, with nearly identical means (40.50 yr in the north, 40.38 yr in the south), but show more extended lowamplitude tails extending ± 10 yr on either sides of the mean, twice the span of the period distributions for the hemispheric toroidal proxies (viz. Fig. 12). These tails also translate into larger standard deviations, σ = 3.3 yr for the north and σ = 4.6 yr for the south. Of the Waldmeierlike correlations of Fig. 9, only the mean growth rate vs. cycle maxima remains significant (see Table 3), the other being effectively erased by the higher fluctuation levels characterizing the poloidal proxies.
6. Discussion
Some of the correlations and fluctuations patterns uncovered in our simulation carry interesting information regarding the operation of the magnetic cycle developing therein.
Generally speaking, correlations between rise time, growth rate, amplitude and period as the cycle fluctuates are indicative of selfsimilarity in the unfolding of individual magnetic cycles. These patterns are present in our full disk proxies, as well as in individual hemispheres. These correlation patterns can also be found in stochastically forced meanfield dynamo models (e.g. Karak & Choudhuri 2011 or Pipin & Sokoloff 2011). Our simulation achieves the same thing through dynamicallydriven fluctuations.
The GnevyshevOhl pattern is often associated with the possible presence of a steady fossil magnetic field present in the deep solar interior, systematically offsetting the cycle produced by dynamo action in the convection zone (see Charbonneau et al. 2007, and discussion therein). In our simulation, there is no fossil field introduced in the stably stratified fluid layers underlying the convecting layers. Instead, in the low dissipation environment of this tachoclinelike layer, any amplitude bias introduced by a higherthanaverage cycle takes a long time to dissipate locally; a larger fraction of magnetic flux building up in subsequent cycle is thus spent cancelling the flux of the preceding cycle, resulting in lowerthanaverage cycle. Then, the next following cycle has less oppositepolarity magnetic flux to reverse, and so will end up with a higherthanaverage amplitude once again. The presence of a GnevyshevOhllike pattern is thus a direct consequence of the low dissipation characterizing this simulation.
Similar trends, correlations and distributions are obtained for toroidal proxy integration domain extending to higher and lower latitudes, and including most of SCZ. In particular, the GnevyshevOhl pattern and bimodal distributions of hemispheric cycle amplitudes are both robust with respect to such alternate definitions of our sunspot number proxies.
The timing of polarity reversals is very well synchronized across hemispheres in our simulation, suggesting good hemispheric coupling; on the other hand, hemispheric cycle amplitudes based on our toroidal proxy are essentially uncorrelated, indicative of poor hemispheric coupling. How can these two conflicting facts be reconciled?
In classical meanfield dynamo models of the solar cycle, only magnetic diffusion operates to couple the hemispheres. Turbulent diffusion being essentially a free parameter in such models, it is usually possible to produce the desired level of coupling by proper adjustment of this parameter. In convection simulations aiming at the solar parameter regime, magnetic dissipation is kept as small as possible, but other coupling mechanism are present. In particular, most of the simulations of the type considered here produce, in the lowlatitude portion of the convection zone, elongated convective rolls parallel to the rotation axis and stacked in longitude (see Miesch & Toomre 2009, Sect. 2.2). These rolls have alternating sense of twist, and corresponding alternate direction of flow along their axis, so as to yield the same kinetic helicity in each hemisphere. As a consequence plasma is being continuously exchanged across the equatorial plane, mixing any magnetic field located at low latitudes on either side of the equator, and thus achieving coupling. Examination of the simulation also reveals significant poloidal magnetic flux threading from one polar cap to the other, along the tachocline and through the equatorial plane, which can also provide crosshemispheric coupling of the dipole moment.
Under this view, the anomaly is thus the low correlation observed in the toroidal proxy amplitude. The regions of strong internal toroidal fields are located at relatively high latitudes in our simulation (viz. Fig. 1), and thus should suffer little crossequatorial mixing. The fact that cycles retain the same overall morphology, and in particular the same latitudinal distributions, would argue against the excitation of a quadrupolar (or higherorder) mode, as seems to be happening in some of the similar simulations reported upon recently in Augustson et al. (2014), for example.
Alternately, amplitude asymmetries could instead result from local processes operating independently in each hemisphere. The observed bimodality in cycle amplitude offers a hint at a possible explanation, namely that two dynamical attractors may exist in the simulation’s phase space. We have explored this possibility by producing various phase space trajectories based on global proxies of the magnetic fields as well as largescale flows developing in the simulation. Morphologically speaking, trajectories associated with lowamplitude cycles are similar to their higher amplitude counterparts.
Another possibility is that the lowamplitude cycles results from the onset of a global instability of the strong toroidal magnetic field building up in the tachocline. In some portions of parameter space, the EULAGMHD simulations such as the one studied here, undergo a transition to a nonaxisymmetric largescale dynamo mode, still characterized by regular polarity reversals on a similar multidecadal cadence, but symmetric about an axis inclined with respect to the rotation axis. Analyses carried out to date indicate that this destabilization may be due to the growth of a MHD instability, the exact nature of which having not yet been identified with confidence (Lawson et al., in prep.). It is possible that the family of low amplitude cycles results when the buildup of this instability occasionally perturbs the axisymmetric cycle, while failing to destabilize it altogether. Because this instability can develops on single toroidal field bands, there is no reason to expect that failed onsets should happen simultaneously in both hemispheres. This nonaxisymmetric dynamo mode could then be the soughtafter second attractor in the simulation’s phase space. Whether this could also explain the similar bimodality uncovered by Usoskin et al. (2014) remains to be seen.
7. Conclusions
We have presented in this paper a EULAGMHD simulation of solar convection developing a largescale axisymmetric magnetic field component undergoing reasonably solarlike polarity reversals. Forty such reversals, occurring at a regular ~40 yr cadence, span 1650yr of total simulation time. We compile from the simulation output various measures of these magnetic cycles, which we then compare to observational measures and proxies of solar activity, namely the sunspot number and areas time series, and the polar flux inferred from faculae counts.
Many observationallyinferred solarlike fluctuations patterns in amplitude and duration of our simulated cycles are recovered, including Waldmeierlawlike correlations between cycle parameters (rise time, period, amplitude), hints of long, Gleissberglike periodicities in cycle amplitude, as well as the GnevyshevOhl pattern of extended sequences of alternating higher and lowerthanaverage successive cycle amplitudes.
Cycles in both hemispheres are wellcoupled with regards to timing of polarity inversions, and show a cumulative lag of only ≃10% of the cycle period after 39 cycles, based on toroidal proxy, and 4% of mean period based on our surface poloidal field proxy. Such good synchrony also characterizes the solar activity record. On the other hand, our simulated cycles show weak hemispheric correlation in cycle amplitude, with a much stronger correlation observed in the solar cycle. Moreover, cycle amplitudes show a bimodal distribution, with a low amplitude mode accounting for ≃20% of cycles; these weak cycles never occur simultaneously in both hemispheres, which suggests some some level of interaction between hemispheres. No significant correlation could be found between the occurrence of these lowamplitude cycle and hemispheric lag in previous, current or subsequent cycle.
It is noteworthy that solarlike fluctuations patterns in cycle amplitude, duration and hemispheric lag can be produced in a simulation which does not incorporate the contribution of active region decay to the reversal of the surface dipole moment. Here the dipole moment building up in the simulation is entirely driven by the turbulent electromotive force. Of course this does not preclude additional contributions to the flux budget of the polar caps, such as provided by active region decays, but clearly this regeneration mechanism for the poloidal component is not required to sustain global dynamo action in this simulation. This has a number of implications for solar cycle prediction; in our simulation, variations of the amplitude and duration of successive cycles are driven stochastically by turbulence on small spatial scales, operating within the solar convection zone. If this is the dominant source of variability, then one would expect very limited predictive capabilities. Observational evidence pointing in the other direction are presented in MuñozJaramillo et al. (2013a,b) where the amplitude of a cycle seems to be directly related to the amount of polar flux resulting from active regions decay during the previous minimum. Indeed we could be faced with a scenario where the amount of polar flux is a superposition of two contributions, one coming from the inner dynamo and the other from active regions decay. If a fraction of the magnetic flux associated with the surface dipole finds its way to the tachocline, e.g. through submergence and transport by a largescale meridional flow, then some modulation of the cycle amplitude can still materialize (see Charbonneau & Barlet 2011 for specific examples). This is a matter that we are exploring in ongoing analyses of the simulation. A crucial question in dynamo modelling thus remains: how important is active region decay to the operation of the global solar cycle?
Online material
Appendix A: Complete correlation tables
In this Appendix we present the complete Pearson’s correlation tables, between all simulated proxies, with the significance level P color coded in the background of each cell for an easier reading. Values with r> 0.5 are displayed in bold format.
Fig. A.1 Color coding used to show different values in the correlation significance level P. Significance below 70% has a white background. 
Fig. A.2 Complete correlation table for the full disk proxy B_{F}, including correlations with a full disk poloidal proxy, A_{F}, which were not mentioned in the main text or included in Table 1. 
Fig. A.3 Complete correlation table between cycle measures in the northern hemisphere (cf. Table 2 and Sect. 5.4). 
Fig. A.4 Complete correlation table between cycle measures in the southern hemisphere (cf. Table 2 and Sect. 5.4). 
Fig. A.5 Complete correlation table between cycle measures in the northern vs. southern hemispheres. 
Fig. A.6 Complete correlation table between lags and cycle measures in the northern and southern hemispheres. 
Solar Influences Data Center at http://sidc.oma.be
Acknowledgments
We wish to acknowledge the contribution of Mihai Ghizaru in the initial design of the EULAGMHD simulation used in this paper, Mark Miesch for coining the name millennium simulation, as well as Nicolas Lawson for his punctual contributions to various exploratory analysis of the simulation output. We are also thankful to Andrés MuñozJaramillo for kindly providing us the SSA and PF observational data. Computing time for this project was granted through CalculCanada and simulations carried out on the computing platforms of CalculQuébec. D. Passos acknowledges support from the Fundação para a Ciência e Tecnologia (FCT) grant SFRH/BPD/68409/2010, CENTRAIST and the University of the Algarve for providing office space. P. Charbonneau is supported by the Natural Sciences and Engineering Research Council of Canada.
References
 Augustson, K., Brun, A. S., Miesch, M. S., & Toomre, J. 2014, ApJL, submitted [arXiv:1310.8417] [Google Scholar]
 Beaudoin, P., Charbonneau, P., Racine, E., & Smolarkiewicz, P. K. 2012, Sol. Phys., 282, 335 [Google Scholar]
 Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2010, ApJ, 711, 424 [Google Scholar]
 Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R., & Schüssler, M. 2007, ApJ, 659, 801 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, R., & Schüssler, M. 2008, ApJ, 685, 1291 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P. 2010, Liv. Rev. Sol. Phys., 7, 3 [Google Scholar]
 Charbonneau, P., & Barlet, G. 2011, J. Atm. Sol. Terrestr. Phys., 73, 198 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P., & Smolarkiewicz, P. K. 2013, Science, 340, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, P., Beaubien, G., & StJean, C. 2007, ApJ, 658, 657 [NASA ADS] [CrossRef] [Google Scholar]
 Cossete, J.F., Charbonneau, P., & Smolarkiewicz, P. K. 2013, ApJ, 777, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Dikpati, M., Gilman, P., & de Toma, G. 2008, ApJ, 673, L99 [NASA ADS] [CrossRef] [Google Scholar]
 Eddy, J. A. 1976, Science, 192, 1189 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Fan, Y. 2009, Liv. Rev. Sol. Phys., 6, 4 [Google Scholar]
 Ghizaru, M., Charbonneau, P., & Smolarkiewicz, P. K. 2010, ApJ, 715, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Hathaway, D. H. 2010, Liv. Rev. Sol. Phys., 7, 1 [Google Scholar]
 Käpylä, P. J., Korpi, M. J., Brandenburg, A., Mitra, D., & Tavakol, R. 2010, Astron. Nachr., 331, 73 [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., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Karak, B., & Choudhuri, A. 2011, MNRAS, 410, 1503 [NASA ADS] [Google Scholar]
 Miesch, M. S., & Toomre, J. 2009, Ann. Rev. Fluid Mech., 41, 317 [NASA ADS] [CrossRef] [Google Scholar]
 MuñozJaramillo, A., Sheeley, N. R., Zhang, J., & Deluca, E. 2012, ApJ, 753, 146 [NASA ADS] [CrossRef] [Google Scholar]
 MuñozJaramillo, A., DasiEspuig, M., Balmaceda, L., & DeLuca, E. 2013a, ApJ, 767, L25 [Google Scholar]
 MuñozJaramillo, A., Balmaceda, L., & DeLuca, E. 2013b, Phys. Rev. Lett., 111, 041106 [NASA ADS] [CrossRef] [Google Scholar]
 Norton, A. A., & Gallagher, J. C. 2010, Sol. Phys., 261, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Ogurtsov, M., & Lindholm, M. 2011, ISRN Astronomy and Astrophysics, 2011, 640817 [CrossRef] [Google Scholar]
 Passos, D., Charbonneau, P., & Beaudoin, P. 2012, Sol. Phys., 279, 1 [Google Scholar]
 Piotrowski, Z. P., Smolarkiewicz, P. K., Malinowski, S. P., & Wyszogrodzki, A. A. 2009, J. Comput. Phys., 228, 6268 [NASA ADS] [CrossRef] [Google Scholar]
 Pipin, V. V., & Sokoloff, D. D. 2009, Phys. Scr., 84, 065903 [Google Scholar]
 Prusa, J. M., Smolarkiewicz, P. K., & Wyszogrodzki, A. A. 2008, Comp. Fluids, 37, 1193 [Google Scholar]
 Racine, E., Charbonneau, P., Ghizaru, M., Bouchat, A., & Smolarkiewicz, P. K. 2011, ApJ, 735, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Schatten, K. 2005, Geophys. Res. Lett., 32, 21106 [NASA ADS] [CrossRef] [Google Scholar]
 Smolarkiewicz, P. K., & Charbonneau, P. 2013, J. Comput. Phys., 236, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Smolarkiewicz, P. K., & Margolin, L. G. 2007, Studies in geophysics, in Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics, eds. F. F. Grinstein, L. G. Margolin, & W. Rider (Cambridge University Press), 413 [Google Scholar]
 Svalgaard, L., Cliver, E. W., & Kamide, Y. 2005, Geophys. Res. Lett., 32, 1104 [NASA ADS] [CrossRef] [Google Scholar]
 Temmer, M., Rybak, J., Bendik, P., et al. 2006, A&A, 447, 735 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin, I. G. 2013, Liv. Rev. Sol. Phys., 10, 1 [Google Scholar]
 Usoskin, I. G., Hulot, G., Gallet, Y., et al. 2014, IPGP AAL Rev., in press [Google Scholar]
 Weiss, N. O. 2010, Astron. Geophys., 51, 3.09 [CrossRef] [Google Scholar]
All Tables
Comparison between the Pearson’s linear correlation coefficients, r, obtained between several quantities measured at cycles N and N + 1 for our simulated B_{F} (39 cycles) and the observed SSN (23 cycles).
Comparison between the correlation coefficients obtained between several quantities for the millennium simulation toroidal proxies B_{N} and B_{S}, and observational SSA data for the north and south.
Correlations between the two simulation proxies, B vs. A for both hemispheres and for the observational data, sunspot areas vs. polar flux timeseries.
All Figures
Fig. 1 Panel A) Time latitude diagram of the zonallyaveraged toroidal field taken at r/R_{⊙} = 0.718, which marks the interface between the convecting layers and underlying stably stratified region. Panel B) Zonallyaveraged radial field extracted from the simulation’s subsurface layer, at r/R_{⊙} = 0.938. The color scale encodes the magnetic field strength, values in Tesla. Panel C) Evolution of the total magnetic energy integrated over the whole computational domain. The simulation uses a hydrostatic initial state, subjected to small perturbations in fluid velocity and magnetic field. Convection reaches a statistically stationary state after about a year, but the buildup of a stable, cyclic largescale axisymmetric magnetic component sets in after about fifty years, here first in the southern hemisphere. Dipolar parity is maintained here through 20 full magnetic cycles. 

In the text 
Fig. 2 Meridional plot of the large scale, zonallyaveraged toroidal field (panel A)) and zonallyaveraged radial field (panel B)), extracted at t = 534 yr, a cycle maximum indicated by the vertical dotted black line in Fig. 1. The color scale on the right is in Tesla. The bounding boxes represent the integration regions for the proxies computation, individually for the northern and southern hemispheres. The dashed circular arc corresponds to r/R_{⊙} = 0.718, the interface between the convecting layers and underlying convectively stable fluid. 

In the text 
Fig. 3 Time series of the fulldisk toroidal field proxy B_{F}. The vertical dashed line represents the beginning of the sequence used for the foregoing analysis, which spans 39 sunspotlike cycles. Cycles are numbered from one minimum to the next, in keeping with the numbering convention of sunspot cycles. 

In the text 
Fig. 4 Histograms of frequency counts for the cycle periods for our fulldisk toroidal proxy B_{F} (panel A)) and the smoothed monthly SSN (panel B)). The vertical dashed lines represent the distributions means at 40.5 and 11 years correspondingly. Although the means differ by a factor of almost four, the two distributions show similar standard deviations and skewness (see text). 

In the text 
Fig. 5 Panels A) and B) Power spectrum for B_{F} and for SSN respectively. The dashed vertical lines mark important frequencies in the spectra and are defined, from left to right as ω_{1}, ω_{2}, ω_{3}, ω_{4} and ω_{5} (the latter is not highlighted for the SSN). 

In the text 
Fig. 6 Histograms of frequency counts for the cycle peak amplitude for A) our fulldisk toroidal proxy B_{F} and B) the SSN. The vertical dashed lines represent the distributions means. 

In the text 
Fig. 7 The GnevyshevOhl pattern in our fulldisk toroidal proxy. Panel A) shows the sequence of cycle maxima (black), and the corresponding 121 running mean thereof (dashed); subtracting the latter yields the detrended sequence plotted in panel B). The GO pattern shows up here as a sustained regular alternance between positive and negative detrended values, e.g. between simulated cycles 16 and 30 (thick segments). Panel C) shows a similarly detrended cycle amplitude sequence for the monthly smoothed SSN; here the GO pattern holds uninterrupted from sunspot cycle 9 to 21 (thick segments). 

In the text 
Fig. 8 Solar cycle proxies built from the millennium simulation. Black and red lines represent magnetic activity on the northern and southern hemispheres correspondingly. Panel A) shows the toroidal proxies, B_{N} and −B_{S}. In panel B) we present the poloidal (dipolar) proxies, A_{N} and −A_{S}. The southern proxies are plotted with a negative sign to better illustrate the coupling differences between both hemispheres. 

In the text 
Fig. 9 Waldmeiertype relationships for the magnetic proxies built from the millennium simulation. A) Period vs. Maxima; B) Mean Growth Rate vs. Maxima; C) Period vs. Rise Time; and D) Rise Time vs. Maxima. The similar anticorrelation observed in Period vs. Maxima (in A)) and Rise Time vs. Maxima (in D)) is a consequence of the good correlation between Period and Rise Time (in C)). 

In the text 
Fig. 10 Histograms of the cycle peak amplitudes, based on the toroidal proxies, measured in the A) north and B) south hemispheres. Although the distributions are built from a relatively small sample of cycle maxima, both hemispheres show a clear hint of bimodality, in the form of a subgroup of lowamplitude cycles. Vertical dashed lines indicate mean values 

In the text 
Fig. 11 Correlation between northern and southern maxima A) and period B) in the time series of toroidal proxies B_{N} and B_{S}. In the two panels the thick blue cross indicates the mean values ±one standard deviation. In A) the horizontal and vertical dashed lines, drawn at B_{N} = 0.6 and B_{S} = 0.6, delineate the cycles belonging to the two modes of the distributions in Fig. 10. The diagonal light gray line indicates a perfect correlation (B_{N} = B_{S}). 

In the text 
Fig. 12 Histograms for the periods measured in the A) north and B) south hemispheres. Signs of bimodality are apparent again here, although not as welldefined in the northern hemisphere as in the south. However, these subgroups of longperiod cycles do not map onto the subgroups of lowamplitude cycles in Fig. 10 (see text). 

In the text 
Fig. 13 The power spectrum for A) B_{N} and for B) B_{S}. The dashed vertical lines mark important frequencies in the spectra and are defined, from left to right as ω_{1}, ω_{2}, ω_{3}, ω_{4} and ω_{5}. 

In the text 
Fig. 14 Time lag measured between the two hemispheres for minima (blue, dashed line) and maxima (black, thick line). Negative values indicate that the northern hemisphere cycle leads its southern counterpart. The numbers inside the squares represent the simulated cycle number. The mean lag values are −2.6 yr for minima and −2.3 yr for maxima, amounting to about 6% of the mean cycle duration. The largest individual lag values approach 10 years for cycles 3, 20 and 33, i.e. 25% of the mean cycle duration. 

In the text 
Fig. A.1 Color coding used to show different values in the correlation significance level P. Significance below 70% has a white background. 

In the text 
Fig. A.2 Complete correlation table for the full disk proxy B_{F}, including correlations with a full disk poloidal proxy, A_{F}, which were not mentioned in the main text or included in Table 1. 

In the text 
Fig. A.3 Complete correlation table between cycle measures in the northern hemisphere (cf. Table 2 and Sect. 5.4). 

In the text 
Fig. A.4 Complete correlation table between cycle measures in the southern hemisphere (cf. Table 2 and Sect. 5.4). 

In the text 
Fig. A.5 Complete correlation table between cycle measures in the northern vs. southern hemispheres. 

In the text 
Fig. A.6 Complete correlation table between lags and cycle measures in the northern and southern hemispheres. 

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.