© 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, well-documented 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 so-called 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 pre-telescopic 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 thermally-driven convection, it is perhaps not surprising to observe cycle-to-cycle 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 mean-field and mean-field-like 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 self-consistently 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 large-scale magnetic fields undergoing more or less regular polarity reversals, some even producing sustained solar-like 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 convectively-driven 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 simulation-based proxies built to allow comparison with sunspot data. In Sect. 4 we examine the statistical properties of our pseudo-sunspot cycles, and the degree to which they exhibit known patterns in sunspot cycle data, namely the so-called Waldmeier and Gnevyshev-Ohl 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 EULAG-MHD code, the MHD generalization described in Smolarkiewicz & Charbonneau (2013) of the robust multi-scale flow solver EULAG (Prussa et al. 2008).

EULAG-MHD belongs to the class of so-called Implicit Large-Eddy 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 multi-scale 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 half-period 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⁷ 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 large-scale, axisymmetric magnetic component building up in the simulation have taken place at a regular cadence. This amounts to 40 sunspot-like cycles, which present variations in amplitude and periodicity but remained generally well-synchronized across both hemispheres. We thus define the large-scale magnetic field to be associated with our simulated solar cycles as the zonal average of the total magnetic field.

Figure 1A shows time-latitude diagrams of the toroidal (zonally-oriented) large-scale 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 mid-latitudes 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 large-scale 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 time-latitude diagram of the near surface radial magnetic component (here at 0.94 r/R_⊙). It is characterized by a well-defined 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 zonally-averaged toroidal field taken at r/R⊙ = 0.718, which marks the interface between the convecting layers and underlying stably stratified region. Panel B) Zonally-averaged 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 large-scale 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.

Open with DEXTER

Fig. 2 Meridional plot of the large scale, zonally-averaged toroidal field (panel A)) and zonally-averaged 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.

Open with DEXTER

Figure 1C shows a time series of total magnetic energy, integrated over the full simulation domain. The cycle of the large-scale axisymmetric magnetic field shows up quite prominently as a regular oscillation with well-defined period and amplitude, both showing relatively modest cycle-to-cycle fluctuations.

This simulation presents other features (not shown here) that are reasonably solar-like and similar to those found in previous data sets produced with EULAG-MHD, 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 full-disk toroidal field proxy BF. The vertical dashed line represents the beginning of the sequence used for the foregoing analysis, which spans 39 sunspot-like cycles. Cycles are numbered from one minimum to the next, in keeping with the numbering convention of sunspot cycles.
Open with DEXTER

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 large-scale, (zonally averaged) axisymmetric toroidal field ⟨ B_φ ⟩ over a thin layer immediately beneath the base of the convective layers, where the large-scale 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₁ ≃ 0.67R_⊙, r₂ ≃ 0.72R_⊙, θ₁ ≃ 33°, θ₂ ≃ 73°, and r₃ ≃ 0.91R_⊙, r₄ ≃ 0.95R_⊙, θ₃ ≃ 62°, θ₄ ≃ 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ñoz-Jaramillo 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 semi-automatically 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 one-month temporal cadence as our computational proxies. We retrieved the SSN data from the SIDC¹ 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 re-analyse 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 full-disk toroidal proxy BF (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).
Open with DEXTER

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.

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 long-term modulation of the cycle amplitude. In the case of the SSN, the best-known such long periodicity is the so-called 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 full-disk toroidal proxy B_F (A) and SSN (B).

	Fig. 5 Panels A) and B) Power spectrum for BF 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).
Open with DEXTER

	Fig. 6 Histograms of frequency counts for the cycle peak amplitude for A) our full-disk toroidal proxy BF and B) the SSN. The vertical dashed lines represent the distributions means.
Open with DEXTER

The primary peak, labelled ω₃, 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, ω₄ and ω₅, that correspond to the first and second harmonics of the main peak of B_F. For the SSN, ω₅ 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 ω₁ and ω₂ which correspond to 279 yr and 170 yr, for B_F and 89 yr and 53 yr for SSN, respectively. These peaks are not well-defined 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 well-established, 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 so-called Gnevyshev-Ohl rule, whereby even-numbered sunspot cycles tend to have lower amplitudes than their odd-numbered counterparts, which leads to an alternance between higher-than-average and lower-than-average 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 1-2-1 running mean: (3)Figure 7A illustrates this procedure, with the black solid and dashed lines showing the raw and 1-2-1 smoothed sequences of maxima, and panel B the result of the above detrending procedure. The Gnevyshev-Ohl 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 Gnevyshev-Ohl pattern sustained here from cycle 9 to 21.

Fig. 7 The Gnevyshev-Ohl pattern in our full-disk toroidal proxy. Panel A) shows the sequence of cycle maxima (black), and the corresponding 1-2-1 running mean thereof (dashed); subtracting the latter yields the detrended sequence plotted in panel B). The G-O 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 G-O pattern holds uninterrupted from sunspot cycle 9 to 21 (thick segments).

Open with DEXTER

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ñoz-Jaramillo et al. 2013a), in polar fields (Svalgaard et al. 2005; Schatten 2005; Muñoz-Jaramillo 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 cross-hemispheric 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ñoz-Jaramillo 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ñoz-Jaramillo 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, BN and −BS. In panel B) we present the poloidal (dipolar) proxies, AN and −AS. The southern proxies are plotted with a negative sign to better illustrate the coupling differences between both hemispheres.

Open with DEXTER

For each hemisphere, we segment the proxy time series into cycles, following the same strategy and numbering convention introduced earlier for the full-disk 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.

As in the case of the analysis carried out for the full-disk proxies (viz. Table 1), we observe both similarities and differences in the levels of correlation obtained for the simulation-based 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, color-coded 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 full-disk 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 short-term prediction of cycle peak amplitude, also present in our full-disk 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 Waldmeier-type 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)).
Open with DEXTER

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 large-scale 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 low-amplitude cycles well separated from the group of high-amplitude 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 long-term solar activity from the abundances of cosmogenic radioisotopes also show a similar bimodal distribution, characterized by a secondary population of low-activity 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 low-amplitude cycles. Vertical dashed lines indicate mean values
Open with DEXTER

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 cycle-to-cycle 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 low-amplitude 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 low-amplitude 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 cross-hemispheric influence. It also explains why the bimodality so clear in Fig. 10 – again, small-N 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 BN and BS. 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 BN = 0.6 and BS = 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 (BN = BS).

Open with DEXTER

Turning to proxy amplitude at cycle minima, no cross-hemispheric correlation is found (r = −0.07 with P = 32%). For the SSA time series the cross-hemispheric 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 long-period cycles separated from a well-defined main group by a gap around 42.5 yr (see Fig. 11B). We note that this subgroup of long-duration cycles does not map onto the low-amplitude group of cycles in the peak amplitude distribution; in fact, examination of Fig. 9A reveals that the low-amplitude 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 well-defined in the northern hemisphere as in the south. However, these subgroups of long-period cycles do not map onto the subgroups of low-amplitude cycles in Fig. 10 (see text).
Open with DEXTER

	Fig. 13 The power spectrum for A) BN and for B) BS. The dashed vertical lines mark important frequencies in the spectra and are defined, from left to right as ω1, ω2, ω3, ω4 and ω5.
Open with DEXTER

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.

Open with DEXTER

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 (ω₃, ω₄ and ω₅), are practically indistinguishable between both hemispheres, the lower frequencies show some noticeable differences. For B_N, ω₁ and ω₂ 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 solar-like long periodicities.

5.2. Lag analysis

The northern and southern hemisphere cycles in Fig. 8 are clearly well-synchronized, 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 large-scale 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 cross-hemispheric 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 39-member 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⁵ 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.

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 self-regulating 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. Gnevyshev-Ohl pattern

The Gnevyshev-Ohl pattern of alternating higher- and lower-than-average 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 G-O pattern, while a shorter, out-of-phase G-O pattern is sustained from cycle 5 to 9 in the north, with the consequence that the full-disk proxy only shows a brief G-O episode from cycle 1 to 4 (viz. Fig. 7). A similar situation occurs at the end of the simulation, leaving only a 5-cycle-long G-O pattern in the full disk proxy. On the other hand, the 15-cycle-long G-O 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 in-phase contribution from cycles 22 to 27. The G-O 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ñoz-Jaramillo 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ñoz-Jaramillo 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 low-amplitude 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 low-amplitude 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 Waldmeier-like 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 self-similarity 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 mean-field dynamo models (e.g. Karak & Choudhuri 2011 or Pipin & Sokoloff 2011). Our simulation achieves the same thing through dynamically-driven fluctuations.

The Gnevyshev-Ohl 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 tachocline-like layer, any amplitude bias introduced by a higher-than-average 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 lower-than-average cycle. Then, the next following cycle has less opposite-polarity magnetic flux to reverse, and so will end up with a higher-than-average amplitude once again. The presence of a Gnevyshev-Ohl-like 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 Gnevyshev-Ohl 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 mean-field 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 low-latitude 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 cross-hemispheric 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 cross-equatorial 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 higher-order) 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 large-scale flows developing in the simulation. Morphologically speaking, trajectories associated with low-amplitude cycles are similar to their higher amplitude counterparts.

Another possibility is that the low-amplitude 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 EULAG-MHD simulations such as the one studied here, undergo a transition to a non-axisymmetric large-scale dynamo mode, still characterized by regular polarity reversals on a similar multi-decadal 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 non-axisymmetric dynamo mode could then be the sought-after 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 EULAG-MHD simulation of solar convection developing a large-scale axisymmetric magnetic field component undergoing reasonably solar-like 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 observationally-inferred solar-like fluctuations patterns in amplitude and duration of our simulated cycles are recovered, including Waldmeier-law-like correlations between cycle parameters (rise time, period, amplitude), hints of long, Gleissberg-like periodicities in cycle amplitude, as well as the Gnevyshev-Ohl pattern of extended sequences of alternating higher- and lower-than-average successive cycle amplitudes.

Cycles in both hemispheres are well-coupled 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 low-amplitude cycle and hemispheric lag in previous, current or subsequent cycle.

It is noteworthy that solar-like 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ñoz-Jaramillo 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 large-scale 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?

Acknowledgments

We wish to acknowledge the contribution of Mihai Ghizaru in the initial design of the EULAG-MHD 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ñoz-Jaramillo for kindly providing us the SSA and PF observational data. Computing time for this project was granted through Calcul-Canada and simulations carried out on the computing platforms of Calcul-Québec. D. Passos acknowledges support from the Fundação para a Ciência e Tecnologia (FCT) grant SFRH/BPD/68409/2010, CENTRA-IST 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

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.
Open with DEXTER

	Fig. A.2 Complete correlation table for the full disk proxy BF, including correlations with a full disk poloidal proxy, AF, which were not mentioned in the main text or included in Table 1.
Open with DEXTER

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

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

	Fig. A.5 Complete correlation table between cycle measures in the northern vs. southern hemispheres.
Open with DEXTER

	Fig. A.6 Complete correlation table between lags and cycle measures in the northern and southern hemispheres.
Open with DEXTER

All Tables

All Figures

	Fig. 1 Panel A) Time latitude diagram of the zonally-averaged toroidal field taken at r/R⊙ = 0.718, which marks the interface between the convecting layers and underlying stably stratified region. Panel B) Zonally-averaged 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 large-scale 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.
Open with DEXTER
In the text

	Fig. 2 Meridional plot of the large scale, zonally-averaged toroidal field (panel A)) and zonally-averaged 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.
Open with DEXTER
In the text

	Fig. 3 Time series of the full-disk toroidal field proxy BF. The vertical dashed line represents the beginning of the sequence used for the foregoing analysis, which spans 39 sunspot-like cycles. Cycles are numbered from one minimum to the next, in keeping with the numbering convention of sunspot cycles.
Open with DEXTER
In the text

	Fig. 4 Histograms of frequency counts for the cycle periods for our full-disk toroidal proxy BF (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).
Open with DEXTER
In the text

	Fig. 5 Panels A) and B) Power spectrum for BF 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).
Open with DEXTER
In the text

	Fig. 6 Histograms of frequency counts for the cycle peak amplitude for A) our full-disk toroidal proxy BF and B) the SSN. The vertical dashed lines represent the distributions means.
Open with DEXTER
In the text

	Fig. 7 The Gnevyshev-Ohl pattern in our full-disk toroidal proxy. Panel A) shows the sequence of cycle maxima (black), and the corresponding 1-2-1 running mean thereof (dashed); subtracting the latter yields the detrended sequence plotted in panel B). The G-O 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 G-O pattern holds uninterrupted from sunspot cycle 9 to 21 (thick segments).
Open with DEXTER
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, BN and −BS. In panel B) we present the poloidal (dipolar) proxies, AN and −AS. The southern proxies are plotted with a negative sign to better illustrate the coupling differences between both hemispheres.
Open with DEXTER
In the text

	Fig. 9 Waldmeier-type 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)).
Open with DEXTER
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 low-amplitude cycles. Vertical dashed lines indicate mean values
Open with DEXTER
In the text

	Fig. 11 Correlation between northern and southern maxima A) and period B) in the time series of toroidal proxies BN and BS. 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 BN = 0.6 and BS = 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 (BN = BS).
Open with DEXTER
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 well-defined in the northern hemisphere as in the south. However, these subgroups of long-period cycles do not map onto the subgroups of low-amplitude cycles in Fig. 10 (see text).
Open with DEXTER
In the text

	Fig. 13 The power spectrum for A) BN and for B) BS. The dashed vertical lines mark important frequencies in the spectra and are defined, from left to right as ω1, ω2, ω3, ω4 and ω5.
Open with DEXTER
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.
Open with DEXTER
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.
Open with DEXTER
In the text

	Fig. A.2 Complete correlation table for the full disk proxy BF, including correlations with a full disk poloidal proxy, AF, which were not mentioned in the main text or included in Table 1.
Open with DEXTER
In the text

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

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

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

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