© C. Thibeault et al., Published by EDP Sciences 2023

1 Introduction

Solar and stellar dynamos are fundamentally nonlinear magnetofluid systems: the magnetic field is amplified and its energy grows through the work done by flows against the Lorentz force. As a consequence of this energy transfer, the dynamo-generated magnetic field unavoidably quenches the inductive flows powering it. That the saturation of solar and stellar magnetic cycles is achieved through this nonlinear interaction is generally agreed upon, but the details of the process are not. This is primarily because of the multiscale, turbulent nature of convective flows pervading the outer envelopes of the sun and solar-type stars, which power large-scale flows, also contributing to magnetic induction. Analysis of global 3D MHD simulation of solar/stellar convection reveals that the overall energy budget of the magnetic field and flows is indeed quite complex, and involve many intertwined channels operating concurrently; see, e.g., the extensive discussion in Brun et al. (2022).

In geometrically and dynamically simplified dynamo models, amplitude-limiting nonlinearities are typically introduced in a more or less ad hoc manner (see Sect. 4.2 in Charbonneau, 2020). Depending on the choice of nonlinearity made, a wide variety of behaviors can be produced, ranging from simple saturation at a constant cycle amplitude, period-doubling cascades, all the way to chaotic modulation and quasi-periodic modulations on timescales much longer than the primary cycle. For a good selection of representative examples, see, Küker et al. (1999), Brooke et al. (2002), Charbonneau et al. (2005); and in particular Karak (2023), and references therein.

This complexity, however, vanishes in situations where the dynamo is very weakly supercritical, in the sense that induction barely offsets dissipation. The ratio of these two processes defines the dynamo number, with dynamo action setting in when the said dynamo number exceeds some critical value that depends on the details of the model. This first bifurcation from the trivial fixed point B = 0 to a constant amplitude limit cycle beyond the critical dynamo number is believed to occur as a simple Hopf bifurcation, whatever the exact form of the nonlinearity, a situation shared by many classes of physical instabilities. In other words, the system behavior near this first bifurcation is generic (see, e.g., Tobias et al., 1995; Arlt & Weiss, 2014, and references therein).

Near a bifurcation point, a dynamo system becomes extremely sensitive to any external perturbation, and this is especially true of the first bifurcation from the trivial solution. A variety of studies have amply demonstrated that stochastic perturbations introduced in a critical or very weakly supercritical dynamo can produce many solar-like features, including marked cycle amplitude fluctuations, (anti)correlations between cycle amplitude and duration, mixed-parity modes, cross-equatorial activity, and on-off intermittency (see, e.g., Moss et al., 1992; Hoyng et al., 1994; Ossendrijver et al., 1996; Usoskin et al., 2009; Cameron & Schüssler, 2017; Karak, 2023, and references therein). Indeed, on the basis of an extensive set of simulations carried out using two distinct simplified dynamo models, Cameron & Schüssler (2017) have argued that weak supercriticality and stochastic forcing are all that is needed to reproduce the observed temporal spectrum of solar activity, independently of the nonlinearity at play. In particular, their model produces a flat spectrum of significant amplitude at long timescales (see Figs. 3 and 9 in Cameron & Schüssler, 2017), due to the slow response of the model to imposed random fluctuations. While their proposal does require some fine-tuning of model parameters, it is demonstrably generic with respect to the chosen nonlinearity.

However attractive the Cameron & Schüssler (2017) proposal might be in its simplicity, there is certainly no a priori reason to believe that the solar dynamo is very weakly supercritical. Moreover, solar activity reconstructions based on cosmogenic radioisotopes reveal patterns that are hard to reconcile with a structureless spectrum at periods longer than the primary 11-year cycle (Arlt & Weiss, 2014). These include intermittently recurring quasi-periodicities such as the so-called Suess-deVries cycle (timescale ∼ 210 years), the Hallstatt cycle (∼2400 years), as well as periods of strongly suppressed or enhanced activity known respectively as Grand Minima and Grand Maxima, the latter clustering around peaks of the Hallstatt cycle (see Usoskin, 2023, and references therein). Such modulation and supermodulation patterns are more typical of dynamos operating in the strongly supercritical regime (e.g., Weiss & Tobias, 2016; Beer et al., 2018).

In this paper, we focus on variability materializing in moderately to strongly supercritical dynamo models of the Babcock-Leighton variety, and in particular on the interactions of nonlinearity with weak stochastic forcing as a means of generating intermittency and Grand Maxima. Towards this end, and in the spirit of the Cameron and Schüssler (2017) simulations, in what follows we make use of two distinct models of Babcock-Leighton dynamos: the dynamical system model of Wilmot-Smith et al. (2006), and the mean-field-like two-dimensional model of Charbonneau et al. (2005). The former can be reduced to a single nonlinear second-order ordinary delay-differential equation, and thus allows a very efficient exploration of parameter space; while the latter, described by two nonlinearly coupled partial differential equations in two spatial dimensions, is computationally more demanding but approximates the sun somewhat more closely.

Babcock-Leighton dynamos, in which the regeneration of the large-scale poloidal magnetic component takes place through the decay of bipolar magnetic regions, represent an instance of flux transport dynamos in which the two magnetic source regions are spatially segregated. As discussed in more details in Section 2, this introduces an unavoidable time delay in the dynamo loop, which in itself can lead to period-doubling bifurcations and chaos in conjunction with even the simplest amplitude-limiting algebraic nonlinearities. The complex variability patterns resulting from this interaction of time delays and strong nonlinearity remain largely unexplored in such dynamo models (but do see Charbonneau et al., 2005; Albert et al., 2021; Tripathi et al., 2021).

The remainder of this paper is organized as follows. In Section 2, we briefly review the workings of the Babcock-Leighton mechanisms and briefly describe the general features and mathematical formulation of two distinct dynamo models built upon it and used in the paper. In Section 3, expanding on the results presented by Albert et al. (2021) and Tripathi et al. (2021), we summarize the properties of a vast set of simulations using the zero-dimensional dynamical system model originally developed by Wilmot-Smith et al. (2006), concentrating on the moderately to strongly supercritical regime. In Section 4, we then seek analogs of the behavior so uncovered in the spatially-extended Babcock-Leighton solar cycle model of Charbonneau et al. (2005). Based on these results, we propose in Section 5 a physical scenario for the generation of the Grand Maxima, based on intermittency. We conclude in Section 6 by framing our proposal in the context of other model-based scenarios put forth to explain solar cycle variability on timescales much longer than the period of the primary magnetic activity cycle.

2 Babcock-Leighton dynamos

2.1 The Babcock-Leighton mechanism

The Babcock-Leighton mechanism for the regeneration of the solar surface dipole moment (Babcock, 1961; Leighton, 1964; Sect. 5 in Charbonneau, 2020) relies on the decay of bipolar magnetic regions (BMR) emerging with a systematic pattern of inclination (or tilt) of the line segment joining each polarity centroid with respect the E–W direction. This tilt pattern is known as Joy’s Law (Hale et al., 1919, see also McClintock & Norton 2013, and references therein). The leading polarity being on average closer to the solar equator, it experiences greater cross-equatorial diffusive cancellation, so that an excess of trailing polarity slowly accumulates in polar regions through transport by surface flows, leading to the buildup of a global dipole moment. The dynamo loop is closed by rotational shearing of this global dipole by internal differential rotation generating the toroidal magnetic component which will give rise to BMRs emerging in the subsequent magnetic cycle. This dynamo scenario finds empirical support in the demonstrated good precursor potential of the observed surface dipole near sunspot cycle minimum (e.g., Schatten et al., 1978; Svalgaard et al., 2005; but see also Charbonneau & Barlet, 2011).

An important feature of such dynamo models is that the poloidal and toroidal source regions are spatially segregated so that the timescale associated with the transport process linking them ends up setting the cycle period. This introduces an unavoidable time delay in the dynamo loop, which in itself can produce complex modulation patterns (see Yoshimura, 1978; Durney, 2000; Charbonneau, 2001; Charbonneau et al., 2005; Wilmot-Smith et al., 2006; Jouve et al., 2010).

2.2 The nonlinearity

A key aspect of the Babcock-Leighton mechanism is its dependence on the magnetic field strength in the solar interior, where the magnetic flux ropes giving rise to emerging BMRs undergo buoyant destabilization and subsequent rise across the convection zone (see Fan, 2021, and references therein). Modeling under the thin flux tube approximation has revealed that a Joys Law-like tilt pattern only materializes in a finite range of internal field strengths. Magnetic flux tubes too strongly magnetized rise too rapidly through the convection zone and emerge without the E-W tilt essential to the Babcock-Leighton mechanism. Too weakly magnetized flux tubes fail to destabilize, or are destroyed by turbulent fluid motions within the convection zone. Obsevationally, evidence for tilt quenching is present but weak, and is sensitive to the dataset and data reduction technique used (cf. Dasi-Espuig et al., 2010; McClintock & Norton, 2013; Jha et al., 2020).

In what follows, we capture these behaviors through a generalization of the algebraic nonlinearity introduced by Wilmot-Smith et al. (2006) and also used by Albert et al. (2021). The Babcock-Leighton source term S(B) for the poloidal magnetic component is written as:

$$S\left(r,\theta ;B\right)={\alpha }_{0}f\left(B\right)B g\left(r,\theta \right),$$(1)

where B is the internal toroidal component and the function g sets the spatial dependency of the source term in the spatially extended dynamo models introduced in Section 2.3. The nonlinear dependency of S on B is parameterized via the function f(B) as:

$$f\left(B\right)=\frac{1}{4}\left[1+\mathrm{erf}\left(\frac{B^2-B_1^2}{w_1^2}\right)\right]\left[1-\mathrm{erf}\left(\frac{B^2-B_2^2}{w_2^2}\right)\right],$$(2)

where erf(x) is the error function. The parameters w₁ and w₂ control the sharpness of the transition at the lower (B₁) and upper (B₂) thresholds. Figure 1 shows three examples for specific parameter choices. Note that if $w_1^2<<B_1^2$ and B₁ > 0, f(B) → 0 as B → 0, the dynamo will no longer be self-excited, i.e., it will fail to amplify an arbitrarily weak seed magnetic field.

Figure 1 The nonlinearity defined by equation (2) and used in the two dynamo models described in Section 2.2. Setting w1 = w2 ≃ 1 yields sharp thresholds, while larger values lead to smoother transitions.

2.3 The axisymmetric mean-field dynamo equations

Dynamo models of the solar cycle can be constructed within the general framework of mean-field electrodynamics, by replacing the α-effect a suitable source term modeling the operation of the Babcock-Leighton mechanism (see, e.g., Charbonneau, 2020, Sect. 5.4). Working in spherical polar coordinates (r, θ, ϕ) the large-scale axisymmetric magnetic field B(r, θ, t) is expressed as B = ∇ × [A(r, θ, t)ê_ϕ] + B(r, θ, t)ê_ϕ, which allows to separate the induction equation into:

$$\frac{\partial A_{\varphi }}{\partial t}+\frac{1}{\varpi }\left({\vec{\mathbf{u}}}_{\mathbf{p}}\cdot \nabla \right)\left(\varpi A_{\varphi }\right)=\eta \left({\nabla }^2-\frac{1}{{\varpi }^2}\right)A_{\varphi }+S\left(B\right),$$(3)

$$\frac{\partial B_{\varphi }}{\partial t}+\varpi \nabla \cdot \left(\frac{{\vec{\mathbf{u}}}_{\mathbf{p}}{B}_{\varphi }}{\varpi }\right)=\eta \left({\nabla }^2-\frac{1}{{\varpi }^2}\right)B_{\varphi }+\frac{1}{r} \frac{\partial \eta }{\partial r}\frac{\partial \left(r{B}_{\varphi }\right)}{\partial r}$$(4)

$$+\varpi \left(\nabla \times A\widehat{\varphi }\right)\cdot \left(\nabla \mathrm{\Omega }\right),$$(5)

where Ω(r, θ) is the (axisymmetric) solar internal angular velocity, u_p(r, θ) is the (axisymmetric) meridional flow, η(r) is the net (turbulent) magnetic diffusivity, and S(B) is the Babcock-Leighton source term S(B) defined through equations (1) and (2).

3 A dynamical system model

As a guide to the design and interpretation of the numerical experiments to be reported upon in Sections 4 and 5, we first use a dynamical system model allowing a (relatively) straightforward characterization of the interaction between the Babcock-Leighton nonlinearity and time delays, as well as a computationally efficient exploration of parameter space.

Wilmot-Smith et al.’s (2006) dynamo model is a heuristic spatial averaging of the mean-field dynamo equations (3) and (4), which become ordinary coupled differential equations (ODEs) for the surface poloidal and internal toroidal magnetic amplitudes A(t) and B(t):

$$\frac{\mathrm{d}A}{\mathrm{d}t}={\alpha }_{0}f\left(B\left(t-T_1\right)\right)B\left(t-T_1\right)-\frac{A}{\tau },$$(6)

$$\frac{\mathrm{d}B}{\mathrm{d}t}=\frac{{\mathrm{\Omega }}_0}{L}A\left(t-T_0\right)-\frac{B}{\tau },$$(7)

where L is a length scale over which the angular velocity varies by Ω₀, thus setting the magnitude of the rotational shear. As detailed in Wilmot-Smith et al. (2006), the parameters T₀ and T₁ measure the time delays in the poloidal → toroidal and toroidal → poloidal inductive mechanisms, respectively, producing the dynamo loop. Choosing the magnetic diffusion time τ as a time unit, equations (6) and (7) can be combined in a single second-order nonlinear ODE:

$$\frac{{\mathrm{d}}^{2}B}{\mathrm{d}{t}^2}+2\frac{\mathrm{d}B}{\mathrm{d}t}+B=\mathrm{Nf}\left(B\left(t-q\right)\right)B\left(t-q\right),$$(8)

with the ratio q between the total time delay and diffusion time,

$$q=\frac{T_0+T_1}{\tau },$$(9)

determining whether the dynamo is “advection-dominated” (q < 1) or “diffusion-dominated” (q > 1). The dynamo number N for this model is given by:

$$N=\frac{{\alpha }_0{\mathrm{\Omega }}_0{\tau }^2}{L} .$$(10)

This (relatively) simple dynamical system model thus neatly captures, in simplified form, the non-locality and nonlinearity of the Babcock-Leighton mechanism, as well as the time-delay characterizing the coupling of the two spatially-segregated source regions.

Although equation (8) looks straightforward, the embodied time delay q makes it non-trivial to solve numerically. After testing a few packages for time-delay ODEs, we selected the Python package Pydelay¹ which, we found, offers a robust balance between accuracy and run-time performance, at least for our ODE (8). Empirical bifurcations diagrams are constructed by integrating equation (8) in time from some non-zero initial condition, removing the transient phase, collecting all B² maxima and plotting them against dynamo number N, repeating this process for a sequence of solutions with increasing N.

Not surprisingly, the shape of the bifurcations diagrams depends on the numerical values chosen for the various model parameters, but some important features turn out to be robust with regard to parameter settings. Figure 2 shows four representative examples, with associated parameter values listed in Table 1.

Figure 2 Bifurcation diagrams for the four selected parameter sets, as listed in Table 1, constructed from a series of dynamo solutions with increasing dynamo number N. The insets shown in (B) and (C) are closeups of the multiperiodic/chaotic window shown in the same panel. In this model, without a lower cutoff on the nonlinearity, the limit cycle emerging from the first bifurcation (N/Ncrit = 1) never undergoes further bifurcation or transition to chaos.

When the time delay parameter q is set at values significantly smaller than unity (Fig. 2A; advection-dominated regime), a single dynamo branch emerges from the first bifurcation and persists all the way to the strongly supercritical regime. The co-existing trivial solution B = 0 retains a small basin of attraction at all N, a consequence of the lower operating threshold characterizing the nonlinearity (2) when w₁ ≲ B₁. A conspicuous kink at N/N_crit ≃ 2 appear in the dynamo branch here at q ∼ 0.4. This is the “strong branch” already identified by Albert et al. (2021) in their investigation for distinct defining parameters for their nonlinearity (viz. Table 1, last line). As q approaches unity (the diffusion-dominated regime), period-doubling and transition to chaos sets in over a finite window in N, extending to higher and higher N as q is further increased. The occurrence of this multiperiodic/chaotic window is robust with respect to the adopted values for the nonlinearity parameters, although its extent in N as well as the bifurcation structure itself can both vary significantly (cf. Figs. 2B and 2C). As exemplified by Figure 2B, multiple stable dynamo branches typically survive beyond the multiperiodic/chaotic window, each with a finite-sized basin of attraction, and sometimes undergoing further period doubling and transition to chaos as N is further increased into the strongly supercritical regime N/N_crit ≳ 5. Most of these features remain as q exceeds unity, with the chaotic window increasing in complexity, and becoming delineated by increasingly broad multiperiodic windows (Fig. 2D). The “strong mode”, high amplitude branch stands out clearly in all these diagrams, bifurcating from the main attractor typically at N/N_crit ≃ 2–5 and commanding the dominant basin of attraction at larger N. We have pushed the dynamo number far beyond the right extent of Figure 2, and have not found further bifurcation of this singly-periodic branch.

The co-existence of multiple stable dynamo branches at a given dynamo number N has other interesting consequences if low-amplitude random noise is introduced in the system, as exemplified already by the simulations of Albert et al. (2021). Figure 3A shows two phase space trajectories for the two stable dynamo branches at N/N_crit = 2.99 in Figure 2C (vertical red line). The axes of this 3D phase space are chosen as the squared magnetic field amplitude B²(t), the same amplitude temporally lagged by q, and the time derivative of B²(t). Distinct initial conditions B_i lead the solution to converge on either the high or low branch. The trajectories are distinct and do not intersect, but do approach one another during the minimum phase of their respective magnetic cycles. This opens the possibility of transiting from one trajectory to the other through perturbations introduced by low-amplitude stochastic noise.

Figure 3 (A) Phase space trajectories for two 0D dynamo solutions with N/Ncrit = 2.99 (see red line segment on Fig. 2C). The two solutions differ only in their initial condition, each located in a distinct basin of attraction. (B) A representative time series of the same dynamo solution, now with additive stochastic noise introduced on the RHS of equation (8). The orange lines are added to delineate the low and high amplitude epochs. (C) Probability density function for the duration of sustained high amplitude epochs, as extracted from panel B.

Figure 3B shows an example, in the form of a very long time series of magnetic energy (simply B² in this model) in a N/N_crit = 2.99 solutions akin to Figure 3A, but now with zero-mean additive noise uniformly distributed in δα/α₀ ∈ [−0.5, 0.5] introduced on the RHS of equation (8), with a very short coherence time of 0.18% of a full magnetic cycles. Individual cycles cannot be distinguished on this plot, but the partition between low and high amplitude epochs is obvious, delineated here in orange. Constructing histogram for the duration of high amplitude phases (Fig. 3C) yields an exponential distribution, as expected if the transitions are triggered by a stationary random process. Sensitivity to noise is greatest near the bifurcation point of the strong mode branch, but the transitions can also be stochastically triggered at higher values of the dynamo number N.

This brief survey of some of the nonlinear behaviors observed in the Wilmot-Smith et al. (2006) dynamical system model highlights the rich set of amplitude modulation patterns that can be produced by the interaction between a strong nonlinearity and a time delay in the dynamo loop. It must be pointed out that these complex patterns materialize in dynamo solutions that are not very far into the supercritical regime. Depending on parameter values, period doubling and transition to chaos can already set in for mildly supercritical models (viz. Fig. 2); chaotic modulation is already ubiquitous at N/N_crit = 2–3 (viz. Figs. 2B, 2C, 2D), unless q is much smaller than unity; stochastic excitation is possible at N/N_crit ≃ 3 (Fig. 3).

One may rightfully wonder how much of this carries over to spatially extended dynamo models including a more realistic solar-like internal differential rotation profile, as well as spatiotemporal smoothing by magnetic diffusion. The limited set of two-dimensional mean-field-like Babcock-Leighton dynamo solutions presented in Charbonneau et al. (2005) suggest that time delay is indeed a robust mechanism for the production of deterministic cycle amplitude modulations in this context. We pursue along these lines by searching for counterparts of behaviors surveyed in this section in the mean-field-like dynamo model described in Section 2 herein.

4 Supercritical solutions in 2D mean-field-like dynamo model

We now turn to numerical solutions of the two-dimensional axisymmetric kinematic dynamo equations (3) and (4). Following Charbonneau et al. (2005), the dynamo equations are cast in dimensionless form using the magnetic diffusion time τ_D = R²/η_T as a time unit, where η_T is the (turbulent) magnetic diffusivity within the convection zone, and the solar radius R as a length unit. Three dimensionless grouping then emerge:

$$C_{\alpha }=\frac{{\alpha }_{0}R}{{\eta }_{\mathrm{T}}}, C_{\mathrm{\Omega }}=\frac{{\mathrm{\Omega }}_0{R}^2}{{\eta }_{\mathrm{T}}}, \mathrm{R}m=\frac{u_{0}R}{{\eta }_{\mathrm{T}}},$$(11)

where α₀ sets the magnitude of the Babcock-Leighton source term (viz. Eq. (1)), and u₀ is the surface mid-latitude meridional flow speed, fixed at a solar-like value u₀ = 10³ cm s⁻¹ in all simulations. The dynamo number (D) for this model is given by the product of C_α and C_Ω:

$$D=\frac{{\alpha }_0{\mathrm{\Omega }}_0{R}^3}{{\eta }_{\mathrm{T}}^2},$$(12)

effectively identical to equation (10) for the 0D model. The parameterized meridional circulation and differential rotation profiles are solar-like the latter including a tachocline straddling the interface between the convective envelope and an underlying fluid layer where the magnetic diffusivity drops by a factor Δη from its envelope value η_T. See Dikpati & Charbonneau (1999) and Charbonneau et al. (2005) for more details.

This spatially extended 2D dynamo model is far more demanding computationally than the dynamical system model of Wilmot-Smith et al. (2006), so we opted to restrict our parameter space to the Babcock-Leighton nonlinearity, fixing B₁ = 100, w₁ = 10 and varying B₂, w₂. Inspired by simulations of the stability and buoyant rise of thin flux tubes in the deep convection zone (see Fan, 2021, and references therein), we set the ratios B₂/B₁ ≃ 10 and w₂/w₂ ≫ 1. We compute solutions for three values of envelope magnetic diffusivity, log(η_T) = 11.0, 11.6, and 12.0 (η_T in cm² s⁻¹), spanning the transition from advection-dominated to diffusion-dominated regimes. Numerical solutions are restricted to the Northern hemisphere, and dipolar parity is enforced via the boundary condition on the equatorial plane.

For each set of model parameters listed in Table 2, we compute sequences of dynamo solutions for increasing values of the dynamo number D, and for each sequence we produce a time series of volume-integrated magnetic energy. For each sequence, a bifurcation diagram is then constructed in the same manner as described in Section 3. The nine resulting bifurcation diagrams are shown in Figure 4, organized so that magnetic diffusivity is the same in each column but increases from left to right, and parameters of the Babcock-Leighton nonlinearity are the same in each row.

Figure 4 Bifurcation diagrams for each parameter set are shown in Table 2. Dynamo numbers (D) are normalized to the critical value (Dcrit) for each sequence, as listed in the rightmost column of Table 2, The red and blue lines highlight the simulations chosen as case studies for Figure 5. The panels are organized such that magnetic diffusivity is the same in each column, and increases from left to right.

While the bifurcation diagrams show a wide range of morphologies and branching structures, some common (and robust) features are present in all cases and resemble features characterizing the dynamical system model of Section 3. Transition to growing solutions at the critical dynamo number takes place as a Hopf-type bifurcation to a fixed amplitude limit cycle, with the saturated amplitude growing steadily as the dynamo number increases along each sequence. Further bifurcation takes place at a few times criticality, to either a new singly periodic branch (viz. the 116 sequences), or chaotic regime, either directly (110 sequences) or through a period-doubling cascade (120 sequences). In the more strongly supercritical regime (D/D_crit ≳ 5, say), most sequences revert to a high-amplitude single branch, undergoing further period-doubling and transition to chaos in the 116 sequence, but remaining stable up to the highest dynamo numbers simulated in the case of the 120 sequence. This is the 2D model’s counterpart to the “strong mode” identified in the 0D model of Section 3.

Figure 5 shows magnetic energy time series, for two selected members in three sequences of Figure 4, as indicated by the the vertical colored line segments. Figure 5A shows two chaotic dynamo solutions from the B110 sequence, one moderately supercritical (D/D_crit = 1.4, in blue), the other extracted from a distinct, more strongly supercritical chaotic window (red). The former is a typical chaotic oscillation, while the second shows a quasi-periodic modulation of periods of ≈0.18 τ_D superimposing itself on a chaotic modulation.

Figure 5 Time series of magnetic energy (red and blue) and mid-latitude toroidal magnetic component at r/R = 0.7 (dotted lines), for samples of numerical solutions chosen from the set of solutions shown in Figure 4. All use the same parameter values for the Babcock-Leighton nonlinearity, and differ only in magnetic diffusivity and dynamo numbers, as indicated. The inset label of each figure represents the set of parameters used for the simulation (see Table 2), and the colors of the line map to the vertical line highlights from Figure 4.

Figure 5B shows two solutions from the B116 sequence, one on each side of but very close to the prominent bifurcation at D/D_crit = 5.4 (see central panel on Fig. 4). Below this bifurcation, the cycle shape is almost sinusoidal in magnetic energy; while above it the rising phase is much steeper than the descending phase, the peak amplitude has nearly doubled, and, more surprisingly, the cycle period also seems to have doubled. In fact, it has not, as can be seen from the toroidal field time series (extracted at r/R = 0.7 and θ = π/4) plotted as a dotted line. Upon crossing the bifurcation, the oscillation simply develops a strong bias in magnetic polarity, to the extent that the half-cycle associated with the weaker polarity does not manage to produce a distinct peak in the volumetric energy time series (more on this further below).

Finally, in the strongly diffusive solution pair of Figure 5C, the first solution shows a regular cycle periodically undergoing brief episodes of strongly reduced amplitudes. The second, with a dynamo number a mere ≃5 percent higher, shows a constant amplitude. Both solutions are also characterized by strong magnetic polarity asymmetries, as for the bottom solution of Figure 5B.

In Figure 6, we further select three solutions from Figure 5 and display them now as a time-latitude diagram of the toroidal magnetic component, extracted at the core-envelope interface r/R = 0.7 (top diagrams in Figs. 6A, 6B, and 6C). The B110 advection-dominated dynamo solution (Figs. 6A) powers itself primarily from the low-latitude radial shear in the tachocline, the more diffusive B116 solution (Fig. 6B) feeds on the latitudinal shear in the bottom convective envelope, while the even more diffusive B120 solution (Fig. 6C) is powered primarily by the radial shear in the high-latitude region of the tachocline. In this solution, the latitudinal shear in the bottom half of the convective envelope drives a secondary, longer period dynamo oscillatory mode, which interferes with the primary high-latitude mode. The resulting beating phenomenon causes the periodic drop-in magnetic energy so conspicuous on Figure 5C (top time series), as well as on the time-latitude-diagram in Figure 6C. In all cases, the advective effect of the equatorward meridional flow drives equatorward propagation at mid to low latitudes, an effect the more pronounced the lower the magnetic diffusivity. The sustained asymmetry in magnetic polarity mentioned previously is clearly apparent in the B116 dynamo solution of Figure 6B.

Figure 6 Northern hemisphere time-latitude “butterfly” diagrams for three sample solutions taken from Figure 5, all using the same Babcock-Leighton nonlinearity but differing in their envelope magnetic diffusivity and dynamo number (viz. Table 2). For each simulation shown, the top panel shows the toroidal magnetic field at R/R⊙ = 0.7, and the bottom panel shows the source term (|f(Bϕ)Bϕ|) at the same depth. Note that the color scale for the toroidal magnetic field differs in each panel.

We also plot In Figure 6, below each time-latitude diagram of the toroidal magnetic component, a time-latitude diagram of the magnetic part of the Babcock-Leighton source term, (|f(B_ϕ) × B_ϕ|), as given by equation (1), also extracted at r/R = 0.7. These illustrate a fundamental dynamical difference distinguishing the “strong mode” from the primary dynamo mode emerging from the first critical Hopf bifurcation. For the solution of Figure 6A, residing on the primary branch, the source term is strongest at epochs of peak toroidal field, as the latter remains too weak to push beyond the upper threshold B₂ of the Babcock-Leighton nonlinearity (viz. Eq. (2) and Fig. 1). In panel B, however, the source term is strongly suppressed when the toroidal fields approaches and exceeds the upper threshold B₂. This leads to a strong quenching of the poloidal source term, and thus a weaker subsequent cycle. This weaker cycle then yields a stronger source term, producing a stronger subsequent cycle, and so on, leading to a persistent Gnevyshev-Ohl-like 2-cycle modulation and marked asymmetry in magnetic polarity. The fact that the positive magnetic polarity is here dominant is a consequence of the initial condition, with dominance of the negative polarity being equally likely.

Figure 7 depicts the “sampling” of the Babcock-Leighton nonlinearity f(B) along fixed-latitude temporal cuts along the two time-latitude diagrams of Figures 6A and 6B, at the latitude of peak toroidal field strength in each dynamo solution. The light gray histograms display the sampling of toroidal field values along the cuts, and the color shades the corresponding values of f(B), the intensity of shading depicting the proportion of the time spent at a given point along the non-linearity curve. For the B110 chaotic solution of Figure 7A, the distributions are nearly symmetrical about B_ϕ = 0, and most cycle maxima of either polarity, corresponding to the two peaks in the histogram, occur while f(B) is of order unity. In contrast, for the strongly supercritical “strong mode” dynamo solution B116 (D/D_crit = 5.5), a marked and systematic asymmetry readily distinguishes the positive and negative magnetic polarity half-cycles, with cycle maxima occurring while f(B) ≃ 0. Note in particular how, on Figure 7B, the solutions extend far in the fully-quenched extent of the Babcock-Leighton nonlinearity, in contrast to the B110 solution of Figure 7A. These two sample solutions jointly highlight the two distinct dynamical modes accessible in these supercritical dynamo solutions: a primary mode mostly residing in the [B₁, B₂] interval of the Babcock-Leighton nonlinearity, and a second “strong” mode spending a significant fraction of its magnetic cycle in the strongly-quenched range (B_ϕ > B₂) of the nonlinearity.

Figure 7 Sampling of the Babcock-Leighton nonlinearity f(B) for the two time-latitude diagrams of Figures 6A and 6B, along temporal cuts at latitudes 25° and 45°, respectively. The histograms display the sampling of toroidal field values along the cuts, and the color shading indicates the fraction of time spent at a given point along the non-linearity curve (dotted line). Note the marked magnetic polarity asymmetry in panel B.

5 Intermittency and Grand Maxima

Having demonstrated the co-existence of dynamo modes characterized by distinct nonlinear dynamics, we now consider the possibility of transit from one to the other via the agency of stochastic forcing, as was shown possible in the simpler 0D dynamical system model of Section 3. Towards this end, we adopt the procedure described in Charbonneau & Barlet (2011), which consists of adding a fluctuation δα to the poloidal source term amplitude parameter (α₀ in Eq. (1)). This fluctuation is extracted from a zero-mean uniform distribution and its value is held fixed during a set coherence time τ, at the end of which it is reset to a new value.

Figure 8 shows the result of such an experiment. Here additive stochastic fluctuation at the level δα/α₀ = ±0.2 is introduced in the Babcock-Leighton source term, with coherence time τ/τ_D = 1.22 × 10⁻⁴, corresponding to 0.26% of the magnetic cycle period in the unperturbed dynamo solution.

Figure 8 (A) Time series of magnetic energy (blue) and the toroidal magnetic field at mid-latitude and at R = 0.7 R⊙ (red). (B) Time series of smoothed full-cycle amplitude Ak (see text for definition). The gray shading corresponds to the time span plotted in (A). The red and green dots represent cycles higher and lower (respectively) than the Ak-threshold chosen to delineate epochs of strong and normal mode of dynamo operation. The choice of threshold value is based on a break in the slope of the |Ak| histogram, plotted in (C). Panels D and E show histograms of strong mode durations and inter-event waiting times, equivalent to the duration distribution for normal dynamo operation.

The top panel shows a portion of the time series of magnetic energy (blue), and toroidal magnetic field (red) extracted at r/R = 0.7 and 45° latitude. This is a small segment from a much longer simulation, spanning 4111 full magnetic cycles and 213 diffusion times. In light of our analysis of the preceding section, we define a measure of magnetic polarity asymmetry, on the basis of which we seek to distinguish epochs spent in normal and strong dynamo modes. Working off the toroidal field time series plotted on red in Figure 8A, we integrate temporally over each successive full magnetic cycle to produce a temporal sequence, to which we apply a 1-4-6-4-1 averaging filter. The resulting smoothed time series displays the evolution of full-cycle polarity asymmetry, which we denote A_k, where k runs over 4111 full magnetic cycles. The results of this procedure are plotted in Figure 8B. Cycles where A_k fluctuates about zero, plotted in red, indicate polarity-symmetric cycles, corresponding to the normal dynamo mode, while A_k departing significantly from zero for extended time periods, plotted in green, is indicative of strong mode operation.

Figure 8C shows a histogram of |A_k| values for all cycles in the simulation. A conspicuous break of slope at |A_k| ≃ 65,000 suggests the use of this numerical value to delineate strong mode operation from the primary dynamo mode. These are indicated by orange horizontal lines at ±65,000 in Figure 8B. It is then a simple matter to construct histograms of strong mode epochs duration and inter-event waiting times (the latter equivalent to duration of normal mode epochs). Both histograms, plotted in Figures 8D and 8E, are exponential in form, consistent with transitions between the two dynamo modes being driven by a stationary Poisson process.

The fraction of time spent in strong mode versus normal mode depends on parameter values for the nonlinearity and dynamo model, with, typically, the strong mode dominating in the strongly supercritical regime D/D_crit ≫ 1. However, based on the relatively limited numerical explorations carried out so far, in the moderately supercritical regime the ability of stochastic noise to trigger transitions between the two modes appears to be a robust feature.

6 Discussion and conclusion

In this paper, we have focused on the behavior of dynamo models of the solar cycle based on the Babcock-Leighton mechanism for the regeneration of the Sun’s poloidal magnetic component. Towards this end, we have used two distinct dynamo models: the dynamical system model of Wilmot-Smith et al. (2006), and the spatially-extended two-dimensional mean-field-like model of Charbonneau et al. (2005). Both models are kinematic, in the sense that they implement temporally steady large-scale differential rotation and meridional circulation, and also use the same non-local amplitude-limiting nonlinearity on the poloidal source term.

As their respective dynamo numbers are increased, both dynamo models exhibit rich bifurcation structures, moreover showing remarkable structural similarities, in turn suggesting robustness with respect to model formulation and details. The key underlying drivers of this complex dynamical behavior are an amplitude-limiting nonlinearity including both upper and lower thresholds, as well as the time delay inherent to this class of flux transport dynamos, in which inductive source regions are spatially segregated.

Of particular importance is the co-existence, over wide ranges of dynamo numbers, of multiple dynamo branches each with a finite basin of attraction. The trivial solution B = 0 retains a finite basin of attraction even in the supercritical regime, because of the lower threshold of magnetic field strength characterizing the Babcock-Leighton mechanism. The primary dynamo branch, emerging as the critical dynamo number is exceeded, typically undergoes a period-doubling transition to chaos until a new bifurcation gives rise to the “strong” dynamo branch, which rapidly reaches much higher amplitude than the original branch from which it bifurcates. In the strongly supercritical regime, this strong mode branch commands the largest basin of attraction, and can undergo further bifurcations and transition to chaos. Nonetheless, over wide portions of the model’s parameter space, the original branch remains stable beyond the bifurcation of the strong branch, and retains a finite basin of attraction, as does the trivial solution B = 0.

In the strongly supercritical regime of our two-dimensional models, we sometimes observe quasiperiodic cycle amplitude modulation on long timescales, i.e., many times the primary cycle period; this is typically observed in non-kinematic dynamos where magnetic backreaction quenches large-scale inductive flows (see, e.g., Tobias, 1997; Bushby, 2006; Simard & Charbonneau, 2020; and references therein), the long timescale modulation resulting from the time required to re-establish the large-scale flows after the magnetic amplitude has fallen in response to the quenching of the large-scale flows. In contrast, in our purely kinematic model long timescale modulation in the strongly supercritical regime typically results from the interaction of distinct dynamo modes subject to different time delay dynamics because of their localization within the modeled convection zone. The observed long timescale modulations are then akin to a beating phenomenon.

Introducing stochasticity into the two dynamo models can lead to transitions between the distinct dynamo branches co-existing at a given dynamo number. Working also with the Wilmot-Smith et al. (2006) dynamical system model but in the mildly subcritical regime, Tripathi et al. (2021) show that the combination of stochastic fluctuation in dynamo number and stochastic additive noise can generate intermittency, leading to Grand Minima-like episodes of strongly suppressed cyclic activity when the solution enters the attraction basin of the trivial solution and collapses to B = 0. Cyclic activity resumes when a favorable realization of the stochastic additive noise pushes the solution back into the attraction basin of the primary dynamo branch. This scenario for Grand Minima is much along the lines of those suggested by Ossendrijver (2000) and Charbonneau et al. (2004), and is consistent with solar activity reconstructions based on cosmogenic radioisotopes, in which Grand Minima stand out as a dynamically distinct dynamo state (see, e.g., Usoskin et al., 2014).

Also working with the Wilmot-Smith et al. (2006) model, Albert et al. (2021) show how small periodic modulations of the dynamo number acting jointly with low amplitude additive white noise can trigger transitions between dynamo branches when operating close to the bifurcation point where the strong branch separates from the main attractor (viz. Fig. 4 herein). Albert et al. (2021) identify the strong mode with “normal” cyclic dynamo behavior, and the “weak mode” with Maunder-like periods of strongly suppressed cyclic activity.

Our proposal combines and extends the ideas put forth by Albert et al. (2021) and Tripathi et al. (2021). We associate Grand Maxima of solar activity (Usoskin, 2023) with the strong mode, normal cyclic behavior with the primary attractor emerging from the first bifurcation, and Grand Minima episodes with the trivial solution branch. Depending to some extent on model parameters, for moderately supercritical solutions (D/D_crit ≃ 3–6) our simulations indicate that these three branches can co-exist over an extended range of dynamo numbers (viz. Fig. 2B) and each retain a finite-sized basin of attraction. It would be particularly interesting to introduce non-kinematic effects in the 2D model, to ascertain to what extent such effect could lead to additional long-timescale modulations, and, acting in conjunction with stochastic forcing, lead to clustering of Grand Minima episodes in troughs of this long timescale modulations.

Acknowledgments

The work reported upon in this paper was funded by the Natural Sciences and Engineering Research Council of Canada (grant number RGPIN/05278-2018; PC), the Fonds de Recherche du Québec Nature et Technologies (graduate fellowship number 287277; CT), and a Hydro-Québec doctoral merit scholarship (CT). This research was also facilitated by support from the International Space Science Institute (ISSI Team 474). PC also acknowledges stimulating exchanges with A. Ferriz-Mas on the subject matter of this paper. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

	Figure 1 The nonlinearity defined by equation (2) and used in the two dynamo models described in Section 2.2. Setting w1 = w2 ≃ 1 yields sharp thresholds, while larger values lead to smoother transitions.
In the text

Figure 2 Bifurcation diagrams for the four selected parameter sets, as listed in Table 1, constructed from a series of dynamo solutions with increasing dynamo number N. The insets shown in (B) and (C) are closeups of the multiperiodic/chaotic window shown in the same panel. In this model, without a lower cutoff on the nonlinearity, the limit cycle emerging from the first bifurcation (N/Ncrit = 1) never undergoes further bifurcation or transition to chaos.

In the text

Figure 3 (A) Phase space trajectories for two 0D dynamo solutions with N/Ncrit = 2.99 (see red line segment on Fig. 2C). The two solutions differ only in their initial condition, each located in a distinct basin of attraction. (B) A representative time series of the same dynamo solution, now with additive stochastic noise introduced on the RHS of equation (8). The orange lines are added to delineate the low and high amplitude epochs. (C) Probability density function for the duration of sustained high amplitude epochs, as extracted from panel B.

In the text

Figure 4 Bifurcation diagrams for each parameter set are shown in Table 2. Dynamo numbers (D) are normalized to the critical value (Dcrit) for each sequence, as listed in the rightmost column of Table 2, The red and blue lines highlight the simulations chosen as case studies for Figure 5. The panels are organized such that magnetic diffusivity is the same in each column, and increases from left to right.

In the text

Figure 5 Time series of magnetic energy (red and blue) and mid-latitude toroidal magnetic component at r/R = 0.7 (dotted lines), for samples of numerical solutions chosen from the set of solutions shown in Figure 4. All use the same parameter values for the Babcock-Leighton nonlinearity, and differ only in magnetic diffusivity and dynamo numbers, as indicated. The inset label of each figure represents the set of parameters used for the simulation (see Table 2), and the colors of the line map to the vertical line highlights from Figure 4.

In the text

Figure 6 Northern hemisphere time-latitude “butterfly” diagrams for three sample solutions taken from Figure 5, all using the same Babcock-Leighton nonlinearity but differing in their envelope magnetic diffusivity and dynamo number (viz. Table 2). For each simulation shown, the top panel shows the toroidal magnetic field at R/R⊙ = 0.7, and the bottom panel shows the source term (|f(Bϕ)Bϕ|) at the same depth. Note that the color scale for the toroidal magnetic field differs in each panel.

In the text

Figure 7 Sampling of the Babcock-Leighton nonlinearity f(B) for the two time-latitude diagrams of Figures 6A and 6B, along temporal cuts at latitudes 25° and 45°, respectively. The histograms display the sampling of toroidal field values along the cuts, and the color shading indicates the fraction of time spent at a given point along the non-linearity curve (dotted line). Note the marked magnetic polarity asymmetry in panel B.

In the text

Figure 8 (A) Time series of magnetic energy (blue) and the toroidal magnetic field at mid-latitude and at R = 0.7 R⊙ (red). (B) Time series of smoothed full-cycle amplitude Ak (see text for definition). The gray shading corresponds to the time span plotted in (A). The red and green dots represent cycles higher and lower (respectively) than the Ak-threshold chosen to delineate epochs of strong and normal mode of dynamo operation. The choice of threshold value is based on a break in the slope of the |Ak| histogram, plotted in (C). Panels D and E show histograms of strong mode durations and inter-event waiting times, equivalent to the duration distribution for normal dynamo operation.

In the text