Understanding solar cycle variability

TL;DR

This paper presents a normal-form dynamo model that explains solar cycle variability, including grand minima, by combining dynamo theory with empirical data, providing a unified understanding of solar magnetic activity fluctuations.

Contribution

It introduces a normal-form model constrained by observations that captures solar cycle variability and grand minima, validated against Babcock-Leighton dynamo models.

Findings

Model reproduces solar activity variations from decades to millennia.

Successfully explains occurrence of grand minima.

Aligns with results from Babcock-Leighton dynamo models.

Abstract

The level of solar magnetic activity, as exemplified by the number of sunspots and by energetic events in the corona, varies on a wide range of time scales. Most prominent is the 11-year solar cycle, which is significantly modulated on longer time scales. Drawing from dynamo theory together with empirical results of past solar activity and of similar phenomena on solar-like stars, we show that the variability of the solar cycle can be essentially understood in terms of a weakly nonlinear limit cycle affected by random noise. In contrast to ad-hoc `toy models' for the solar cycle, this leads to a generic normal-form model, whose parameters are all constrained by observations. The model reproduces the characteristics of the variable solar activity on time scales between decades and millennia, including the occurrence and statistics of extended periods of very low activity (grand minima). Comparison with results obtained with a Babcock-Leighton-type dynamo model confirms the validity of the normal-mode approach.

Equations

$$\frac{\mathrm{d}X}{\mathrm{d}t}-(\beta +i{\omega }_0)X+({\gamma }_r+i{\gamma }_i)|X{|}^{2}X=0,$$

$$\mathrm{d}X=\left(\beta +i{\omega }_0-({\gamma }_r+i{\gamma }_i)|X{|}^2\right)X\mathrm{d}t+\sigma X\mathrm{d}{W}_c=0,$$

$$a(\theta ,t)=\frac{1}{\sin \theta }{\int }_0^{\theta }\sin \theta R_{\odot }^2{B}_r{|}_{R=R_{\odot }}\mathrm{d}\theta ,$$

$$b(\theta ,t)={\int }_{R_b}^{R_{\odot }}{B}_{\varphi }r\mathrm{d}r,$$

$$\mathrm{d}a$$

$$\mathrm{d}b$$

$$a_{\mathrm{S}}(\theta ,t)=\frac{{\alpha }_0}{1+b^2/b_c^2}\cos \theta \sin \theta b(\theta ,t),$$

Full text

**

The level of solar magnetic activity, as exemplified by the number of sunspots and by energetic events in the corona, varies on a wide range of time scales. Most prominent is the 11-year solar cycle, which is significantly modulated on longer time scales. Drawing from dynamo theory together with empirical results of past solar activity and of similar phenomena on solar-like stars, we show that the variability of the solar cycle can be essentially understood in terms of a weakly nonlinear limit cycle affected by random noise. In contrast to ad-hoc ‘toy models’ for the solar cycle, this leads to a generic normal-form model, whose parameters are all constrained by observations. The model reproduces the characteristics of the variable solar activity on time scales between decades and millennia, including the occurrence and statistics of extended periods of very low activity (grand minima). Comparison with results obtained with a Babcock-Leighton-type dynamo model confirms the validity of the normal-mode approach.

Sun: magnetic fields, Sun: activity

1 Introduction

Apart from its 11-year (quasi)periodicity, the most striking property of the solar activity record is the marked variability of the cycle amplitudes (cf. top panels of Fig. 1), including extended intervals of very low or particularly high activity (grand minima and maxima, see Usoskin, 2017). Understanding the nature of the variability is a prerequisite for sensible attempts to predict future activity levels. Therefore, we need to clarify to what extent randomness, intrinsic periodicities apart from the 11-year cycle, and nonlinearities of the underlying dynamo process generating the solar magnetic field contribute to the observed long-term variability of solar activity.

There exists a rich literature describing attempts to understand the variability in the framework of hydromagnetic dynamo theory, a full review of which is beyond the scope of this paper (see, e.g., Tobias, 2002; Charbonneau, 2010, 2014). Such studies can be roughly divided into two approaches: (1) nonlinear dynamics and deterministic chaos, and (2) random fluctuations of dynamo excitation. The nonlinear dynamics approach typically considers the bifurcation structure of low-order dynamical systems (see reviews by Weiss, 1990; Tobias et al., 1995; Lopes et al., 2014), but models based on nonlinear PDEs have also been investigated frequently (e.g., Schmitt & Schuessler, 1989; Tobias, 1997; Bushby, 2006). Studies assuming random fluctuations reach from the minimalistic ‘model’ of Barnes et al. (1980) to detailed considerations of the mode structure of stochastically excited dynamos (Hoyng, 1993; Hoyng & van Geffen, 1993; Ossendrijver & Hoyng, 1996). Most of these studies adopt an $\alpha\Omega$-type dynamo approach with random fluctuations of the $\alpha$-effect thought to result from the non-stationary nature of solar convection (e.g. Choudhuri, 1992; Moss et al., 1992; Ossendrijver et al., 1996). There are also studies combining nonlinear dynamo models with random fluctuations, either in the framework of low-order dynamical systems (e.g. Mininni et al., 2001; Passos & Lopes, 2011) or assuming detailed dynamo models (e.g. Charbonneau & Dikpati, 2000; Mininni & Gómez, 2002; Moss et al., 2008; Lemerle & Charbonneau, 2017). Grand minima can also result from of ‘on-off intermittency’ due to the interaction of two spatially separated dynamos (Platt et al., 1993; Schmitt et al., 1996; Passos et al., 2014). Recently, Olemskoy et al. (2013) have estimated the fluctuating source term of a Babcock-Leighton-type dynamo considering the observed scatter of the tilt angles of sunspot groups. Assuming weakly supercritical dynamo excitation, the simulated long-term evolution of solar activity exhibits grand minima whose statistics are consistent with the actual solar record (see also Kitchatinov & Olemskoy, 2016). Global 3D-MHD simulations of convection and magnetic field in a rotating spherical shell typically show strong variability of dynamo-generated magnetic field and may also provide cyclic fields (e.g. Augustson et al., 2015; Passos & Charbonneau, 2014; Käpylä et al., 2016; Hotta et al., 2016; Fan & Fang, 2016). However such simulations remain far from the Sun in terms of various parameters, and so the variability seen in the simulations is not yet directly relatable to that of the Sun.

Although no definite quantitative model of the global dynamo exists so far, the basic ingredients are uncontroversial. The systematic properties of sunspot groups (Hale et al., 1919) indicate that they originate from a reservoir of organized East-West orientated (toroidal) field in the solar convection zone. This field is generated by winding up a poloidal field (such as a dipole field aligned with the rotation axis) by the differential rotation of the Sun, so that its axisymmetric component dominates. The poloidal field is (re)generated against the effect of Ohmic decay by the collective effect of loops formed from the toroidal field by convective flows and/or magnetic buoyancy. The loops become twisted owing to the Coriolis force and thus acquire a systematic meridional component (Parker, 1955; Babcock, 1961; Steenbeck & Krause, 1966). This interplay of toroidal and poloidal magnetic field leads to a 22-year magnetic cycle and an 11-year cycle of sunspot activity.

Here we consider these basic ingredients of solar dynamo theory together with solar and stellar observations to elucidate the nature of the variability of the solar cycle. We argue that the dynamo works near its excitation threshold (marginal state), so that a normal-form model representing the essence of the underlying dynamo process can be used. This model is generic in the sense that it is valid for any weakly nonlinear system in the vicinity of a supercritical Hopf bifurcation, independent of the nature of the nonlinearity. The very setup of the model and its parameters are all constrained by observations.

2 Generic normal-form model

Irrespective of the details of most models for the solar dynamo proposed so far, the corresponding systems of equations for the magnetic field components show qualitatively similar behaviour near the onset of dynamo action, which is governed by a control parameter (often called ‘dynamo number’) involving differential rotation, Coriolis effect, and magnetic diffusion. As long as the control parameter remains below a critical value, the stationary solution has zero magnetic field and the dynamo is not excited. When the critical value of the control parameter is exceeded, the system exhibits a periodic solution. This behaviour is generic for almost all models of the global solar dynamo (technically named $\alpha\Omega$ dynamos). In the language of dynamical systems, the periodic solution emerges from a (supercritical) Hopf bifurcation: a fixed point becomes unstable and spawns a limit cycle (Guckenheimer & Holmes, 1983; Tobias et al., 1995).

The study of solar-like stars has revealed that the level of magnetic activity systematically declines with decreasing stellar rotation rate (e.g., Reiners, 2012; Reiners et al., 2014). In particular, the relatively slow rotation rate of the Sun appears to put it near to the threshold for which global dynamo action ceases (van Saders et al., 2016; Metcalfe et al., 2016). This allows us, irrespective of the nature of the nonlinearity that limits the amplitude of the cycle, to describe the solar dynamo generically by the normal form of a weakly nonlinear system near a Hopf bifurcation, which is independent of the nature of the nonlinearity (Arnol’d, 1972; Guckenheimer & Holmes, 1983), viz.

[TABLE]

where $X$ is a complex variable. Its real and imaginary components are related to the toroidal and poloidal components of the magnetic field in a manner that depends on the kind of nonlinearity considered (e.g., back reaction on the differential rotation or quenching of the regeneration term for the poloidal field). We may therefore consider the real or the imaginary part of $X$ to represent the cyclically varying field magnitude. For simplicity, we take the activity level as quantified by the sunspot number, SSN, to be proportional to the absolute value of the field and write ${\rm SSN}=|{\rm Re}(X)|$, thus scaling $|{\rm Re}(X)|$ in units of the sunspot number. This quantity represents the 11-year sunspot cycle, in contrast to the 22-year magnetic cycle reflected in ${\rm Re}(X)$ and ${\rm Im}(X)$.

The quantity $\beta$ in Eq. (1) is the linear growth rate of the cycle amplitude (related to the dynamo excitation, i.e., supercriticality of the dynamo number) and $\omega_{0}$ is the (magnetic) cycle frequency in the limit of zero amplitude. The parameters $\gamma_{r}$ and $\gamma_{i}$ generically represent the nonlinearity of the system. They determine the cycle amplitude, $|X|=\sqrt{\beta/\gamma_{r}}$, and the nonlinear cycle frequency, $\omega=\omega_{0}-\gamma_{i}\beta/\gamma_{r}$.

In order to study the variability of the cycle amplitude, we need to account for the randomness inherent to the dynamo process. There is strong evidence that the dynamo is of Babcock-Leighton type (Wang & Sheeley, 2009; Muñoz-Jaramillo et al., 2013). This means that magnetic flux connected to the polar fields (axial dipole) is the relevant poloidal field for the generation of the toroidal field by differential rotation in the convection zone (Cameron & Schüssler, 2015). The polar fields result from the emergence of bipolar magnetic regions (sunspot groups) with a systematic average tilt relative to the solar East-West direction together with the subsequent transport of their magnetic flux across the surface by differential rotation, supergranulation, and meridional flow (Mackay & Yeates, 2012; Wang, 2016). Consequently, any scatter in the properties of the bipolar regions (e.g., emergence latitude, tilt angle) or in the flux transport process introduces a corresponding scatter in the resulting poloidal dipole field and, therefore, in the amplitude the subsequent cycle (Charbonneau & Dikpati, 2000; Wang & Sheeley, 2009; Jiang et al., 2014). Since Babcock-Leighton dynamos can be put into the general mathematical framework of $\alpha\Omega$-dynamos (Stix, 1974), we can describe this effect by random scatter in the term generating the poloidal field, which enters the normal form model in the form of multiplicative noise. Eq. (1) thus transforms into a stochastic differential equation, viz.

[TABLE]

where we take $W_{c}$ to represent a complex Wiener process (random walk with uncorrelated, Gaussian distributed increments) with a variance of unity after 11 years. The real parameter $\sigma$ then corresponds to the standard deviation of the cycle amplitudes due to the noise in the generation process.

Values for all five parameters entering Eq. (2) are constrained by empirical results. The linear growth rate, $\beta$, can be estimated from the time scale for the recovery from a grand minimum of very low activity. The duration of the Maunder minimum of about 70 years implies a recovery time of the order of decades, so that we take $\beta=1/50$ year*-1* as a reference value. The activity cycle appears to have persisted at low amplitude and with unchanged period during the Maunder minimum (Beer et al., 1998; Vaquero et al., 2015), so that we take the (magnetic) cycle frequency to be unaffected by the weak nonlinearity, $\omega=\omega_{0}=2\pi/22\,$year*-1*, and thus $\gamma_{i}=0$. In order to fix the parameter $\gamma_{r}$, we need to determine the (average) cycle amplitude, $\sqrt{\beta/\gamma_{r}}$, in terms of the sunspot number. For a cycle of sinusoidal shape, the amplitude is given by $\pi/2$ times its mean level. Using the average of the Group Sunspot Numbers (Hoyt & Schatten, 1998) since 1700, the end of the Maunder minimum, we obtain a cycle amplitude of 64 and thus $\gamma_{r}=4.9\cdot 10^{-6}$ year*-1*.

The remaining parameter to be fixed is the noise level, $\sigma$, which we relate to the scatter of the tilt angles of bipolar magnetic regions. Surface flux transport simulations show that the observed Gaussian scatter of the tilt angles leads to a fluctuation level of 30%–40% of the polar dipole field (Jiang et al., 2014; Muñoz-Jaramillo et al., 2013). This implies a degree of fluctuations of the cycle amplitudes owing to the linear correlation between the polar field at the end of a cycle and the amplitude of the next cycle (Wang & Sheeley, 2009). This level is consistent with the value of about 35% for the scatter of the cycle amplitudes determined from the sunspot numbers observed since 1700 amounts and also with the scatter shown by the 10-year sampled sunspot numbers reconstructed from the cosmogenic isotopes (Usoskin et al., 2016). Since the limit cycle is an attractor, we have a slight damping of the random perturbations during a cycle, so that we have to use a somewhat higher value, $\sigma=0.40$, in order to compensate for this effect and reproduce the observed cycle-to-cycle fluctuation of 35%. The effect on the results of varying the parameters $\beta$ and $\sigma$ is discussed further below.

We used the Euler-Maruyama method (Kloeden & Platen, 1992) to perform Monte-Carlo simulations of Eq. (2) with the parameter values given above and a time step of one day. Random numbers were generated with the routine ‘random.normal’ from the Python package ‘numpy’ version 1.10.4. The convergence of the numerical result was checked against the analytic solution of the normal form without noise. Covering 10,000 years for each realization, we obtained the amplitude variability over a large range of time scales. Results for one such realization are illustrated in the middle panels of Fig. 1. Panel C covers a period of 450 years, exhibiting extended periods of very low activity (grand minima). Panel D shows the 10-year averages for the whole time interval of 10,000 years. The main features of the actual and the simulated time series are qualitatively similar: there is considerable variability of the cycle amplitudes with occasional grand minima (and maxima). For a more quantitative comparison, Fig. 2 gives the non-cumulative (left panels) and cumulative (right panels) distributions of the lengths of grand minima and of the waiting times between the end of one grand minimum and the start of the next from the simulations (black lines) in comparison to the empirical distributions (blue lines) derived from the cosmogenic isotope record (Usoskin et al., 2016). The model distributions are well approximated by exponentials, which would be expected for a Poisson process (see also Olemskoy et al., 2013). Overall, the distributions appear to be consistent with each other. Whether the empirical events in the tail of the waiting time distribution represent a significant deviation from an exponential distribution (cf. Moss et al., 2008) cannot be decided owing to the small number (20) of grand minima in the isotope record. Note that the normal-form model generically shows ongoing low-amplitude cycles during grand minima, such as indicated during the Maunder minimum (Beer et al., 1998; Vaquero et al., 2015).

Further quantitative comparison between model and data is obtained by considering power spectra as shown in Fig. 3. Panel A corresponds to the historical sunspot record (yearly resolution, green line) and to the reconstruction based on cosmogenic isotopes (10-year resolution, blue line) in comparison to the median of 1000 realizations of 10,000 years each of the normal-form model (black line). The model results are consistent with the observed spectra, demonstrating that a weakly nonlinear limit cycle driven by random noise is sufficient to explain the observed variability of solar activity on a wide range of time scales. A single realization of the model is considered in panel B of Fig. 3, which shows spectra for 450 years with yearly resolution (orange line) and for 10,000 years with 10-year resolution (purple line). These spectra are qualitatively similar to the observed spectra shown in panel A, including some apparent (but spurious) long-term periodicities.

We also considered the effect of varying the parameters $\beta$ (linear growth rate, related to the supercriticality of the dynamo) and $\sigma$, the level of the random forcing. They represent two competing effects: the forcing tends to drive the solution away from the limit cycle while perturbations of the stable limit cycle decay with a timescale $(2\beta)^{-1}$. Changing these parameters therefore varies the relative importance of the two effects. In addition to the reference value, $\beta=1/50$ year*-1*, suggested by the recovery of the Sun from the Maunder minimum, we considered values representing a much smaller growth rate, $\beta=1/250$ year*-1*, and a markedly higher growth rate, $\beta=1/10$ year*-1*. Likewise, we took values for the random forcing of $\sigma=0.$ (no perturbations) and $\sigma=1.$ (very strong perturbations), in addition to the reference value of $\sigma=0.4$ indicated by the observed variability of the solar cycle maxima.

For the cases with $\beta=1/50$ year*-1* and $\beta=1/10$ year*-1* we performed 1000 simulations of 11,000 years length each. We excluded the first 1,000 years from each realization in order to stay clear of the transients related to the arbitrary initial conditions. For the very weakly excited cases with $\beta=1/250$ year*-1*, we run 1000 simulations of 21,000 years length each, excluding the first 11,000 years from the analysis. To illustrate the results, a stretch of 450 years of one arbitrarily chosen realization for each pair of parameters is shown in Fig. 4, while 10-year running averages covering 10,000 years are given in in Fig. 5. In the cases without random forcing (bottom rows of the figures), the solutions are rectified sine functions with amplitude $\sqrt{\beta/\gamma_{r}}$. For the observationally well constrained reference value of the random forcing, $\sigma=0.4$, the cycles show variations in amplitude and phase (middle rows in Figs. 4 and 5). The value of $\beta$ determines the rate of occurrence and the lengths of grand minima: since the recovery time from perturbations is $(2\beta)^{-1}$, low values of $\beta$ lead to long extended minima, while high values suppress their occurrence. In the case of very strong random forcing (top rows of Figs. 4 and 5), the perturbations lead to high cycle-to-cycle fluctuations and to very long grand minima for low $\beta$. From these results, we would expect that very slowly rotating stars, which are even nearer to critical dynamo excitation, show very long quiescent phases, while stars with higher excitation (faster rotators) display very strong cycle-to-cycle fluctuations and, if at all, only rarely show short grand minima of activity.

Fig. 3 shows median power spectra corresponding to 1000 realizations for each of the nine parameter combinations. The case without perturbations (bottom row) obviously shows the spectra of a rectified sine wave. For growing strength of the random forcing, the 11-year peak becomes less prominent and broader as the amplitude and phase fluctuations of the cycles grow stronger. For growing dynamo excitation (higher $\beta$), the long-period tail of the spectrum flattens as a result of the shortening and increasing suppression of the grand minima.

3 Babcock-Leighton-type dynamo model

To demonstrate that the normal-form approach actually fits in the context of a more detailed nonlinear dynamo model, we consider an updated version of the Babcock-Leighton dynamo (Babcock, 1961; Leighton, 1969; Wang et al., 1991), which reproduces key features of the solar cycle (Cameron & Schüssler, 2017). The model comprises the essential ingredients of a flux-transport dynamo (Charbonneau, 2010, 2014). It considers the axisymmetric part of the magnetic field and is based on the evolution equations for the azimuthal component of the vector potential (determining the poloidal field) at the solar surface,

[TABLE]

and the radially integrated toroidal flux per radian,

[TABLE]

both as functions of colatitude, $\theta$, and time. Here, $R_{\odot}$ is the solar radius, $R_{b}$ the radial location of the bottom the solar convection zone. $B_{r}$ and $B_{\phi}$ are the radial and azimuthal components, respectively, of the magnetic field. Introducing fluctuations in the source term for the poloidal field, the evolution equations given in Cameron & Schüssler (2017) become stochastic differential equations, viz.

[TABLE]

and

[TABLE]

Here $\Omega_{R_{\odot}}(\theta)$ and $\Omega_{R_{\rm NSSL}}(\theta)$ are, respectively, the angular rotation rates at the solar surface and at the base of the near-surface shear layer, for which we take the same profiles as in Cameron & Schüssler (2017). $\eta_{R_{\odot}}=350$ km${}^{2}\cdot$s*-1* and $\eta_{0}$ are, respectively, the magnetic diffusivities at the solar surface and in the bulk of the convection zone. $U(\theta)$ represents the poleward meridional flow at the surface, for which we use the profile given by Hathaway & Rightmire (2011, see their Eqs. 9–11). $V_{0}\sin(2\theta)$ refers to the equatorward return flow affecting the toroidal field in the bulk of the convection zone; we take $V_{0}=2$ m$\cdot$s*-1*, which roughly corresponds to the speed of equatorward propagation of the activity belts. $a_{\rm S}$ is the source of the poloidal field resulting from the emergence of bipolar magnetic regions, which we write as

[TABLE]

where we have introduced a nonlinearity with parameter $b_{c}$. The quantities $\alpha_{0}$, $\eta_{0}$, $b_{c}$ and $\sigma^{*}$ are the free parameters of the model.

The random forcing is considered as a two-dimensional Wiener process, $W(t,\theta)$, which depends on both latitude and time and has a variance of 1 radian*-1* after 11 years. The strength of the forcing is determined by the parameter $\sigma^{*}$. Since the noise parameter, $\sigma$, in the normal-form model is independent of the growth rate (dynamo excitation), for consistency we take $\sigma^{*}\alpha_{0}={\rm const.}$ to represent a fixed noise level. In the computations, we used a second-order centered difference scheme for the spatial derivatives with 180 grid points in colatitude, and advanced the solution in time using the Euler-Maruyama method with a time step of 1 day. The numerical results for mildly supercritical, nonlinear dynamo action with this model ($\eta_{0}=65$ km${}^{2}\cdot$s*-1*, $\sigma^{*}\alpha_{0}=0.17$ m$\cdot$s*-1*, $b_{c}=10^{24}$ Mx, and $\alpha_{0}=5.7\alpha_{\mathrm{crit}}=2.5$ m$\cdot$s*-1* where $\alpha_{\mathrm{crit}}$ is the critical value of $\alpha$ for the onset of dynamo action with the other parameters as stated), are shown in the bottom panels of Figs. 1 and 3. The measure of the activity shown here is the integrated subsurface toroidal field corresponding to the dipole mode of the dynamo $|\int_{0}^{90}b\mathrm{d}\theta-\int_{90}^{180}b\mathrm{d}\theta|$, scaled to a similar level as the sunspot number. They are consistent with the results from observations and from the generic noisy normal-form model.

We also carried out a parameter study with this dynamo model. There is no simple relationship between the four free parameters of the Babcock-Leighton model ($\alpha_{0}$, $\eta_{0}$, $b_{c}$ and $\sigma^{*}$) and the four parameters of the normal-form model ($\omega_{0}$, $\beta$, $\gamma_{r}$ and $\sigma$). We have chosen parameters so that the solutions have a period of about 22 years (i.e. an 11-year activity cycle) and kinematic growth rates and levels of noise similar to those of the cases presented in the previous section for the normal-form model. This specifies three of the four degrees of freedom. The fourth constraint is equivalent to specifying the amplitude of the limit cycle. In this illustrative study we have kept $b_{c}=10^{24}$ Mx and expect solutions of this order of magnitude – the exact amplitude of the limit cycle will however also depend on the other parameters in a non-trivial way.

The results of the parameter study are illustrated in Figs.7–9, which show the corresponding quantities in the same format as in Figs. 4–6 for the normal-form model. For cases with a very low kinematic growth rate of 1/250 year*-1* (left column of the figures) we took $\alpha_{0}=1.8\alpha_{\mathrm{crit}}=1.66$ m$\cdot$s*-1* and $\eta_{0}=75$ km${}^{2}\cdot$s*-1*. For the cases with a growth rate of 1/50.8 year*-1* (middle column) we chose $\alpha_{0}\bf=5.7\alpha_{\mathrm{crit}}=2.5$ m$\cdot$s*-1* and $\eta_{0}=65$ km${}^{2}\cdot$s*-1*. Finally, for the case with a high growth rate of 1/5.2 year*-1* (right column) we used $\alpha_{0}=167\alpha_{\mathrm{crit}}=20$ m$\cdot$s*-1* and $\eta_{0}=30$ km${}^{2}\cdot$s*-1*. For the purpose of a qualitative comparison, it is unnecessary to perform a tedious fine tuning of the parameters in order to exactly match the growth rates to those of the normal-form model. The cases without random forcing are shown in the bottom row of Figs. 7–9. The forcing for the cases given in the middle row reproduces the observed variability of the sunspot maxima since 1700 for the reference case (central panel of the figures). For the cases with very strong forcing (top row of the figures) we multiplied the reference value by the same factor 2.5 as in the case of the normal-form model. In all cases, we carried out 150 simulations covering a time of 11,000 years, omitting the first 1,000 years from the analysis in order to exclude initial transients.

Comparing the results for both models given in Figs. 4–6 (normal-form model) and Figs. 7–9 (dynamo model), respectively, we find that they are similar for the cases with low and medium growth rates, including the ”solar” reference models (central panels of the figures). Because we have varied the two of the model parameters for the dynamo cases with growth rates of 1/250 years, 1/50.8 years and 1/5.2 years the amplitudes changes between the different is difficult to interpret. While the individual realizations obviously differ in detail between the normal form and dynamo models, the average spectra are very similar. For the high growth rate, however, the models yield clearly different results. The corresponding growth time is shorter than the cycle period, which invalidates the normal-form approach. Furthermore, the Babcock-Leighton model has entered a strongly nonlinear and presumably chaotic regime. While this indicates the limits of that model, the similarity in the other cases demonstrates its validity for not too strong dynamo excitation. This is the realistic case for the Sun, but we expect more rapidly rotating and very active stars to be in a more strongly nonlinear or even chaotic regime (see, e.g., Tobias et al., 1995).

4 Conclusion

Our results suggest that the variability of the solar cycle amplitudes between decadal and millennial timescales can be understood in terms of a weakly nonlinear and noisy limit cycle. This approach is motivated by observational results and its parameters are constrained by observations as well. It represents the generic model for the fundamental mode of a weakly excited $\alpha\Omega$-dynamo, such as the observationally well supported Babcock-Leighton approach. Owing to its simplicity, the model does not cover the possible forcing of higher dynamo modes (such as the quadrupole mode leading to hemispheric asymmetry) and also does not account for the deviations of the solar cycle from a (rectified) sinusoidal shape. On the other hand, our results show that the long-term variability of solar activity is consistent with fluctuations due to a stochastic process, such as random scatter in the tilt angles of bipolar magnetic regions and sunspot groups. No intrinsic periodicities apart from the 11-year cycle are required to understand the variability, although the possible existence of such periodicities cannot be strictly excluded by our analysis.

The data for the Group Sunspot Numbers were obtained from the SILSO data base maintained by the Royal Observatory of Belgium, Brussels (http://sidc.oma.be/silso/groupnumberv3). I. Usoskin kindly provided the sunspot number reconstruction from the cosmogenic isotope record presented in Usoskin et al. (2016).