Long-term modulation of solar cycles

TL;DR

This paper reviews the long-term variability of solar activity over millennia, discussing grand minima and maxima, quasi-periodicities, and the stochastic nature of solar cycles, supported by cosmogenic isotope data and dynamo models.

Contribution

It provides an up-to-date overview of long-term solar activity variability, integrating observational constraints with theoretical dynamo models and emphasizing the stochastic components of solar cycles.

Findings

Solar cycles exhibit significant modulation in length, shape, and magnitude.

Long-term solar activity shows quasi-periodicities like Gleissberg, Suess/de Vries, and Hallstatt cycles.

Solar cycles contain a random component with no long-term phase locking.

Abstract

Solar activity has a cyclic nature with the ~11-year Schwabe cycle dominating its variability on the interannual timescale. However, solar cycles are significantly modulated in length, shape and magnitude, from near-spotless grand minima to very active grand maxima. The ~400-year-long direct sunspot-number series is inhomogeneous in quality and too short to study robust parameters of long-term solar variability. The cosmogenic-isotope proxy extends the timescale to twelve millennia and provides crucial observational constraints of the long-term solar dynamo modulation. Here, we present a brief up-to-date overview of the long-term variability of solar activity at centennial--millennial timescales. The occurrence of grand minima and maxima is discussed as well as the existing quasi-periodicities such as centennial Gleissberg, 210-year Suess/de Vries and 2400-year Hallstatt cycles. It is shown that the solar cycles contain an important random component and have no clock-like phase locking implying a lack of long-term memory. A brief yet comprehensive review of the theoretical perspectives to explain the observed features in the framework of the dynamo models is presented, including the nonlinearity and stochastic fluctuations in the dynamo. We keep gaining knowledge of the processes driving solar variability with the new data acquainted and new models developed.

Equations

$$R_{\mathrm{W}}=k\cdot (10\cdot G+S),$$

$$S=\frac{E(\sigma {{}_{\varphi }}^2)}{E(\sigma {{}_P}^2)}$$

$$S_c=\frac{E(\sigma {{}_{\varphi }}^2)}{E(\sigma {{}_P}^2)}=\frac{N^2}{2(N+1{)}^2}$$

$$S_r=\frac{E(\sigma {{}_{\varphi }}^2)}{E(\sigma {{}_P}^2)}=\frac{N(N+3)}{15(N+1)}$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{v}\times \mathbf{B}-\eta \nabla \times \mathbf{B}),$$

$$\rho\left[\frac{\partial\boldsymbol{v}}{\partial t}+(\boldsymbol{v}\cdot\mbox{\boldmath$\nabla$})\boldsymbol{v}\right]=-\mbox{\boldmath$\nabla$}P+\boldsymbol{J}\times\boldsymbol{B}+\mbox{\boldmath$\nabla$}\cdot(2\nu\rho S)+\boldsymbol{F},$$

$$\frac{\partial A}{\partial t}+\frac{1}{s}(\boldsymbol{v_{m}}\cdot\mbox{\boldmath$\nabla$})(sA)=\eta_{t}\left(\nabla^{2}-\frac{1}{s^{2}}\right)A+\alpha B,$$

$$\frac{\partial B}{\partial t}+\frac{1}{r}\left[\frac{\partial(rv_{r}B)}{\partial r}+\frac{\partial(v_{\theta}B)}{\partial\theta}\right]=\eta_{t}\left(\nabla^{2}-\frac{1}{s^{2}}\right)B+s(\boldsymbol{B_{p}}\cdot\mbox{\boldmath$\nabla$}){\rm\Omega}+\frac{1}{r}\frac{d\eta_{t}}{dr}\frac{\partial(rB)}{\partial r},$$

Full text

1 Introduction

Sun is a magnetically active star whose activity is a result of the magnetic dynamo process operating in the Sun’s convection zone (see, e.g., Karak et al, 2014; Charbonneau, 2020). Solar surface magnetic activity varies cyclicly with the main period of about 11 years (called the Schwabe cycle) or, considering inversion of the sign of its magnetic polarity, the 22-year Hale cycle. More details can be found in an extensive review by Hathaway (2015). The physics of the dynamo mechanism is currently believed to be reasonably well understood. However, solar cyclicity is far from being a regularly ticking clock and experiences essential long-term variability at timescales longer than the Schwabe cycle. The solar cycles are not perfectly regular and vary in length, shape, and strength/intensity, or even can enter periods of almost inactive state, called grand minima of solar activity (e.g., Usoskin, 2017).

The standard index quantifying solar activity is related to sunspot numbers which are available from 1610 AD onward with the quality degrading backwards in time, as discussed in Section 2. On one hand, this 410-year-long series exhibits a great deal of variability covering the range from an almost spotless period of the Maunder minimum between 1645 – 1715 AD (Eddy, 1976) to an epoch of very active Sun between 1940 – 2009 called the Modern grand maximum (Solanki et al, 2004; Usoskin et al, 2007). This great variability raises important questions, answers to which can put crucial observational constraints on the solar/stellar dynamo theory:

• •

Do the changes between the Maunder minimum and the Modern grand maximum cover the full possible range of solar variability?

• •

Do the grand minima and maxima represent special states of the solar dynamo or simply represent the tails of the distribution?

• •

How typical are these changes?

• •

Do the grand minima episodes appear periodically or randomly?

• •

What physical processes drive such changes?

The four-century-long sunspot number series is not sufficiently long to answer these questions, and a much longer dataset is needed to form a basis for the answers. Fortunately, solar activity can be reliably reconstructed from indirect natural proxy data (cosmogenic radioisotopes) on the timescale of 10 – 12 millennia, during the period of the Holocene with a stable warm climate on Earth, as discussed in Section 3. This reconstruction extends the solar-activity dataset by a factor of about 25 making it possible to perform a thorough statistical analysis of solar variability as discussed in Section 4, while statistical properties of the solar-cycle modulation are summarized in Section 5. In Section 6, we discuss the implications of the long-term solar variability for the solar dynamo theory and our present level of understanding of the related physics.

2 Direct Sunspot number series since 1610

Sunspots have been more or less systematically studied since 1610, soon after the invention of the telescope. Thousands of observational records and drawings exist in archives as being continuously recovered and analyzed (e.g., Vaquero and Vázquez, 2009; Arlt and Vaquero, 2020). The most recent and continuously updated database of raw sunspot-group observation is collected at the HASO (Historical Archive of Sunspot Observations, http://haso.unex.es/haso – Vaquero et al, 2016).

Despite numerous observational records, it was noticed only in the middle of the 18th century by the Danish astronomer Christian Horrebow and finally confirmed in the early 19th century by the German observer Heinrich Schwabe, that the number of sunspots varies cyclicly with about 10-year period. This cycle was later shown to be of about 11 years mean length and appears to be a fundamental feature of solar activity and is now called the Schwabe cycle. More details of the sunspot number measurements and reconstructions can be found elsewhere in this volume or in comprehensive reviews by Hathaway (2015) and Usoskin (2017).

2.1 Wolf sunspot series $R_{\rm W}$ and International sunspot number $R_{\rm I}$

Following the discovery of the solar cycle, Rudolf Wolf from Zürich Observatory founded a synthetic index called the sunspot number presently known as Wolf or Zürich sunspot number $R_{\rm W}$ (WSN) defined as

[TABLE]

where $G$ and $S$ are the numbers of sunspot groups and all sunspots, including those in groups, respectively, visible on the solar disc during a given day by the primary observer whose quality scaling factor $k$ is set to reduce his/her counts to the reference observer with $k\equiv$1. Obviously, the sunspot number is not the same as the number of spots, and for a single sunspot, $R_{\rm W}$=11 assuming $k$=1. This series, constructed by R. Wolf in 1861 using his own and recovered earlier observations, formally covered the period since 1749 (solar cycle SC #1 in Wolf’s numbering), but in fact, it was more or less reliable only since the 1820’s when H. Schwabe started his observations. Later it was extended back to 1700 with unreliable data. The compilation of the $R_{\rm W}$ was continued at Zürich by Wolf’s successors Wolfer, Brunner, Waldmeier and Koeckelenbergh until 1981 when the formation of the sunspot series was transferred to the Royal Observatory of Belgium (Clette et al, 2007).

Until 1981, the $R_{\rm W}$ was constructed considering the observation of only one primary observer for each day, all other observations were discarded. This series could not, till now, be revisited or redone because of the lack of original raw data. Accordingly, when several apparent inhomogeneities were found in the standard Wolf sunspot series (Leussu et al, 2013; Clette et al, 2014; Lockwood et al, 2014), only step-wise corrections to the old series could be done (Clette et al, 2014; Clette and Lefèvre, 2016). This ‘corrected’ sunspot series is known as the International sunspot number series version 2.0, $R_{\rm I}$(2.0), and is available at the SILSO (Sunspot Index and Long-term Solar Observations, https://www.sidc.be/silso/datafiles) formally since 1700. The $R_{\rm I}$(2.0) is shown in Figure 1a along with the standard Zürich sunspot cycle numbering.

Although the update of the series was through several adjustments of scaling jumps, an important effort is currently done by the community to restore and digitize old raw data (Clette et al, 2021) so that it will be possible to redo the sunspot number series from scratch increasing its reliability and assessing realistic uncertainties.

2.2 Group sunspot number series $G_{\rm N}$

Since the sunspot number (Equation 1) includes both numbers of sunspot groups (weighted by a factor of 10) and individual sunspots, it is sensitive to the quality of observations. This was addressed by Hoyt and Schatten (1998) who noticed that sunspot groups are defined more reliably than individual spots and created the group sunspot number series $G_{\rm N}$ which is simply the number of sunspot groups $G$ on the solar disc corrected for the observer’s quality. This series is shown in Figure 1b. Sometimes it is scaled up to match the values typical for $R_{\rm W}$. However, contrary to $R_{\rm W}$, $G_{\rm N}$ is based on the average of all available observations for each day, not only the primary ones. Another principal difference between $R_{\rm W}$ and $G_{\rm N}$ is that Hoyt and Schatten (1998) created and published a full database of raw data they used to construct the $G_{\rm N}$ series. Accordingly, this series can be completely redone as a whole, without limitation to the ‘correction factors’.

It was recognized that the original $G_{\rm N}$ underestimated solar activity during the 19-th century (Clette et al, 2014), and several efforts have been made to revisit it using different methodologies and inter-calibrations (e.g., Svalgaard and Schatten, 2016; Usoskin et al, 2016b; Chatzistergos et al, 2017; Willamo et al, 2017). One of the reconstructions is also shown in Figure 1b. However, these new series often moderately disagree with each other illustrating the problem of compiling a homogeneous series from individual raw datasets (Muñoz-Jaramillo and Vaquero, 2019). It is presently impossible to decide between different reconstructions of the group sunspot series, but the zoo of those gives a clue of what the related uncertainties are, and presently they are bounded by the series of Svalgaard and Schatten (2016) from the top and from below by Hoyt and Schatten (1998).

3 Cosmogenic-isotope-based reconstructions of long-term solar variability

The sunspot number series covers ca. 410 years in the past with the quality degrading back in time (Muñoz-Jaramillo and Vaquero, 2019) and principally cannot be extended before the 17-th century because of the lack of instrumental data. Unaided (naked-eye) observations of sunspots do not provide systematic quantitative information on solar activity (Usoskin, 2017). There are some other proxy-based indices of solar activity, such as geomagnetic or heliospheric activity, and radio-emission of the Sun, but they all are based on scientific measurements and typically do not go beyond the middle of the 19-th century. Fortunately, there is one solar-activity proxy which can help in reconstructing solar variability on the multi-millennial timescale. This is related to cosmogenic radioisotopes which are produced and preserved in dateable archives in a natural way.

3.1 Method of cosmogenic isotopes

Solar surface magnetic activity and hot corona create the solar wind which is a supersonic outflow of solar coronal plasma permanently emitted from the Sun (see, e.g., Vidotto, 2021). Because of its high conductivity, solar wind drags away the solar magnetic field which appears ‘frozen’ in the solar-wind plasma. This wind radially expands forming the heliosphere, a region of about 200 astronomical units across which is totally controlled (in the magnetohydrodynamical sense) by the solar wind and magnetic field (e.g., Owens and Forsyth, 2013). The heliosphere makes an obstacle for charged highly energetic particles of galactic cosmic rays (GCRs) which permanently bombard it isotropically with nearly constant flux. Inside the heliosphere, cosmic rays are affected by four major processes, viz. scattering and diffusion on magnetic irregularities, convection by expanding solar wind, adiabatic cooling, and large-scale drifts. All these processes are ultimately driven by solar activity leading to the solar modulation of cosmic-ray flux near Earth so that the cosmic-ray flux is stronger when solar activity is weak and vice-versa (e.g., Potgieter, 2013). Thus, knowing the modulated flux of GCRs at a moment in time, one can assess the level of solar activity slightly before that (within one year – Koldobskiy et al, 2022). Of course, there were no scientific cosmic-ray detectors in the distant past, but there is a natural cosmic-ray monitor – cosmogenic radioisotopes.

Cosmogenic radioisotopes are unstable nuclides, which cannot survive from the time of the solar-system formation, and whose main source is related to nuclear reactions caused by cosmic rays in the Earth’s atmosphere (Beer et al, 2012). After production in the atmosphere by GCR, nuclides can be stored in natural independently dateable archives, such as tree trunks, polar ice cores, lake/marine sediments, etc. Accordingly, the flux of GCR can be estimated in the past by measuring the abundance of such isotopes in the archives, forming the only quantitative proxy of solar activity over long timescales (see more details in Beer, 2000; Usoskin, 2017). The most important cosmogenic isotopes are 14C ‘radiocarbon’ (half-life 5730 years) measured in dendrochronologically dated tree rings and 10Be ($\approx 1.4\cdot 10^{6}$ years) measured in glaciologically dated ice cores.

Conversion between the measured isotope concentration and production by cosmic rays requires a knowledge of the isotope’s transport and deposition processes which are currently well modelled (e.g., Roth and Joos, 2013; Heikkilä et al, 2013; Golubenko et al, 2021). Additionally, it needs to be corrected for the changing geomagnetic field (e.g., Pavón-Carrasco et al, 2018), and the resulting variability can be attributed to solar activity. The conversion from the cosmic-ray modulation to the heliospheric properties (open solar flux) and then to the pseudo-sunspot numbers is done via a chain of physics-based models making it possible to reconstruct solar activity and the related uncertainties (see, e.g., Usoskin, 2017; Wu et al, 2018).

3.2 Holocene ($\approx$12 kyr) decadal reconstruction

While the idea of the use of cosmogenic-isotope data as a proxy to solar activity has been discussed since long (Stuiver, 1961; Lal and Peters, 1962), first approaches were empirical as based on timescale separation of the cosmogenic data: timescales longer than 500 years were thought to be caused by changes in the large-scale geomagnetic field, while shorter time scales – by solar activity (Damon and Sonett, 1991). That approach made it possible to identify grand solar minima (Eddy, 1976; Stuiver and Braziunas, 1989) but was unable to provide a quantitative reconstruction of solar activity because both factors are important at the centennial timescales. A full reconstruction of solar activity from cosmogenic-isotope data became possible only after the development of models of cosmic-ray-induced atmospheric cascades (Masarik and Beer, 1999). The first quantitative reconstruction of solar activity using a physics-based approach was made by Usoskin et al (2003) on the millennial time scale (see also Solanki et al, 2004). Later the reconstructions were extended to the Holocene (the present period of stable warm climate lasting for about 12 millennia) using different cosmogenic isotopes (e.g., Vonmoos et al, 2006; Steinhilber et al, 2012; Usoskin et al, 2016b). The most recent and accurate multi-millennial solar-activity reconstruction by Wu et al (2018) is based on a multi-proxy Bayesian approach providing also realistic uncertainties. It is shown in Figure 2. One can see that solar activity varies essentially between the grand minima, visible at sharp dips down to 10 – 20 (in sunspot number, SN), and grand maxima when SN exceeds 60, while most of the time the solar-activity level remains moderate at SN$=40\pm 10$ (see more detail in Usoskin et al, 2014). The results of an analysis of the solar-activity variability are reviewed in Section 4.

Because of the low time resolution of the cosmogenic-isotope throughout the Holocene (typically decadal – see, e.g., Reimer et al, 2020), reconstructions of solar activity are also usually limited to the 10-year resolution being thus unable to resolve individual solar cycles. Long-term reconstructions of solar activity are limited to the Holocene timescale because of the stable climate so that the standard models of the isotope atmospheric transport and deposition can apply. However, for the ice-age-type of climate, the properties of the atmospheric transport are quite uncertain including the large-scale atmospheric and ocean circulation, which prevents quantitative assessment of solar activity. At present, there is no model which is able to handle this in a satisfactory manner, but progress is expected in the future.

3.3 $\approx$100 solar cycles reconstructed

Thanks to the recent technological progress, high-precision measurements of annual 14C concentrations have been performed with the annual resolution for the last millennium (Brehm et al, 2021). It allowed us to make, by applying the physics-based model, the first reliable reconstruction of individual solar cycles beyond the epoch of telescopic observations (Usoskin et al, 2021) as shown in Figure 3. Four known grand minima are seen – Oort, Wolf, Spörer and Maunder minima, and between the minima, there are clear solar cycles of variable amplitude. In this way, 85 individual solar cycles have been reconstructed from 14C of which 35 cycles are reasonably and well resolved, 21 are poorly and 29 are not reliably resolved, mostly during the grand minima of activity. Overall, including both direct solar observations and proxy-based reconstructions, we now have information on 96 solar cycles of which 50 are well resolved, thus nearly tripling the extent of the solar-cycle knowledge and doubling the number of well-defined cycles.

The extended statistic made it possible to perform a primary analysis of the solar-cycle parameters. The length of the well-defined cycles was $10.8\pm 1.4$ years which is in good agreement with $11.0\pm 1.1$ years known for the ISN dataset. The statistical significance of the Waldmeier rule (solar-cycle height is inversely correlated with the length of the ascending phase – high cycles rise fast) has been confirmed with the extended dataset, implying its robust nature (Usoskin et al, 2021). However, the Gnevyshev-Ohl rule of even–odd cycle pairing (Gnevyshev and Ohl, 1948; Usoskin et al, 2001) has not been confirmed, nor rejected with the extended data. A more detailed analysis of this new dataset is still pending.

4 Long-term solar activity

With the reconstructed long series, one can investigate properties of solar variability which pose observational constraints crucially important for solar physics but cannot be set by the too short-ranging conventional direct telescopic observations of the Sun. While the 11-year solar cycle forms the main feature of solar activity, the cycles are far from being perfect clock ticks – they vary by both duration and amplitude including periods of greatly suppressed activity, grand minima (see Figure 3). Here we review the most important features of long-term solar variability.

4.1 Long quasi-periodic variations (Gleissberg, Suess/de Vries, Hallstatt cycles)

It is hardly possible to distinguish whether solar variability on a long-term scale (Figure 2) is stochastic/chaotic or (quasi)periodic. Power-spectrum analyses are controversial but generally agree that there are three period ranges with apparent and barely significant variability. An example of the global wavelet power spectrum is shown in Figure 4.

One is the centennial variability, called the Gleissberg cycle, which is not a strict periodicity but a characteristic period range between 60 – 140 years (e.g., Peristykh and Damon, 2003; Ogurtsov, 2004). The Gleissberg cycle is clearly seen in the direct sunspot data but is less pronounced throughout the Holocene.

Another important periodicity is the Suess cycle (called also de Vries cycle in the literature), which has a narrow period range between 200 – 210 years and an intermittent occurrence. It is typically seen as a recurrence of grand minima within clusters of reduced solar activity (Usoskin et al, 2014) as seen, e.g., in Figure 3, but is not readily observed during the epochs of moderate solar activity.

Sometimes, the so-called Eddy millennial cycle is claimed to exist (Steinhilber et al, 2012), but it is unstable and cannot be identified in a significant way (see Figure 4).

Additionally, there exists a very-long cycle with a timescale of 2000 – 2400 years called the Hallstatt cycle (Damon and Sonett, 1991; Vasiliev and Dergachev, 2002; Usoskin et al, 2016a). Because of its length, it cannot be robustly defined in the $\approx$10-kyr time series (see Figure 4). The nature of the Hallstatt cycle is still unclear: it is likely to be ascribed to the Sun (Usoskin et al, 2016a) but geomagnetic or climatic origin cannot be excluded. Longer-scale variability cannot be reliably assessed from the cosmogenic-isotope data, in particular, because of the unresolved discrepancy between 14C and 10Be datasets on the multi-millennial timescale as probably related to the effect of deglaciation (e.g., Vonmoos et al, 2006; Usoskin et al, 2016a; Wu et al, 2018).

4.2 Grand minima and maxima

As seen, e.g., in Figures 2 and 3, solar activity sometimes drops fast, within one–two solar cycles, to the very quiet level with almost no sunspots on the solar surface. These drops are called grand minima of activity. Until the 1970s, the existence of such minima was debated, but Eddy (1976) had convincingly proved that the sunspot activity indeed dropped to almost no sunspots between 1645 – 1715 as confirmed also by other proxies such as auroral displays at mid-latitudes. That grand minimum was called the Maunder minimum. More grand minima have been found later using the cosmogenic-isotope data (e.g., Usoskin et al, 2007; Inceoglu et al, 2015). At present, about 30 grand minima of duration ranging between 40 – 70 (Maunder-type minima) and 100 – 140 years (Spörer-type) each, have been identified during the Holocene occupying about 1/6 of the time. It has been shown that the grand minima correspond to a special state of the solar dynamo (e.g., Usoskin et al, 2014).

Solar activity was abnormally high in the second half of the 20th century compared to the 19th or 21st centuries (Lockwood et al, 1999) but it was unknown whether this high level is unique or typical. Using the cosmogenic-isotope data, it was discovered that the period from the 1940s to 2010 was not unique and there are other similarly high but very rare episodes, that forms the concept of a grand solar maximum (Usoskin et al, 2003; Solanki et al, 2004). Grand maxima represent periods of enhanced solar activity covering at least a few solar cycles. There were about 20 grand maxima over the Holocene which cover $\approx$10 % of the time (Usoskin et al, 2007; Inceoglu et al, 2015), but they are defined not as robustly as grand minima. No apparent clustering in the grand-maxima occurrence or duration has been found, nor do they form a special distribution of solar cycles (Usoskin et al, 2014, 2016a). It is still unknown whether grand maxima make a special mode of the dynamo, similar to grand minima, or just represent a rare tail of the solar-cycle-strength distribution.

5 Statistical properties of the long-term modulation of solar cycles

As historical records show, solar cycles are highly variable in amplitude and length. The validity of theoretical models that attempt to predict this variability depends heavily on whether the cycle exhibits long-term phase stability or whether the phase is subject to a random walk, or a mixture of these. In the first of the two extreme cases, the system has infinite phase memory and in the second case no phase memory at all. Phase stability could be achieved through synchronization processes, such as high-quality torsional oscillations in the solar interior (Dicke, 1970) or the weak tidal forces of planets (e.g., Stefani et al, 2021). Dynamo models generally predict phase progression without memory. An insightful summary of the use of historical observations to explain solar phenomena was given by Vaquero and Vázquez (2009).

The question of the regularities and randomness of solar activity variability has been studied for a long time. For example, statistical methods including those based on the Lyapunov and Hurst exponents or Kolmogorov entropy (e.g., Ostriakov and Usoskin, 1990; Mundt et al, 1991; Carbonell et al, 1994; Ruzmaikin et al, 1994; Lepreti et al, 2021) were inconclusive, implying that a mixture of different components is likely (see more details in Usoskin, 2017; Petrovay, 2020).

Various publications (e.g., Lomb, 2013; Russell et al, 2019; Stefani et al, 2020) claim that the solar cycle is phase stable. However, to answer the question of whether the phase is stable or not, one needs a clear definition of phase stability, an appropriate statistical analysis as well as reliable data on which to apply the analysis. Dicke (1978) and Gough (1978) were among the first to perform a systematic statistical analysis based on telescopic sunspot records. Independently, but using similar concepts, they concluded that the time span of the available data was too small for a clear distinction between the two cases. Later, Gough (1981, 1983, 1988) corrected and modified his earlier analysis without altering the conclusion.

Interestingly, Eddington and Plakidis (1929) analyzed the light-curve variations of long-period variable stars, a problem close to the variability of the solar cycle. By deriving a statistical function to which the processed observational data were fitted, they were able to determine two indicators for the composition of clock-synchronised phase perturbations and random phase perturbations of the light signal.

Weisshaar et al (2023) have revisited Gough’s analysis based on newly available data. For clarity, a brief outline of Gough’s test is given here: From the arithmetic mean of the individual cycle lengths (Gough, 1981), the regular minima or maxima of the hypothetical dynamo or clock cycles and thus the corresponding phase deviations can be determined as the difference to the observed minima or maxima. The basic statistics are the expectation values of the variances of cycle period, $E(\sigma_{P}^{2})$, and phase, $E(\sigma_{\phi}^{2})$. The final statistics is defined as the ratio of the two variances to cancel out the unknown fluctuation amplitude:

[TABLE]

Later, Gough (1983) modified the method by replacing the arithmetic mean of the cycle period with a value that minimizes the variance of the phase deviations, resulting in a more sensitive distinction between the clock regime and the random phase regime. Calculating the expectation values of the variances for the two cases, one obtains the following expressions for $S_{c}$ (clock) and $S_{r}$ (random phase) using the modified method:

[TABLE]

which asymptotically reaches $N\to\infty$, ${S_{c}}\to{1\over 2}$;

[TABLE]

which asymptotically reaches $N\to\infty$, ${S_{r}}\to{N\over 15}$.

The procedure to apply Gough’s test to an observed data set is as follows: The data set is divided into contiguous segments of $N$ cycles each. Then the ratio of the averages of the empirical variances is calculated and compared with the ratio of the expectation values, plotted as functions of $N$ in Figure 6.

Weisshaar et al (2023) augmented the method by determining suitable confidence intervals through Monte Carlo simulations for the clock and the random phase cases, assuming normally distributed variations in cycle length. They applied the test to the extended sunspot record of now 28 cycles, four more than available to Gough. The main improvement is narrower confidence intervals, rejecting the synchronization hypothesis on a 2$\sigma$ level (Figure 5).

Recently, a reconstruction of yearly sunspot numbers from the record of cosmogenic 14C in tree rings for the years 976 until 1888 (Brehm et al, 2021; Usoskin et al, 2021) has extended the number of contiguous cycles available for the analysis to 84. The Gough test confirms the previous result based on the direct sunspot record, in fact strengthening it significantly, since now the synchronization hypothesis can be rejected even on a $>3\sigma$ level (Figure 6).

Weisshaar et al (2023) also applied the method of Eddington and Plakidis (1929) mentioned above to these new data and found, consistent with the analysis discussed here, that the fraction of clock-synchronised perturbations is negligible.

The question may arise how misidentifications of the observed solar cycles can affect the results. If this happens not too common, the nature of the fluctuations (phase stability or migration) is not expected to be changed by this bias. As a test, a lost cycle between more distant minima was “restored” by placing a minimum in between. This did not cause the $S$-values to leave the phase migration confidence interval.

Furthermore, the above-mentioned method of the phase evolution of empirical cycle data is therefore consistent with a random walk (such as provided by a memory-less dynamo process). External synchronization by a ‘clock’ is clearly excluded at a high significance.

6 Implications for the dynamo theory

The solar magnetic cycle is maintained by a dynamo process, operating in the solar convection zone (SCZ). Thus, it is natural to expect that the variations in the solar cycle are caused by some mechanisms in the solar dynamo. Here we identify the causes of the variations in the solar cycle and demonstrate them by presenting results from some illustrative models. Let us first summarise the mechanism of the solar dynamo.

6.1 Introduction to the solar dynamo

There is enough evidence that the solar dynamo is a mechanism in which toroidal and poloidal fields sustain each other through a cyclic loop (e.g., Parker, 1955; Cameron and Schüssler, 2015). In this loop, the toroidal field is generated due to the shearing of the poloidal field by the differential rotation in the deeper CZ. The toroidal field rises to the surface due to magnetic buoyancy to give rise to sunspots or more generally bipolar magnetic regions (BMRs). These BMRs are systematically tilted with respect to their East-West orientations. Due to these tilts, after their decay, BMRs produce a poloidal field. This, the so-called Babcock–Leighton process is clearly identified in the observed magnetic field data on the solar surface (e.g., Mordvinov et al, 2022). The observed correlation between the polar field (or its proxy) at the solar minima and the amplitude of the next cycle (Wang and Sheeley, 2009; Kitchatinov and Olemskoy, 2011; Muñoz-Jaramillo et al, 2013; Priyal et al, 2014) and the flux budgets of the observed and the generated poloidal and toroidal fields (Cameron and Schüssler, 2015) suggest that the Babcock–Leighton process is possibly the main source of the poloidal field in the Sun.

There is however another mechanism through which the poloidal field in the sun can be produced and that is the classical $\alpha$ effect as originally proposed by Parker (1955) and mathematically formulated by Steenbeck et al (1966). In this mechanism, the toroidal field is twisted by the helically rising blobs of plasma in the SCZ. However, this process of lifting and twisting of the field by the convective flow experiences catastrophic quenching due to helicity conservation and thus this process operates when the energy density of the toroidal field is less than the energy density of the convective motion (Sec. 8.7 of Brandenburg and Subramanian, 2005). Therefore, this $\alpha$ effect is unfavourable in the solar convection zone and the obvious option is to consider the observationally supported Babcock–Leighton process for the generation of the poloidal field in the sun.

To study the dynamo action, we need to begin with at least following two fundamental equations of magnetohydrodynamics (MHD).

[TABLE]

[TABLE]

where $\boldsymbol{B}$ and $\boldsymbol{v}$ are the magnetic and velocity fields, respectively, $\eta$ is the magnetic diffusivity, $\rho$ is the density, $P$ is the pressure, \boldsymbol{J}=\mbox{\boldmath\nabla}\times\boldsymbol{B}/\mu_{0}, the current density, $\nu$ is the kinetic viscosity, S_{ij}=\frac{1}{2}(\nabla_{i}v_{j}+\nabla_{j}v_{i})-\frac{1}{3}\delta_{ij}\mbox{\boldmath\nabla}\cdot\boldsymbol{v} is the rate-of-strain tensor, and the term $\boldsymbol{F}$ includes gravitational, Coriolis and any other body forces acting on the fluid. These equations along with the mass continuity and energy equations and equation of state are numerically solved with appropriate boundary conditions in the solar CZ to study the dynamo problem. Broadly there are two approaches for doing this, namely, the global MHD simulations and mean-field modellings. In global MHD simulations, we solve the above MHD equations numerically to resolve the full spectrum of turbulent convection. In mean-field models, we study the evolution of the mean/large-scale quantities by parameterizing the small-scale/fluctuating quantities using suitable approximations.

Global MHD simulations for the Sun are challenging due to extreme parameter regimes, such as high fluid and magnetic Reynolds numbers and large stratification. Despite these, simulations have begun to produce some solar-like features; see Section 6 of Charbonneau (2020). However, due to their computationally expensive nature, these simulations were rarely run for many cycles so that the cycle variabilities can be studied. Passos and Charbonneau (2014) have produced simulations for several cycles and shown long-term modulations (also see Karak et al, 2015, for a simulation at solar rotation rate although ran for not many cycles). Augustson et al (2015) and Käpylä et al (2016) performed MHD convection simulations for the cases of three and five times the solar rotation rate, respectively. They both found an episode of suppressed surface activity, somewhat resembling the solar grand minimum. Although these results of cycle modulations are encouraging, simulations face serious issues when matching with observations, for example, concerning solar observations, simulations (i) produce higher power at the largest length scale, (ii) do not produce BMRs, and (iii) do not produce correct large-scale flows, particularly, they produce a large variation in the differential rotation.

On the other hand, mean-field models are computationally less expensive and easy to analyse their results. Probably due to these reasons, long-term modulations are studied using mean-field dynamo models. Due to the observational facts that the magnetic field at the solar minima and the large-scale velocity field are largely axisymmetric, historically the mean-field models are constructed under axisymmetric approximation. With this approximation, the equations for the poloidal and toroidal fields are written as

[TABLE]

[TABLE]

where $A$ is the potential for the poloidal field (\boldsymbol{B_{p}}=\mbox{\boldmath\nabla}\times(A\boldsymbol{\hat{\phi}}), $B$ is the toroidal field, $s=r\sin{\theta}$, $\boldsymbol{v_{m}}(=v_{r}\boldsymbol{\hat{r}}+v_{\theta}\boldsymbol{\hat{\theta}}$) represents the meridional circulation, $\eta_{t}$ is the turbulent diffusivity which is assumed to depend only on $r$, $\alpha$ is the $\alpha$ effect, and ${\rm\Omega}$ is the angular frequency.

The term $\alpha B$ in Equation (7) is the source for the poloidal field through the $\alpha$ effect. The generation of the poloidal field through the Babcock–Leighton process is also parameterised in the 2D (axisymmetric models) through the same $\alpha B$ term. However, this $\alpha$ operates near the surface of the sun and it has a completely different origin than the $\alpha$ effect which operates in the whole convection zone due to helical convection. In comprehensive 3D dynamo models (Yeates and Muñoz-Jaramillo, 2013; Miesch and Dikpati, 2014; Miesch and Teweldebirhan, 2016; Kumar et al, 2019; Bekki and Cameron, 2022), this $\alpha B$ term is not added in Equation (7), instead, explicit BMRs are deposited whose decay produces a poloidal field. The source for the toroidal field in Equation (8) is due to the ${\rm\Omega}$-effect which is represented by the term: s(\boldsymbol{B_{p}}\cdot\mbox{\boldmath\nabla}){\rm\Omega}. The above equations technically represent the equations for the $\alpha{\rm\Omega}$ dynamo model, in which the generation of the toroidal field through the $\alpha$ effect is assumed to be much less than the generation due to $\Omega$ effect, which is true in the sun; see e.g., Cameron and Schüssler (2015).

6.2 Causes for long-term variations in the solar activity

With the above discussion of the solar dynamo, we now identify the causes of the cycle modulation. As the solar dynamo is nonlinear, it is natural to expect that the modulation in the solar cycle is caused by the back reaction of the flow on the magnetic field. Therefore, we first identify the nonlinearities in the dynamo models and check if they can lead to cycle modulations.

6.2.1 Nonlinearities in the dynamo

As we can see from Equation (6), the magnetic field can alter the flow directly through the Lorentz force. The Lorentz force can come from the mean magnetic field and the mean current (which is popularly known as the Malkus-Proctor effect (Malkus and Proctor, 1975) in the mean-field context) and from the fluctuating magnetic field and the current. The mean magnetic field can also alter the anisotropic convection which is responsible for transporting angular momentum and maintaining differential rotation and meridional flow in the Sun (Kitchatinov et al, 1994b). This effect is also called micro-feedback. When these Lorentz feedbacks of the magnetic fields are included in the flow, we expect a long-term modulation in the flow and the magnetic cycle.

In mean-field models, the magnetic feedback is captured by considering a direct Lorentz force of the mean magnetic field in the zonal flow (e.g., Bushby, 2006) and/or by a quenching term in the $\Lambda$ effect (e.g., Küker et al, 1999). Cycle modulations in these systems can generally happen in two ways. In the first one, the magnetic energy of the primary mode (the equatorial symmetry or antisymmetric) can oscillate due to the energy exchange between the flow and the magnetic field via the nonlinear Lorentz feedback. In this case, a considerable amount of modulation in the differential rotation is observed. In the second case, a small magnetic perturbation on the differential rotation can slowly change one dominant dynamo mode into another. In this case, the magnetic field parity can change (between equatorially symmetric (quadrupole) and antisymmetric (dipole)) without producing a large change in the differential rotation. These two mechanisms are respectively coined as Type II and I modulations. Mean-field models have demonstrated that nonlinear back reaction of magnetic field on large-scale flow through these types of modulations can induce a variety of modulation patterns in the cycle amplitude, including grand minima and parity modulations which do not leave a strong imprint in differential rotation (e.g., Beer et al, 1998; Knobloch et al, 1998; Bushby, 2006; Weiss and Tobias, 2016). Both types of modulation can arise in a model, however, as the observed differential rotation shows a tiny variation over the solar cycle, we expect the Type II modulation is less likely to occur in the Sun. Even for Type I modulation, a detailed comparison of the magnetic field and the flows in these models with the observations is missing (also see Section 7 of Charbonneau, 2020, for a discussion on this topic).

Next is the meridional flow, which is the second important large-scale flow in the Sun. As it arises due to a slight imbalance between the non-conservative centrifugal and buoyancy forces, we expect its large variation. In fact, the global simulations find a large variation in the meridional flow despite a small variation in the differential rotation (Karak et al, 2015). In Babcock–Leighton type dynamo models, meridional circulation plays a crucial role in transporting the field on the surface from low to high latitudes and down to the deeper CZ where the shear produces a toroidal field. The toroidal field is transported to the low latitudes via the equatorward return flow and possibly causes the equatorward migration of the sunspot belt. Thus, in these models, meridional circulation largely regulates the cycle period (Dikpati and Charbonneau, 1999; Karak and Choudhuri, 2011). It also affects the strength of the field as a weak meridional circulation allows the field to advect slowly and gives more time for diffusion (Yeates et al, 2008). Karak (2010) showed that when a variable meridional flow is used in a high diffusivity dynamo model to match the observed solar cycle periods, the amplitudes of the cycles are also modelled up to some extent (also see Karak and Choudhuri, 2011; Hazra et al, 2015, for modelling various aspects of solar cycle using variable meridional flow). In an extreme case, a largely reduced meridional circulation can trigger a Maunder-like grand minimum. In reality, how large the variation in the meridional flow occurred in the past remains uncertain. However, it is obvious that any changes in the flow can lead to modulation in the solar cycle.

Turbulent transport as parameterized by, for example, the turbulent diffusivity, $\Lambda$ effect, and heat diffusion are also nonlinear because the Lorentz force of the small-scale as well as the large-scale dynamo-generated fields act on the small-scale turbulent flows. However, due to limited knowledge in the turbulence theory for solar parameter regions, we do not have a satisfactory model for the magnetic field-dependent form of the turbulent transport parameters; however, see Ruediger and Kichatinov (1993) and Kitchatinov et al (1994a) respectively, for the magnetic field-dependent forms of $\alpha$ and $\eta$ based on the quasi-linear approximation.

Finally, the toroidal to poloidal part of the dynamo loop involves some nonlinearities. When the generation of poloidal field is due to the classical $\alpha$ effect, there is a well-known $\alpha$ quenching of the form $1/\left(1+(B/B_{\rm eq})^{2}\right)$ with $B_{\rm eq}$ being the equipartition field strength. However, this type of $\alpha$ quenching tries to make a stable cycle rather than producing irregularity in the cycle. In the Babcock–Leighton dynamo, the generation of the poloidal field from the toroidal one also involves several nonlinearities. Here we discuss the following three potential candidates for these.

• •

Flux loss due to magnetic buoyancy

The magnetic buoyancy as proposed by Parker (1955) plays a critical role in the emergence of BMRs on the solar surface. As the shearing of the poloidal field due to differential rotation intensifies the strength of the toroidal field, there comes a point where the magnetic energy density of the toroidal flux tubes becomes greater than the kinetic energy of the local convective plasma inside the CZ, as a result, the flux tubes become buoyant and start rising through the CZ, eventually giving birth to the sunspots. Following this process, the strength of the magnetic field gets locally reduced as a part of it rises due to buoyancy and the flux tube becomes inefficient to produce further sunspots for some time (however see a counter-argument by Rempel and Schüssler, 2001). The sharp rise in the flux loss once the toroidal field strength exceeds a certain value clearly indicates a nonlinear mechanism in the solar dynamo. Incorporating this mechanism of toroidal flux loss due to buoyancy in a simple manner, Biswas et al (2022) showed that this nonlinear process plays a critical role in limiting the growth of the solar dynamo which is a potential mechanism to explain why different solar cycles rise differently depending on their strength but all the solar cycles decay with similar statistical properties (see Figure 7). They found that introducing the flux loss in the dynamo simulations was critical to reproduce the long-term features of the latitudinal distribution of the sunspots (Waldmeier, 1955; Cameron and Schüssler, 2016); also see Cameron and Schüssler (2016) and Talafha et al (2022) for an alternative explanation of the universal decay of the solar cycle using cross-equatorial diffusion.

• •

Latitude quenching

It has been found that when BMRs appear in low latitudes, the leading polarities from both hemispheres get efficiently cancelled at the equator. This leads to the following polarities of the BMRs efficiently getting carried to the poles and contributing to the polar field, see Figure 8. On the other hand, BMRs appearing in the high latitudes do not exhibit efficient cross-hemisphere cancellation and thus do not contribute significantly to the polar field (Jiang et al, 2014; Karak and Miesch, 2018). It is seen that strong cycles produce more BMRs at high latitudes. In other words, the average latitude of the BMRs is high for the strong cycles (Solanki et al, 2008; Mandal et al, 2017). Hence for a strong cycle, most of its BMRs emerging at high latitudes would be less efficient in polar field production and vice versa for the weak cycles. This mechanism, so-called the latitude quenching (Petrovay, 2020) may help to stabilize the growth of the magnetic field in the Sun (Jiang, 2020).

Introducing a latitude-dependent threshold on the BMR emergence condition into a 3D Babcock–Leighton dynamo simulation, Karak (2020) showed that latitude quenching can regulate the growth of a magnetic field when the dynamo is not too supercritical.

• •

Tilt quenching

The tilt angle of BMR plays a crucial role in generating poloidal field in the Sun. For a given latitude, the amount of generated poloidal field increases with the increase of tilt. The thin flux tube model for the sunspot formation suggests that the tilt of the BMR is produced due to a torque acting on the diverging flows produced from the apex of the rising flux tube which forms the BMR (D’Silva and Choudhuri, 1993; Fan et al, 1994). Thus, if the magnetic field of the sunspot-forming flux tube is strong, then it will rise quickly and the Coriolis force will get less time to induce tilt. In a strong cycle, the toroidal magnetic field is strong and the number of BMRs with strong magnetic field tends to be high (Jha et al, 2020). Thus, we expect the mean tilt in that cycle to be smaller. A lesser tilt will produce less poloidal field and the next cycle will be weak. Hence, this may be a potential mechanism for stabilizing the growth of the magnetic cycle through the reduction of tilt which is known as the tilt quenching.

The observational evidence of tilt quenching is limited. Dasi-Espuig et al (2010); Jiao et al (2021) showed that there is a statistical anti-correlation between the cycle-average tilt of the sunspots with the cycle strength (Figure 9a). On the other hand, Jha et al (2020) examined the variation of BMR tilt with the strength of its magnetic field within a cycle. They found a non-monotonous dependence of the tilt with the BMR field strength as seen in Figure 9(b). For weak field strengths, the tilt first increases, however at sufficiently strong field strengths, the BMR tilt starts to decrease.

6.2.2 Stochastic effects in the dynamo

The solar convection zone is turbulent and thus the turbulent quantities (such as $\alpha$ effect) are subject to fluctuate around their means. Hoyng (1993) showed that as there are finite numbers of convection eddies along the longitudes in the sun, the fluctuations of the turbulent transport coefficients can be larger than their means. There is a long history including the stochastic noise in the $\alpha$ effect in the mean-field dynamo models. Most of these studies find long-term modulations in the cycle and grand minima in a certain parameter range of the dynamo number (Choudhuri, 1992; Ossendrijver and Hoyng, 1996; Ossendrijver et al, 1996; Gómez and Mininni, 2006; Brandenburg and Spiegel, 2008; Moss et al, 2008).

In Babcock–Leighton dynamo also stochastic fluctuations are unavoidable. The toroidal to poloidal part of this model primarily involves stochastic fluctuations due to the following effects.

• •

Scatter around Joy’s law

Observations find that the tilt “statistically” increases with the increase of latitude, which is known as Joy’s law. However, a large number of BMRs do not follow this relation (so-called non-Joy), as seen by a huge scatter around the mean trend in Figure 10. In fact, there are many BMRs which are of anti-Hale type. These anti-Hale and non-Joy BMRs, having opposite tilts (negative in the northern hemisphere) are responsible for generating opposite polarity field (with respect to the expected polarity) and lead to large fluctuations in the polar field (Jiang et al, 2014; Hazra et al, 2017; Nagy et al, 2017; Mordvinov et al, 2022).

• •

Variations in the BMR eruption rates

There are spatial and temporal variations in the BMR eruptions. BMRs near the equator are much more efficient in generating poloidal field in the Sun because for them the leading polarity can easily connect with the opposite polarity flux from the opposite hemisphere (Cameron et al, 2013; Jiang et al, 2014; Karak and Miesch, 2018; Karak, 2020; Mordvinov et al, 2022). Thus variation in the latitudinal position can produce variation in the generated poloidal field. Next, the rate of BMR eruption is not the same—there is a distribution. Thus, the rate of generation of the poloidal field is not the same (Karak and Miesch, 2017). Furthermore, the flux contents of the BMR has also a distribution and thus a wrongly tilted BMR with high flux can disturb the polar field in the sun considerably (Nagy et al, 2017).

In summary, the randomness involved in the BMR properties (originated due to the turbulent nature of the convection) produces variation in the poloidal field. Although the sun produces thousands of spots in a cycle, only a few spots are produced (on average) per day. This leads to variations in the polar field comparable to its mean value. In the next section, we shall demonstrate some illustrative results from stochastically driven Babcock–Leighton dynamo models.

6.3 Babcock–Leighton dynamo models for the long-term variation

As discussed above, the generation of the poloidal field in the Babcock–Leighton dynamo models involves some randomness. Thus, in axisymmetric dynamo models, these randomnesses were captured by adding a noise term in the poloidal source (e.g., Charbonneau and Dikpati, 2000). Long-term modulations, including Gnevyshev-Ohl/Odd-Even rule (Charbonneau, 2001; Charbonneau et al, 2007) and grand minima (Charbonneau et al, 2004; Choudhuri and Karak, 2009; Passos et al, 2012, 2014) are naturally produced in these models. Variations within the cycle, like the amplitude-period anti-correlation (Charbonneau and Dikpati, 2000; Karak, 2010) and Waldmeier effect (Karak and Choudhuri, 2011; Biswas et al, 2022) are also reproduced. Karak et al (2018) showed that a large variation in the Babcock–Leighton process can change the polar field abruptly and this can lead to double peaks in the following cycle. While in most of the studies, the level of fluctuations was tuned to produce the observed variation of the solar cycle including a reasonable number of grand minima, Choudhuri and Karak (2012) and Olemskoy and Kitchatinov (2013) made some estimate of the fluctuations in the Babcock–Leighton process from observations. Choudhuri and Karak (2012) found the correct frequency of grand minima as observed in the cosmogenic data for the last 11,000 years. Olemskoy and Kitchatinov (2013) showed that the statistics of grand minima are consistent with the Poisson random process, indicating the initiation of grand minima to be independent of the history of the past minima.

In recent years, cycle modulations were, in particular, produced by including the variations in the BMR properties in two comprehensive models, namely, 2$\times$2D (Lemerle and Charbonneau, 2017) and 3D dynamo models (Karak and Miesch, 2017). In Figure 11, we show cycles from the 3D dynamo model presented by Karak and Miesch (2017). As seen in Figure 11(a), the variation in the BMR emergence rate and the flux distribution produce little variation in the solar cycle. When the variation around Joy’s law tilt is included, it produces a large variation, including suppressed magnetic activity like the one seen during Dalton minimum and Maunder minimum as shown in Figure 11b (the regions shaded in green). Here, the grand minima are identified in the same manner as done in the observed data (Usoskin et al, 2007), i.e., the modelled-sunspot data are first binned in 10 years window and smoothed and then a grand minimum is considered when the smoothed data fall below $50\%$ of the average at least for two cycles.

In Figure 12, we present a detailed view of a grand minimum. We find that some of the observed features of the Maunder minimum (hemispheric asymmetry, gradual recovery, slightly longer cycle) are reproduced in this figure. We note that during this grand minimum, some BMRs are still produced, the number of which is a bit larger than that was observed during Maunder minimum (Usoskin et al, 2015; Vaquero et al, 2015; Zolotova and Ponyavin, 2016; Carrasco et al, 2021). However, we should keep in mind that the observations during Maunder minimum were limited (due to the poor resolving power of the 17th-century telescopes) to detect the small BMRs (e.g., Vaquero and Vázquez, 2009); only big sunspots were detected. In our Babcock–Leighton dynamo model, few BMRs erupt which produces a poloidal field at a slow rate through the Babcock–Leighton process and the model emerges from the grand minimum episode. It is the downward magnetic pumping included in our model which helps to reduce the magnetic flux loss through the surface and recovers the model from grand minima (Cameron et al, 2013; Karak and Cameron, 2016).

There have been suggestions that during Maunder-like extended grand minima, the Babcock–Leighton process may not operate due to few observed sunspots, and $\alpha$ effect (Parker, 1955) is the best candidate for this as it efficiently operates in sub-equipartition field strength (Karak and Choudhuri, 2013; Passos et al, 2014; Ölçek et al, 2019). We observe that our model also fails to recover when it enters a deep grand minimum and stops producing BMRs due to the fall of the toroidal field below the threshold for BMR formation. However, this happens very rarely. While it is a critical question to answer what mechanism dominates in recovering the Sun from an extended grand minimum, it is expected that Babcock–Leighton process becomes less efficient during this phase and the $\alpha$ effect certainly helps in recovering the Sun from grand minima.

Dynamo models with stochastic fluctuations also produce grand maxima. Our model presented in Figure 11b also produces a few grand maxima shown by the regions shaded in red. Similar to the grand minima, grand maxima are also computed based on the smoothed sunspot number, but here the threshold is taken as $150\%$ of the long-term mean. Systematic studies of grand maxima using dynamo models are limited (however, see Karak and Choudhuri, 2013; Olemskoy and Kitchatinov, 2013; Inceoglu et al, 2017). Kitchatinov and Olemskoy (2016) showed that at the beginning of the cycle, if the generation of the poloidal field is reversed (say due to the emergence of some wrongly tilted BMRs), then it will amplify the existing polar field, instead of reversing it. This increase in the magnetic field can lead to a grand maximum. Another mechanism of grand maxima was given by Ölçek et al (2019), who showed that when the deep-seated $\alpha$ effect is coupled with the surface Babcock–Leighton source, then these two sources more or less contribute equally to generate a strong poloidal field through a sort of constructive interference.

Finally, for the secular and supersecular modulations (modulations beyond 11-year periodicity, e.g., Gleissberg cycle, Suess/de Vries cycle, Eddy cycle, and 2400-year Hallstatt cycle; Beer et al, 2018), there are limited studies available in the literature. In a simplified \alpha$$\Omega dynamo model coupled with the angular momentum equation, Pipin (1999) found the Gleissberg cycle as a result of the re-establishment of differential rotation after the magnetic feedback on the angular momentum transport. Cameron and Schüssler (2017) modelled the overall power spectrum of solar activity using a generic normal form model for a noisy and weakly nonlinear limit cycle, and Cameron and Schüssler (2019) showed that the long-term modulations beyond the 11-year cycle are consistent with the realization noise, thus casting doubt whether secular and supersecular modulations are connected to the intrinsic periodicities of the solar dynamo.

7 Summary

Herewith, a brief overview is presented of the long-term variability of solar activity at centennial – millennial timescales. The main feature of solar variability is the 11-year quasi-periodic Schwabe cycles, which is however variable per se in both magnitudes, duration and phase. While the direct telescopic observations of the Sun cover roughly four centuries since 1610 and cover a full range of solar-activity levels from the Maunder minimum in the 17th century to the Modern grand maximum in the late 20th century, the quality of the sunspot-number dataset is inhomogeneous and greatly degrades back in time, being quite imprecise before $\approx$1820s. Moreover, it is too short to study the statistical properties of the solar-cycle modulation on a long timescale.

The cosmogenic-isotope method provides quantitative reconstructions of solar activity on the multi-millennial timescale with stable quality throughout ages making it possible to study long-term solar-cycle modulation. Using the decadal data for the Holocene (the last twelve millennia), it is possible to identify specific observed properties of solar variability beyond the Schwabe cycle:

• •

The Sun spends about ${}^{1}/_{6}$ of its time in the grand minimum state, grand minima tend to cluster with a $\approx$210-year recurrence time;

• •

The Sun spends about ${}^{1}/_{10}$ of its time in the grand maximum state, grand maxima appear without any regular pattern;

• •

During the major fraction of time, the Sun is in the cyclic moderate activity state;

• •

Several quasi-periodicities can be found in long-term solar variability, but they are intermittent and barely significant: Centennial Gleissberg cycle which is an oscillation with the characteristic time of 60 – 140 years; 210-year Suess/de Vries cycles manifesting itself through intermittent recurrence of grand minima; About 2400-year Hallstatt cycle whose nature is still unclear; Other long-term cycles, including the millennial Eddy cycle, are insignificant.

A recent reconstruction of the annual sunspot numbers from high-precision radiocarbon data for the last millennium makes greatly extended, nearly tripling, the statistic of solar cycles to 96 individually resolved cycles. In particular, the Waldmeier rule (high cycles rise faster) is statistically confirmed on a larger statistical basis, while the Genvyshev-Ohl rule of the even-odd cycle pairing is not confirmed. The extended statistic of solar cycles has made it possible, for the first time, to answer the question principle to the solar dynamo theory: is the solar cycle phase-locked, implying an intrinsic synchronisation process as proposed by some external clocking mechanisms, or is random and incoherent. The new analysis excludes the phase-locking hypothesis at a high significance level, implying that solar cycles vary randomly.

A brief review of the theoretical perspectives to explain the observed features in the framework of the dynamo models is presented. It is discussed that the nonlinearities in the dynamo, including the effects of the flux loss due to magnetic buoyancy as well as latitude and tilt quenching, help to stabilize the solar dynamo, rather than producing variability in the solar cycle. Primary causes of the solar cycle variability are the stochastic fluctuations in the dynamo which are inherent in different processes such as a large scatter of the BMR’s tilts around Joy’s law, and variability in the BMR eruption rates and locations. On one hand, while modern dynamo models are able to reproduce, with a reasonable ad-hoc tuning of the parameters, the observed features of solar variability, the exact role of those factors is not clear, and some discrepancies between the model results and the data still remain. On the other hand, the progress in the accuracy of models is significant, and we keep gaining knowledge of the processes driving solar variability with the new data acquainted and new models developed.

\bmhead

Acknowledgments IU acknowledges the Academy of Finland (project ESPERA No. 321882). AB and BBK gratefully acknowledge the financial support provided by ISRO/RESPOND (project No. ISRO/RES/2/430/19-20), the Department of Science and Technology (SERB/DST), India through the Ramanujan Fellowship (project No. SB/S2/RJN-017/2018), the International Space Science Institute (ISSI, Team 474), and the computational resources of the PARAM Shivay Facility under the National Supercomputing Mission, the Government of India, at the Indian Institute of Technology Varanasi. This work was performed in the framework of the ISSI workshop “Solar and Stellar Dynamos: A New Era”.

\bmhead

Competing Interests The authors declare they have no conflicts of interest.