Shaken and stirred: When Bond meets Suess-de Vries and Gnevyshev-Ohl

We argue that the most prominent temporal features of the solar dynamo, in particular the Hale cycle, the Suess-de Vries cycle (associated with variations of the Gnevyshev-Ohl rule), Gleissberg-type cycles, and grand minima can be self-consistently explained by double synchronization with the 11.07-years periodic tidal forcing of the Venus-Earth-Jupiter system and the (mainly) 19.86-years periodic motion of the Sun around the barycenter of the solar system. In our numerical simulation, grand minima, and clusters thereof, emerge as intermittent and non-periodic events on millennial time scales, very similar to the series of Bond events which were observed throughout the Holocene and the last glacial period. If confirmed, such an intermittent transition to chaos would prevent any long-term prediction of solar activity, notwithstanding the fact that the shorter-term Hale and Suess-de Vries cycles are clocked by planetary motion.


Introduction
Thanks to the seminal work of Gerard Bond and his collaborators, we now have overwhelming evidence that a significant component of sub-Milankovich climate variability occurs in certain 1-3-kyears "cycles" of abrupt changes of the North Atlantic's surface hydrography (Bond et al., 1997(Bond et al., , 1999, and that those Bond events are closely related to corresponding variations in solar output, evidenced by measurements of the cosmogenic radionuclides 14 C and 10 Be (Bond F. Stefani f.stefani@hzdr.de 1 et al., 2001). While originally identified for the Holocene (Bond et al., 1997) from certain ice drift proxies (in particular volcanic glass from Iceland and hematitestained grains from East Greenland, found in deep-see sediments), many similar events were later revealed by Bond et al. (1999) in the last glacial period, too. Viewed from this angle, the little ice age, comprising in particular the Spörer and Maunder grand minima, appears just as the latest link in the chain of Bond events, and the temperature increase since the end of the Dalton minimum as a rebound from those frosty times.
In solar physics, similar variations with time scales of 1-3 kyears are usually discussed under the notion Eddy cycle and Hallstadt cycle Abreu et al., 2012;Soon et al., 2014;Scafetta et al., 2016;Usoskin et al., 2016). Yet, some caution seems to be appropriate when stretching the very concept of "cycles" from the decadal (Schwabe, Hale) to the millennial time scale, in particular when the underlying 14 C and 10 Be data bases have typical durations of only 10 kyears, or just slightly longer (Kudryavtsev and Dergachev, 2020). To gain more insight into the statistics of those "cycles", we make here the plausible assumption that the established link between solar activity and Bond events, with correlation coefficients of around 0.5 during the Holocene (Bond et al., 2001), extends also to the last glacial. With this proviso, we re-plot in Figure 1(a) the data of time separations (or waiting times) between the 54 Bond events as identified over the last 80 kyears, which we have drawn from Figure 6(c) of Bond et al. (1999). What we observe then is a broad range of time separations between 600 years and 2600 years, with a mean value of 1469 years and a standard deviation of 514 years, according to Bond et al. (1999). However, when focusing only on the first eight time separations in the Holocene, located chiefly around 1500 years, 2400 years and 600 years, it comes as no surprise to find similar periods in Fourier or wavelet analyses (Mayewski et al., 1997;Dima and Lohmann, 2009;Soon et al., 2014). Special attention on the mean 1470-years cycle, and even speculations about its origin from a coincidence of 17 Gleissberg and 7 Suess-de Vries cycles (Braun et al., 2005), seem to be justified in this case.
If, on the contrary, we consider the entire series of Bond events over 80 kyears, we get the same impression as Usoskin et al. (2007) who had argued that the "occurrence of grand minima/maxima is driven not by long-term cyclic variability, but by a stochastic/chaotic process." More quantitatively, the random walk character of this series is analyzed in Figure 1(b) which shows Dicke's ratio N i r 2 i / i (r i − r i−1 ) 2 between the mean square of the residuals r i (i.e., the distances between the actual Bond events and hypothetical events according to a linear fit to the series) to the mean square of the differences r i −r i−1 between two consecutive residuals (Dicke, 1978). Obviously, Dicke's ratio for N Bond events taken into account roughly approaches the theoretical random walk dependence (N + 1)(N 2 − 1)/(3(5N 2 + 6N − 3)), with its limit N/15, while significantly deviating from the corresponding dependence (N 2 − 1)/(2(N 2 + 2N + 3)) for a clocked process, with its limit 1/2, as it had been confirmed previously for the Schwabe cycle (Stefani, Giesecke and Weier, 2019). Figure 1(a) is complemented by another curve of 16 time separations (restricted to an interval of 30 kyears), as it will come out of our numerical model to be described further below. Suffice it to say here that the average of the  Bond et al. (1999) for the data of the deep-sea sediment core VM23-81. Green empty squares denote the corresponding time separations of numerically simulated "Bond events" for the parameters ω 0 = 10000, α c 0 = 15, α p 0 = 50 q p α = 0.2, q c α = 0.8, D = 0.05 and κ(t) = 0.6 + 0.2m(t)). (b) Dependence of Dicke's ratio on N for the two time series from (a), together with the theoretical curves for a random walk and a clocked process. Note that the horizontal positions of the individual points slightly differ between (a) and (b) since the abscissa in (a) shows the centre point for the time separation, while the abscissa in (b) indicates the number of Bond events taken into account.
waiting times, and their broad distribution, are quite similar to those of the 54 real Bond events, and that Dicke's ratio points also in direction of a random walk process, although the statistical significance of those mere 17 numerical "Bond events" is certainly not satisfactory.
In this paper we will try to give a relatively simple and self-consistent answer to the questions of how such an intermittent and random walk like behaviour comes about, and how it can be related to the much clearer periodicities of the Hale, the Suess-de Vries, and the Gleissberg cycle(s). For this purpose, we will start from a rather conventional α − Ω-dynamo model in form of a simple 1D PDE system (with co-latitude as the only spatial variable), to allow for very long simulations of appr. 30 kyears. With the source terms α and Ω appropriately chosen, this dynamo develops typical oscillation periods of 20-30 years. After enhancing these internal "stirring" terms of the dynamo by two external "shaking" terms with the specific periods 11.07 years and 19.86 years (both being related to planetary motions), we will end up with a more or less complete reproduction of the most relevant temporal features of the solar dynamo. But beforehand, we have to clarify the origin of these two external "shaking" terms, and how their use in a solar dynamo model can be justified.
We start with the 11.07-years period. Although the similarity of the Schwabe cycle and the spring-tide period of the Venus-Earth-Jupiter (VEJ) system has been known for a long time (Bollinger, 1952;Takahashi, 1968;Wood, 1972;Condon and Schmidt, 1975), the precise correspondence of these two 11.07-years periodicities was recognized only recently by Hung (2007); Scafetta (2012); Wilson (2013); Okhlopkov (2014Okhlopkov ( , 2016. According to Scafetta (2012), this period, corresponding to 4043 days, results from the resonance condition for the period of VEJ alignments, with the sidereal periods P V = 224.701 days, P E = 365.256 days, and P J = 44332.589 days for Venus, Earth and Jupiter, respectively. Quite often, the action of planetary tidal forces on the Sun is discarded in view of the tiny acceleration in the order of 10 −10 ms −2 (De Jager and Versteegh, 2005;Callebaut, de Jager, and Duhau, 2012), leading to a negligible tidal height h tidal ≈ GmR 2 tacho /(g tacho d 3 ) in the order of 1 mm (m is the planet's mass and d its distance to the Sun). One should note, however, that this tidal height translates (by virtue of the virial theorem) into a non-negligible tidal flow velocity of v ∼ (2g tacho h tidal ) 1/2 ≈ 1 m/s, when taking into account the huge gravity at the tachocline of g tacho ≈ 500 m/s 2 (Öpik, 1972). But even then it is hard to conceptualize how tidal forces could influence the solar dynamo without employing any sort of amplification mechanism. One candidate for such a mechanism was discussed by Wolff and Patrone (2010); Scafetta (2012) who speculated that planetary forces might affect the nuclear burning rate deep in the solar core and that this effect would be promptly felt at the tachocline via the resulting change in g-mode oscillations. While this sounds highly speculative, we mention here a pertinent observation by Kotov and Haneychuk (2020) regarding a possible influence of Jupiter's synodic cycle on the (still controversially discussed) 160 min oscillation of the Sun's radius.
An entire suite of most interesting synchronization models was corroborated by Wilson (2013). First, the 11.07-years periodicity was explained in terms of a VEJ tidal-torque model, in which Jupiter plays a distinguished role by exerting a torque upon the periodically forming Venus-Earth tidal bulges. Second, a gear effect was invoked to modulate the changes in the rotation rates (as driven by the tidal-torque effect) by the 19.86-year periodic Saturn-Jupiter quadratures, leading ultimately to a long-term modulation with 193-year period (which will play a key role further below). Further derivations of a 208-years Suess-de Vries cycle, a 1156-years Eddy cycle, and a 2302-years Hallstatt cycle can also be found in this highly instructive and recommendable paper.
In Wilson (2013) and in the earlier model of Zaqarashvili (1997), the synchronization of the solar cycle was essentially based on some variation of the rotation rate in the tachocline and/or the convective layers of the Sun. This leads to the questions: how can the solar dynamo be synchronized by such a weak variation of Ω? Both from our above estimation of the tidal effect on the velocity, as well as from the helioseismological measurements (Howe, 2009), we can infer that the variation of the Ω effect is certainly not larger than 1 per cent, and very likely much smaller, since a significant portion of the Ω variation can also be attributed to the back-reaction of the self-excited field on the flow (Proctor, 2007). If such a minor change is not sufficient to entrain the entire solar dynamo, can we perhaps take resort to the idea of Abreu et al. (2012) who emphasized that the maximum field strength of flux tubes (which can be stored prior to eruption) is very sensitive to small perturbations by gravitational, tidal and, as we add here, centrifugal forces due to changes of Ω? Although our preliminary efforts  to implement synchronization mechanisms of this sort into a Babcock-Leighton type dynamo model with time delay (Wilmot-Smith et al., 2006), either via a direct 11-07-years variation of the Ω effect or a corresponding variation of the rise condition for flux tubes, have not been very promising so far, we do not exclude that greater success may result from future investigations.
Still another avenue for synchronization was opened by Weber et al. (2015); Stefani et al. (2016), who recognized that the intrinsic helicity oscillations of the current-driven, kink-type Tayler instability (Tayler, 1973;Pitts and Tayler, 1985;Seilmayer et al., 2012;Weber et al., 2013), characterized by an azimuthal wavenumber m = 1, can be entrained by a tide-like (m = 2) perturbation, without (or barely) changing the energy content of the m = 1 mode. This idea of an "energy-efficient" mechanism of helicity synchronization was then first incorporated into a simple ODE solar dynamo model (Stefani et al., 2016(Stefani et al., , 2017, and later into a 1D PDE system with the co-latitude as the only spatial variable (Stefani, Giesecke and Weier, 2019). From the 11.07-years tidal entrainment of the helicity, and the α-effect associated with it, these models produced dipolar fields with an oscillatory 22.14-years Hale cycle, although in some parameter regions quadrupolar and hemispherical fields were observed, too. For the somewhat academic case of a purely 11.07-years periodic α-effect we proved the existence of a massively nonlinear dynamo of the Tayler-Spruit type (Spruit, 2002;Rüdiger and Schultz, 2020), whereas for a more realistic hybrid dynamo, comprising both an externally "shaken" and an internally "stirred" α-term, synchronization was accomplished by parametric resonance.
We would like to point out that meanwhile the empirical evidence for an 11.07-years synchronization is quite impressive, though not accepted (or not even recognized) throughout the solar dynamo community. Not only have the cycle minima from the last 1000 years turned out to be very close to a clockedprocess with 11.07-years periodicity (Stefani, Giesecke and Weier, 2019), but various algae growth data from a 1000-years period in the early Holocene have also shown a phase coherent cycle with basically the same period (Vos et al., 2004).
As longer-term cycles are concerned, it was recently confirmed (Stefani et al., 2020a) that the modulation period of the duration of the Schwabe cycles, as inferred from Schove's maxima data (Schove, 1983), is close to 200 years, a number which is consistent with previous results for the Suess-de Vries cycle relying on historic sunspot observations (Ma and Vaquero, 2020), 10 Be and 14 C data (Muscheler et al., 2007), and various climate related data (Lüdecke, Weiss and Hempelmann, 2015). It was not least the relative sharpness of that Suessde Vries cycle which had motivated many authors (Jose, 1965;Fairbridge and Shirley, 1987;Charvatova, 1997;Landscheidt, 1999;Abreu et al., 2012;Wolff and Patrone, 2010;McCracken, Beer and Steinhilber, 2014;Cionco and Soon, 2015;Scafetta et al., 2016) to search for a link between the solar dynamo and planetary forcings with correspondingly long periods.
Yet, when entering this playing field one can hardly avoid the question why not to consider, first and foremost, the strongest of all planetary influences on the Sun's motion, namely the 19.86-years synodic cycle of Jupiter and Saturn. This cycle governs the orbit of the Sun around the barycenter of the planetary system, comprising vast deflections in the order of the Sun's diameter and velocities of up to 15 m/s (Sharp, 2013;Cionco and Pavlov, 2018). Superposed on that period are minor wiggles stemming mainly from the orbits of Uranus and Neptun, which ultimately leads to a rather complicated motion with another 172-years periodicity, sometimes called "Jose cycle" (Jose, 1965;Charvatova, 1997;Landscheidt, 1999;Sharp, 2013). Still, it is the dominant 19.86-years cycle which has the capacity to produce, in concert with the 22.14-years Hale cycle, a beat period of 19.86 × 22.14/(22.14 − 19.86) = 193 years, as worked out in the above mentioned paper by Wilson (2013) and also by Solheim (2013). This 193-years period is suspiciously close to the Suess-de Vries cycle. Assuming an appropriate, though not yet well understood, coupling effect of the Sun's orbital motion into some change of the stratification of the tachocline, this beat period was numerically found to produce a 193-years modulation of the North-South asymmetry of the dynamo field (Stefani et al., 2020a).
This sequel to Stefani, Giesecke and Weier (2019); Stefani et al. (2020a) is structured as follows: in Section 2 we recapitulate our 1D solar dynamo code, and motivate the choice and structure of its main ingredients, such as the Ω effect, two parts of the α-effect, and the time variation of the loss parameter κ. The results of our simulations over 30 kyears are then presented in Section 3. We will illustrate how a weak time-variation of κ produces a clear 193-years modulation of the North-South-asymmetry, and the Gnevyshev-Ohl effect associated with it. For stronger variations of κ we then obtain intermittent breakdowns of the solar cycle which are indeed reminiscent of grand minima, and clusters thereof. Since those breakdowns occur already in the absence of any noise, we argue here for the onset of an intermittent route to chaos. One of those long runs will be utilized to produce the numerical time series of breakdowns as shown in Figure 1. We will conclude with a number of suggestions for future work, including an urgent call for a better physical and numerical modelling of the two main synchronization mechanisms, which in the present paper can only be implemented in a parameterized form.

Numerical model
Inspired by early work of Parker (1955); Schmalz and Stix (1991); Jennings and Weiss (1991); Roald and Thomas (1997); Kuzanyan and Sokoloff (1997), we will use here the same system of PDSs as in Stefani, Giesecke and Weier (2019); Stefani et al. (2020a), with the solar co-latitude θ as the only spatial coordinate. While appropriately constructed ODE systems have been surprisingly successful in describing different types of dynamo modulations (Knobloch, Tobias and Weiss, 1998), and even supermodulation (Weiss and Tobias, 2016), we consider a 1D PDE model a reasonable intermediate step towards a 2D model (in r and θ) which is presently under development, but which might become numerically costly when aiming at very long dynamo runs over 30 kyears, say.
With the axisymmetric magnetic field split into a poloidal component B P = ∇ × (Ae φ ) and a toroidal component B T = Be φ , the 1D PDE system reads where A(θ, t) represents the vector potential of the poloidal field at co-latitude θ (running between 0 and π) and time t, and B(θ, t) the corresponding toroidal field. The two sources of dynamo action are the helical turbulence parameter α and the radial derivative ω = sin(θ)d(Ωr)/dr of the rotational profile. Here, α and ω denote the non-dimensionalized versions of the dimensional quantities α dim and ω dim , according to α = α dim R/η and ω = ω dim R 2 /η, where R is the radius of the considered dynamo region and η the magnetic diffusivity. The time is non-dimensionalized by the diffusion time, i.e. t = t dim η/R 2 . The boundary conditions at the North and South pole are A(0, t) = A(π, t) = B(0, t) = B(π, t) = 0. The PDE system is solved by a finite-difference scheme using the Adams-Bashforth method. The initial conditions are A(θ, 0) = 0 and B(θ, 0) = s sin(θ) + u sin(2θ), with the chosen pre-factors s = −1 and u − 0.001 denoting symmetric and asymmetric components of the toroidal field 1 . We employ the typical solar θ-dependence of the ω-effect (Charbonneau, 2010) in the form with a plausible, but still moderate, value ω 0 = 10000. The helical source term α comprises, first, a non-periodic part with a constant α c 0 and a noise term ξ(t) (which is defined by the correlator ξ(t)ξ(t + t 1 ) = D 2 (1 − |t 1 |/t corr )Θ(1 − |t 1 |/t corr )), and second, a periodic part where the B-dependent term has the typical resonance-type structure ∼ B 2 /(1+ q p α B 4 ). The expression t 11.07 denotes the dimensionless counterpart of the 11.07year tidal forcing period. With our special choice of the diffusion time R 2 /η = 110.7 years, this amounts to t 11.07 = 0.1. Note that the latitudinal dependence comprises the same smoothing term (although more conveniently written here) as in Stefani, Giesecke and Weier (2019), which avoids a steep jump of α at the equator.
The term κ(t)B 3 (θ, t) in Eq. (1), as originally introduced by Jones (1983); Jennings and Weiss (1991), has been included to account for field losses owing to magnetic buoyancy, on the assumption that the escape velocity is proportional to B 2 . While we openly admit that the spin-orbit coupling of the angular momentum of the Sun around the barycenter into some dynamo relevant parameters remains an open question (for ideas, see Zaqarashvili (1997); Palus et al. (2000); Wilson (2008); Sharp (2013)) we employ in the following a time-variation of the parameter κ(t) proportional to the time series of the angular momentum. Since κ(t) is related to the very sensitive adiabaticity in the tachocline (Abreu et al., 2012), which could be easily influenced by slight changes in the internal rotation profile, its modification by some sort of spin-orbit coupling seems not completely unrealistic.
For the computation of the Sun's orbital angular momentum we utilized the DE431 ephemerides (Folkner et al., 2014) in the time interval between 13199 B.C.-A.D. 17000, of which a ≈1800-years segment is visualized in Figure 2(a). This function is dominated by the 19.86-years synodes of Jupiter and Saturn, to which further contributions, mainly from Uranus and Neptun, are added (see the PSD in Figure 2(c)). Further below, we will use the normalized version m(t) of this angular momentum curve for parametrizing the time-variation of κ(t). For the sake of comparison, we will also assess the results for a modified variant m JS (t) which relies exclusively on the 19.86-years periodic motion of Jupiter and Saturn.

Results
In this section, we present and discuss numerical results for a sequence of parameters similar to those utilized in Stefani et al. (2020a). The main difference is the much longer simulation time of 30 kyears, which is actually the interval for which the orbital angular momentum of the Sun was available to us (Folkner et al., 2014). Such longer simulations will allow for a more systematic study of the typical breakdowns of the 193-years modulated wave, and some preliminary comparisons with the sequence of Bond events. Another difference to Stefani et al. (2020a) is that here we focus chiefly on the noise-free case in order to evidence the intermittent route to deterministic chaos. A few results with noise included will nevertheless be presented at the end of the section and in the Appendix.  -d), whose "patchy" appearances simply result from the large number (∼ 1350) of Hale cycles involved (more details will be shown further below). The right panels (e-h) show the associated power spectral densities (PSD) resulting from Lomb-Scargle analyses of the respective curves in (a-d).

All planets, no noise
The fixed parameters ω 0 = 10000, α c 0 = 15, α p 0 = 50 q p α = 0.2, q c α = 0.8, D = 0 are chosen similar to those of previous studies (Stefani, Giesecke and Weier, 2019;Stefani et al., 2020a), but note the complete absence of noise in the present runs. Moreover, here we set out with a sufficiently large value of α p 0 so that the dynamo is already synchronized to 22.14 years; the parametric resonance phenomenon behind the frequency synchronization, when going over from α p 0 = 0 to some finite value, had been discussed in detail in Stefani, Giesecke and Weier (2019); Stefani et al. (2020a) and will not be re-iterated here. Hence, the only parameter to be varied from top to bottom of Figure 3 is the loss parameter κ.
We start in Figure 3(a,e) with a time-independent value κ(t) = 0.6, which yields a very clean Hale cycle with a period 22.14 years. Already in panel (a) we observe a slight asymmetry between positive values (reaching ≈ 9) and negative values (reaching ≈ −8) that can be attributed to the presence of a mixed mode, in which a weak quadrupolar field component part is added to the dominant dipolar one. This effect is related to the Gnevyshev-Ohl rule (Gnevyshev and Ohl, 1948). Apart from the two minor peaks at the doubled and tripled Hale (b,f): κ(t) = 0.6 + 0.3m(t), as (a,e), but with a modulation of κ with an angular momentum function m(t) according to Figure 2(b). As seen in the spectrum (f), this dipole solution contains a beat period of 193 years, which also appears in (b) as a minor wiggle of B(72 • , t). Note the appearance of Gleissberg-type cycles around 100 years in (f). (c,g): κ(t) = 0.6 + 0.32m(t), as (b,f), but with a slightly increased variation of κ. Note the appearance in (c) of four sudden events, where the positive-negative asymmetry of B(72 • , t) changes sign. The spectrum (g) has become noisy. (d,h): κ(t) = 0.6 + 0.4m(t), as (c,g), but with increased variation of κ, producing significantly more breakdowns.
frequency (which naturally result from the nonlinear terms in the PDE system) the spectrum in (e) is quite smooth and featureless. This is changing in Figure 3(b,f) when the loss parameter is equipped with a time-variation according to κ(t) = 0.6 + 0.3m(t), where m(t) is the normalized orbital angular momentum function from Figure 2(b). The most prominent feature that appears in the PSD (f) is the strong peak at 193 years. In panel (b) this peak manifests itself in form of minor wiggles of the maxima and minima, which indeed correspond to a modulation of the positive-negative asymmetry. We note in passing the occurrence of some peaks at Gleissberg-type periods (around 96 years and 64 years), and two side bands at around 20 and 25 years which are reminiscent of the Wilson gap (Hathaway, 2010). These two features were discussed in more detail in Stefani et al. (2020a), and will not play a particular role in the following.
When increasing the variation of the loss parameter just a little further to κ(t) = 0.6 + 0.32m(t), we observe in panel (c) four sudden jumps of the positivenegative asymmetry. While the main peaks at 22.14 years and 193 years survive, the rest of the spectrum (g) becomes noisy, an effect of spectral leakage due to the segmentation of the entire time domain into 5 parts. While in this case one might still speculate about a regular occurrence of these breakdowns (with a waiting time of appr. 7 kyears), any alleged regularity is clearly lost in our last example (d,h), corresponding to κ(t) = 0.6 + 0.4m(t). Here we observe approximately 11 breakdowns without any clear periodicity. This insinuates an intermittent route to chaos, although we leave the detailed analysis of this transition, for our non-trivial case of double synchronized system, to future work.
For three of the examples from Figure 3, some details are discussed in the following. Figure 4 shows the central interval between 14 and 16 kyears from Figure 3(a), both for B(72 • , t) (a) and for the entire field B(θ, t) as shown here as a contour-plot (b). In panel (a) we observe again the positive-negative asymmetry (the value varies between -8 and +9), which translates into a North-South asymmetry as visible in (b) by the color asymmetry between Northern and Southern regions. Evidently this asymmetry is also connected with a Gnevyshev-Ohl rule. A clear 193-years modulation of the positive-negative asymmetry, and the Gnevyshev-Ohl effect related to it, appears in Figure 5 which illustrates some details of Figure 3(b). For the central interval of Figure 3(d), with κ = 0.6+0.4m(t), Figure 6 shows one typical breakdown of the modulated wave. The resulting disordered state, lasting appr. 500 years, resembles indeed a cluster of grand minima where the dynamo is not switched off (Beer, Tobias and Weiss, 1998)  in another state. Quite interesting here is the appearance of hemispherical fields (around 14900) and quadrupole fields (around 15000), which are reminiscent of sunspot observations within or shortly after the Maunder minimum (Sokoloff and Nesme-Ribes, 1994;Arlt, 2009). The corresponding trajectory in the dipolequadrupole space, as shown in Figure 7, resembles strongly the corresponding behaviour in Figure 5a of Knobloch, Tobias and Weiss (1998) and Figure 4 of Weiss and Tobias (2016). Meanwhile, a similar supermodulation effect has a also been found in a 3D simulation of dynamo action in rotating anelastic convection (Raynaud and Tobias, 2016).

Only Jupiter and Saturn, no noise
We now consider a modification of the angular momentum that enters the timevariation of the loss parameter κ(t). While in the previous subsection the full (normalized) angular momentum curve m(t) was used (violet full line in Figure 2(b)), we consider now the idealized curve m JS (t) as it would result from exclusively taking into account the orbital motion of Jupiter and Saturn (green dashed curve in Figure 2(b)). Figure 8 shows the corresponding numerical results for otherwise the same parameters as used previously in Figure 3. Superficially, the results are very similar, with the most significant differences showing up for the range of Gleissberg-type periods. With the simplified m JS (t)curve, we observe now clean peaks at one half (96.5 years) and one third (64.3 years) of the 193-years beat period, whereas in Figure 3 those peaks were more complicated. In Figure 9 we summarize our present understanding regarding the origin of the different peaks. The PSD for m(t) (violet) is a reproduction from Figure 2(c). It is dominated by the Jupiter-Saturn-peak at 19.86 years, and some other peaks, including the Jupiter-Neptune synode (12.78 years) and the Jupiter-Uranus synode (13.81 years). Both for the field resulting from m(t) and from m JS (t), the two dominant peaks are the Hale period (22.14 years) and the Suess-de Vries period (193 years). While the (blue) field curve for m JS (t) contains basically only one half (96.5 years) and one third (64.3 years) of this Suess-de Vries period, the (green) curve for the full m(t) contains additional peaks in the Gleissberg-region, comprising in particular the peaks at 55.8 years and 82.7 years which are beat periods between the 11.07-years Schwabe cycle and the Jupiter-Uranus synode (13.81 years) and Jupiter-Neptune synode (12.78 years), respectively. Moreover, a few additional peaks, indicated by question marks, seem to be related, e.g., to the Saturn-Neptun (35.87 years) and the Saturn-Uranus (44.36 years) synodes.  Figure 2(b)) which is restricted to the 19.86-years periodic part resulting from the orbital motion of Jupiter and Saturn only.

All planets, noise included
In Figure 10 we assess the role of noise. The parameters are identical to those in Figure 3, except that we use now a finite noise level D = 0.05. Not surprisingly, even without any κ variation (a), the spectrum (e) is already noisy, while the positive-negative asymmetry in (a) is very similar to Figure 3(a). Apart from that, the overall structure turns out to be quite comparable to Figure 3.
If we go beyond the examples of Figure 10 by choosing a still stronger variation κ = 0.6 + 0.5m(t), we end up with Figure 11 which shows now a longer segment between 4-16 kyears. While for such long simulations many details are lost in the contour-plot (b), it highlights the fact that the breakdowns occur at instants where the North-South asymmetry (evidenced in particular by the reddish parts) has acquired a certain critical threshold. This entire behaviour is reminiscent of the supermodulation described by Weiss and Tobias (2016); Raynaud and Tobias (2016).
With Figure 11(c) we add a contour-plot for the full 30-kyears simulation period, just to illustrate the sequence of 17 breakdowns which were used for the numerical curve in Figure 1 (where the time direction is inverted, though). As seen in Figure 1(b), those 17 events seem to obey the same random walk law as the 54 Bond events, although our 30-kyears simulation time is still too Figure 9. Comparison of the Lomb-Scargle PSDs for the angular momentum function m(t) (violet), for B(72 • , t) (green) resulting from κ(t) = 0.6 + 0.3m(t) with the full m(t), and for B(72 • , t) (blue) resulting from κ(t) = 0.6 + 0.2m JS (t) with the reduced m JS (t). The individual peaks of the PSD for m(t) are the same as in Figure 2(c). The Suess-de Vries period of 193 years, as well Gleissberg-type periods 193/2 and 193/3 years emerge already when using only m JS (t). Some more peaks, which appear when the full m(t) is utilized, can be attributed to corresponding peaks in m(t) (see the question marks). However, there are two additional peaks at 55.8 years and 82.7 years which represent beat periods between the 11.07-years Schwabe cycle and the Jupiter-Uranus synode (13.81 years) and the Jupiter-Neptune synode (12.78 years), respectively. WG denotes the Wilson gap between 19.86 and 25 years.
short for a convincing statistics. An interesting feature becomes visible in the wavelet spectrogram of Figure 11(c) for B(72 • , t), which shows the distribution of oscillations with period τ in the vicinity of time t. During the breakdowns, characterized by a reduced energy in the 22.14-years period, we observe a significant increase of the energy in a wide range of τ from 30 to 1000-years. In fact the spectrum seems to become continuous here which reflects a transition to chaos due to nonlinearity. The quantitative difference between "regular" and "chaotic" regimes is highlighted in Figure 12 which shows the wavelet spectral density integrated over the corresponding intervals, separated by the green dashed lines in Figure 11(d). With the energy in the τ -range from 30 to 1000-years being generally increased in the "chaotic" segments by 1-2 orders of magnitude, the particular peak around 200-years is still markedly pronounced. A corresponding behaviour had been reported in Figure 6 of McCracken, Beer and Steinhilber (2013).
In the Appendix, we will carry out a similar analysis as in Figure 11 but for the case of fixed κ = 0.6 without any time variation but stronger noise with D = 0.1. We will then also find breakdown regions but no particular role of the Suess-de Vries cycle. Figure 10. Same as Figure 3 but with noise intensity D = 0.05. Even without any κ variation, the spectrum (e) is already noisy, while the positive-negative asymmetry in (a) is very similar to that in 3(a). Apart from that, the overall structure, including the critical value for the transition to chaos, is very similar as for D = 0.

Summary and open problems
In this paper we have pursued the ambitious program of finding a more or less complete and self-consistent description of the most significant periodicities (and time-scales) of the solar dynamo. First, we recapitulated the basic idea that the 22.14-years Hale cycle is synchronized by the 11.07-years tidal (m = 2) forcing period of the tidally dominant VEJ-system, the importance of which had been pointed out earlier by Hung (2007); Scafetta (2012); Wilson (2013);Okhlopkov (2014Okhlopkov ( , 2016. In various numerical models (Stefani et al., 2016Stefani, Giesecke and Weier, 2019) this (weak) tidal forcing was supposed to trigger a resonant excitation of 11.07-years oscillations of that part of the α-effect which is related to a typical m = 1 instability within or close to the tachocline, such as the Tayler instability or a magneto-Rossby wave. Strong empirical evidence for a phase coherent 11.07-year Schwabe cycle comes both from algae data in the early Holocene (Vos et al., 2004), and from 14 C and 10 B data for the last 600 years (Stefani et al., 2020b). We note in passing that the numerical models produce, via the resonant dependence of the α-effect on the magnetic field strength, also secondary peaks of solar activity, which might be linked to mid-term oscillations (Obridko and Shelting, 2007;Valdés-Galicia and Velasco, 2008; McIntosh et al., 2015; Bazilevskaya et al., 2016;Karak, Mandal and Banarjee, 2018;Frick et al., 2020).
Second, motivated by ideas of Cole (1973); Wilson (2013); Solheim (2013), we argued that the Suess-de Vries cycle emerges as a 193-years beat period between the primary, tidally synchronized 22.14-years Hale cycle and the single strongest component of solar motion around the barycenter of the planetary system, which is governed by the 19.86-years synodic cycle of Jupiter and Saturn. Without a detailed model for spin-orbit coupling at hand, we hypothesized that this coupling would lead to a periodic variation of the field loss parameter κ in the tachocline region. Numerically, the combination of such a κ-variation with the synchronized component of the α-effect led to a modulation of the dynamo wave Figure 12. Comparison of the wavelet power density for the "chaotic" intervals (separated by the green dashed lines in Figure 11(d)) with that for the remaining "regular" intervals with a beat period of 193 years, which manifests itself in a modulation of the North-South asymmetry and, closely related to that, in a change of the dipolequadrupole relation (Knobloch, Tobias and Weiss, 1998;Moss and Sokoloff, 2017) and the Gnevyshev-Ohl rule. We would like to point out that the emergence of this beat period depends critically on the phase stability of the two underlying 11.07-years and 19.86-years processes, or, to put it otherwise: the existence of the long-term Suess-de Vries cycle gives a "backward argument" for the synchronized character of the short-term Hale cycle. As discussed in detail in Stefani et al. (2020a) (but only shortly touched upon in this paper), a stronger variation of κ leads to the occurrence of a Wilson gap (Wilson, 1987) with two side peaks at 19.86 years and 25 years. Some aspects of this behaviour were illustrated in Figure 9.
The main focus of this paper was, however, on the intermittent occurrence of grand minima, and clusters thereof. We have observed irregular breakdowns of the 193-years modulated dynamo wave, preferably at instants where the North-South asymmetry reaches a certain critical level. This threshold effect fits well to a corresponding observation of Tlatov (2013) for the Maunder and Dalton minima. Since those irregular breakdowns are already observed in the absence of any noise, they seem to be connected with an intermittent transition to (deterministic) chaos, similar as in the supermodulation concept developed by Weiss and Tobias (2016). For an appropriately chosen time-variation of κ (and some weak noise), we obtained a waiting time distribution with a similar mean value and standard deviation as inferred from the 54 Bond cycles observed over the last 80 kyears. Such an intermittent transition to chaos would hamper any long-term predictability of solar activity (and the climatic changes connected with it), even if the planetary clocking of the shorter-term Hale and Suess-de Vries cycles could be confirmed.
Based on a conventional α − Ω-dynamo, our model thus required only the two synchronization periods 11.07 years and 19.86 years, related to tidal forcing and the strongest component of solar orbital motion, respectively, to produce essentially all relevant periods, and "periods", of the solar dynamo. While this appears promising, we conclude with a discussion of remaining problems and "missing links" in the theory.
First, we have to admit that, up to present, the basic synchronization mechanism for the helicity of an m = 1 mode by some tide-like (m = 2) forcing has been evidenced only for the paradigmatic case of a non-rotating, full cylinder (Stefani et al., 2016). Preliminary attempts for a hollow (more "tachoclinic") cylinder were promising, although not entirely conclusive. Neither rotation nor stratification were implemented yet. A test of the same concept in the simplified framework of an m = 1 buoyancy instability of toroidal flux rings, as developed by Ferriz Mas, Schmitt, and Schüssler (1994), could be very helpful and instructive in this respect. Interestingly, helicity synchronization with an m = 2 forcing was numerically observed for the physically different, but topologically similar case of an m = 1 large scale circulation of Rayleigh-Bénard convection in a cylinder (Galindo, 2020). A liquid metal experiment to confirm this effect is presently under way (Stepanov and Stefani, 2019;Jüstel et al., 2020).
As this suggests a generic and robust character of the helicity synchronization mechanism, there is good hope to apply the same principle also to the m = 1 magneto-Rossby waves which were recently discussed by Dikpati et al. (2017); Marquez-Artavia, Jones, and Tobias, (2017);McIntosh et al. (2017);Zaqarashvili (2018). Unstable shallow-water modes had been shown earlier (Dikpati et al., 2009) to produce kinetic helicity, which is concentrated in the neighborhood of toroidal flux bands and migrates with them toward the equator as the solar cycle progresses. It should not be too complicated to implement a tide-like m = 2 forcing into those shallow-water models in 2D (here: in θ and longitude φ) with the aim to identify a similar helicity synchronization mechanism as found for the Tayler instability.
As a next step one could think about combining such a 2D model (in θ, φ) with another 2D model (in θ, r) as it was utilized in various Babcock-Leighton type dynamo models (Guerrero and de Gouveia Dal Pino, 2007;Jouve et al., 2008). Presently we are working on an enhancement of the latter type of models in order to assess the direction of the butterfly diagram and to see whether a weak synchronized part of α (with an amplitude of less then 1 m/s, as limited by the argument ofÖpik (1972)) is indeed sufficient for synchronizing the dynamo.
Turning to the spin-orbit coupling based on the 19.86-years orbital motion, and the resulting 193-years beat period, we first have to ask ourselves: does this beat period indeed correspond to the Suess-de Vries cycle? Since typical periods of this cycle between 190 until 210 years have been discussed in the literature, 193 years sounds not too bad in this respect. It also seems to be supported by a recent result of Ma and Vaquero (2020) who had found a 195-years period for the strong and clearly expressed Suess-de Vries cycle in the relatively "quiet" interval between 800 and 1340 A.D.
But even if the equivalence of the 193-years modulation with the Suess-de Vries could once be confirmed, we would still be left with the problem to explain the coupling of the (mainly) 19.86-years periodic orbit into some internal, dynamo relevant motion. To the best of our knowledge, there have been no serious attempts to tackle this problem in its full beauty, including a realistic orbital motion, the 7 degree inclination of the Sun's rotation axis, and an alleged non-sphericity of the tachocline with a prolateness that might even vary with the magnetic field strength (Dikpati and Gilman, 2001). In this respect, we recall a relatively recent result on the similar problem of precession, for which a significant braking of the (solid body like) rotational profile was observed already for weak precessional forcing (Giesecke et al., , 2019Meunier and Albrecht, 2020). If such a braking effect would also occur in the solar spin-orbit coupling problem, it could indeed result in a change of the very sensitive adiabaticity (Abreu et al., 2012), and the associated κ parameter as used here. Admittedly, this complex problem has to be left for future studies.
With the main focus laid on the synchronized helicity of m = 1 instabilities, we should not overlook alternative synchronization mechanisms based on axisymmetric (m = 0) instabilities, the relevance of which had been discussed by several authors (Dikpati et al., 2009;Rogers, 2011). The recent detection of a double-diffusive helical magnetorotational instability for flows with positive shear (Mamatsashvili et al., 2019) might be an interesting candidate in this respect, as it depends quite sensitively on the ratio of toroidal to poloidal field which, in turn, could be easily influenced by variations of κ.
In summary, given its skill in reproducing the various solar cycles, and "cycles", with only mild (and indeed unessential) parameter fitting, our model looks like a reasonable choice for Occam's razor to point toward. Yet, we cannot completely rule out that we were maliciously mislead by all those regularities and coincidences, that in reality the Schwabe cycle results from a peculiar self-synchronization mechanism (Hoyng, 1996), and that less "astrological" explanations will also be found for the long-term rhythms of our home star. particular peak at 193 years in the wavelet spectrogram (d). It seems that the transition to disorder/chaos is different in that case, and that the existence of a second frequency is decisive in this respect.