Issue 
J. Space Weather Space Clim.
Volume 13, 2023
Topical Issue  Space Climate: Longterm effects of solar variability on the Earth’s environment



Article Number  32  
Number of page(s)  13  
DOI  https://doi.org/10.1051/swsc/2023029  
Published online  22 December 2023 
Research Article
Nonlinearity, time delay, and Grand Maxima in supercritical BabcockLeighton dynamos
^{1}
Département de physique, Université de Montréal, C.P. 6128 Succ. Centreville, Montreal, QC H3C 3J7, Canada
^{2}
Department of Physics, McGill University, 3600 rue University, Montreal, QC H3A 2T8, Canada
^{*} Corresponding author: christianthibeaultnb@gmail.com, christian.thibeault@umontreal.ca
Received:
25
May
2023
Accepted:
15
November
2023
The physical origin of centennial and millennialscale variations in solar activity remains illunderstood. Although stochastic fluctuations of the solar dynamo are unavoidable in view of the turbulent nature of the solar convection zone, the quasiperiodic long timescale modulations revealed by the cosmogenic radioisotope records are suggestive of a deterministic process. In this paper, we investigate the nonlinear behavior of two solar cycle models based on the BabcockLeighton mechanism, with particular emphasis on deterministic amplitude modulation patterns materializing in the moderately to strongly supercritical dynamo regimes. Although formulated quite differently, both models show common long timescale modulation patterns arising from the interaction between the timedelay dynamics inherent to these flux transport dynamos, with the threshold nonlinearity characterizing the BabcockLeighton mechanism of poloidal field regeneration. In particular, we demonstrate the existence of multiple coexisting dynamo branches in the supercritical regime, each retaining a finitesized basin of attraction over a substantial range in dynamo number. The transition from one branch to another is shown to be possible via the introduction of lowamplitude stochastic noise with short coherence time. On this basis, we propose a novel physical scenario potentially accounting for the occurrence of both Grand Minima and Maxima of solar activity.
Key words: Solar activity / Solar dynamo / Space Climate
© C. Thibeault et al., Published by EDP Sciences 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 dynamogenerated 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 solartype stars, which power largescale 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, amplitudelimiting 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, perioddoubling cascades, all the way to chaotic modulation and quasiperiodic 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 solarlike features, including marked cycle amplitude fluctuations, (anti)correlations between cycle amplitude and duration, mixedparity modes, crossequatorial activity, and onoff 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 finetuning 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 11year cycle (Arlt & Weiss, 2014). These include intermittently recurring quasiperiodicities such as the socalled SuessdeVries 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 BabcockLeighton 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 BabcockLeighton dynamos: the dynamical system model of WilmotSmith et al. (2006), and the meanfieldlike twodimensional model of Charbonneau et al. (2005). The former can be reduced to a single nonlinear secondorder ordinary delaydifferential 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.
BabcockLeighton dynamos, in which the regeneration of the largescale 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 perioddoubling bifurcations and chaos in conjunction with even the simplest amplitudelimiting 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 BabcockLeighton 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 zerodimensional dynamical system model originally developed by WilmotSmith 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 spatiallyextended BabcockLeighton 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 modelbased scenarios put forth to explain solar cycle variability on timescales much longer than the period of the primary magnetic activity cycle.
2 BabcockLeighton dynamos
2.1 The BabcockLeighton mechanism
The BabcockLeighton 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 crossequatorial 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; WilmotSmith et al., 2006; Jouve et al., 2010).
2.2 The nonlinearity
A key aspect of the BabcockLeighton 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 Lawlike 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 EW tilt essential to the BabcockLeighton 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. DasiEspuig 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 WilmotSmith et al. (2006) and also used by Albert et al. (2021). The BabcockLeighton source term S(B) for the poloidal magnetic component is written as:
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:
where erf(x) is the error function. The parameters w_{1} and w_{2} control the sharpness of the transition at the lower (B_{1}) and upper (B_{2}) thresholds. Figure 1 shows three examples for specific parameter choices. Note that if and B_{1} > 0, f(B) → 0 as B → 0, the dynamo will no longer be selfexcited, 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 w_{1} = w_{2} ≃ 1 yields sharp thresholds, while larger values lead to smoother transitions. 
2.3 The axisymmetric meanfield dynamo equations
Dynamo models of the solar cycle can be constructed within the general framework of meanfield electrodynamics, by replacing the αeffect a suitable source term modeling the operation of the BabcockLeighton mechanism (see, e.g., Charbonneau, 2020, Sect. 5.4). Working in spherical polar coordinates (r, θ, ϕ) the largescale axisymmetric magnetic field B(r, θ, t) is expressed as B = ∇ × [A(r, θ, t)ê_{ϕ}] + B(r, θ, t)ê_{ϕ}, which allows to separate the induction equation into:
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 BabcockLeighton 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 BabcockLeighton nonlinearity and time delays, as well as a computationally efficient exploration of parameter space.
WilmotSmith et al.’s (2006) dynamo model is a heuristic spatial averaging of the meanfield 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):
where L is a length scale over which the angular velocity varies by Ω_{0}, thus setting the magnitude of the rotational shear. As detailed in WilmotSmith et al. (2006), the parameters T_{0} and T_{1} 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 secondorder nonlinear ODE:
with the ratio q between the total time delay and diffusion time,
determining whether the dynamo is “advectiondominated” (q < 1) or “diffusiondominated” (q > 1). The dynamo number N for this model is given by:
This (relatively) simple dynamical system model thus neatly captures, in simplified form, the nonlocality and nonlinearity of the BabcockLeighton mechanism, as well as the timedelay characterizing the coupling of the two spatiallysegregated source regions.
Although equation (8) looks straightforward, the embodied time delay q makes it nontrivial to solve numerically. After testing a few packages for timedelay ODEs, we selected the Python package Pydelay^{1} which, we found, offers a robust balance between accuracy and runtime performance, at least for our ODE (8). Empirical bifurcations diagrams are constructed by integrating equation (8) in time from some nonzero initial condition, removing the transient phase, collecting all B^{2} 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/N_{crit} = 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; advectiondominated regime), a single dynamo branch emerges from the first bifurcation and persists all the way to the strongly supercritical regime. The coexisting 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_{1} ≲ B_{1}. 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 diffusiondominated regime), perioddoubling 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 finitesized 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 singlyperiodic branch.
The coexistence of multiple stable dynamo branches at a given dynamo number N has other interesting consequences if lowamplitude 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^{2}(t), the same amplitude temporally lagged by q, and the time derivative of B^{2}(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 lowamplitude stochastic noise.
Figure 3 (A) Phase space trajectories for two 0D dynamo solutions with N/N_{crit} = 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^{2} in this model) in a N/N_{crit} = 2.99 solutions akin to Figure 3A, but now with zeromean additive noise uniformly distributed in δα/α_{0} ∈ [−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 WilmotSmith 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 solarlike internal differential rotation profile, as well as spatiotemporal smoothing by magnetic diffusion. The limited set of twodimensional meanfieldlike BabcockLeighton 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 meanfieldlike dynamo model described in Section 2 herein.
4 Supercritical solutions in 2D meanfieldlike dynamo model
We now turn to numerical solutions of the twodimensional 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^{2}/η_{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:
where α_{0} sets the magnitude of the BabcockLeighton source term (viz. Eq. (1)), and u_{0} is the surface midlatitude meridional flow speed, fixed at a solarlike value u_{0} = 10^{3} cm s^{−1} in all simulations. The dynamo number (D) for this model is given by the product of C_{α} and C_{Ω}:
effectively identical to equation (10) for the 0D model. The parameterized meridional circulation and differential rotation profiles are solarlike 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 WilmotSmith et al. (2006), so we opted to restrict our parameter space to the BabcockLeighton nonlinearity, fixing B_{1} = 100, w_{1} = 10 and varying B_{2}, w_{2}. 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_{2}/B_{1} ≃ 10 and w_{2}/w_{2} ≫ 1. We compute solutions for three values of envelope magnetic diffusivity, log(η_{T}) = 11.0, 11.6, and 12.0 (η_{T} in cm^{2} s^{−1}), spanning the transition from advectiondominated to diffusiondominated 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 volumeintegrated 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 BabcockLeighton 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 (D_{crit}) 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. 
Defining parameters for solution sequences.
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 Hopftype 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 perioddoubling cascade (120 sequences). In the more strongly supercritical regime (D/D_{crit} ≳ 5, say), most sequences revert to a highamplitude single branch, undergoing further perioddoubling 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 quasiperiodic modulation of periods of ≈0.18 τ_{D} superimposing itself on a chaotic modulation.
Figure 5 Time series of magnetic energy (red and blue) and midlatitude 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 BabcockLeighton 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 halfcycle 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 timelatitude diagram of the toroidal magnetic component, extracted at the coreenvelope interface r/R = 0.7 (top diagrams in Figs. 6A, 6B, and 6C). The B110 advectiondominated dynamo solution (Figs. 6A) powers itself primarily from the lowlatitude 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 highlatitude 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 highlatitude mode. The resulting beating phenomenon causes the periodic dropin magnetic energy so conspicuous on Figure 5C (top time series), as well as on the timelatitudediagram 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 timelatitude “butterfly” diagrams for three sample solutions taken from Figure 5, all using the same BabcockLeighton 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 timelatitude diagram of the toroidal magnetic component, a timelatitude diagram of the magnetic part of the BabcockLeighton 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_{2} of the BabcockLeighton 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_{2}. 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 GnevyshevOhllike 2cycle 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 BabcockLeighton nonlinearity f(B) along fixedlatitude temporal cuts along the two timelatitude 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 nonlinearity 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 halfcycles, with cycle maxima occurring while f(B) ≃ 0. Note in particular how, on Figure 7B, the solutions extend far in the fullyquenched extent of the BabcockLeighton 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_{1}, B_{2}] interval of the BabcockLeighton nonlinearity, and a second “strong” mode spending a significant fraction of its magnetic cycle in the stronglyquenched range (B_{ϕ} > B_{2}) of the nonlinearity.
Figure 7 Sampling of the BabcockLeighton nonlinearity f(B) for the two timelatitude 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 nonlinearity curve (dotted line). Note the marked magnetic polarity asymmetry in panel B. 
5 Intermittency and Grand Maxima
Having demonstrated the coexistence 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 (α_{0} in Eq. (1)). This fluctuation is extracted from a zeromean 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} = ±0.2 is introduced in the BabcockLeighton source term, with coherence time τ/τ_{D} = 1.22 × 10^{−4}, 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 midlatitude and at R = 0.7 R_{⊙} (red). (B) Time series of smoothed fullcycle amplitude A_{k} (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 A_{k}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 A_{k} histogram, plotted in (C). Panels D and E show histograms of strong mode durations and interevent 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 14641 averaging filter. The resulting smoothed time series displays the evolution of fullcycle 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 polaritysymmetric 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 interevent 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 BabcockLeighton 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 WilmotSmith et al. (2006), and the spatiallyextended twodimensional meanfieldlike model of Charbonneau et al. (2005). Both models are kinematic, in the sense that they implement temporally steady largescale differential rotation and meridional circulation, and also use the same nonlocal amplitudelimiting 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 amplitudelimiting 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 coexistence, 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 BabcockLeighton mechanism. The primary dynamo branch, emerging as the critical dynamo number is exceeded, typically undergoes a perioddoubling 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 twodimensional models, we sometimes observe quasiperiodic cycle amplitude modulation on long timescales, i.e., many times the primary cycle period; this is typically observed in nonkinematic dynamos where magnetic backreaction quenches largescale 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 reestablish the largescale flows after the magnetic amplitude has fallen in response to the quenching of the largescale 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 coexisting at a given dynamo number. Working also with the WilmotSmith 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 Minimalike 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 WilmotSmith 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 Maunderlike 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 coexist over an extended range of dynamo numbers (viz. Fig. 2B) and each retain a finitesized basin of attraction. It would be particularly interesting to introduce nonkinematic effects in the 2D model, to ascertain to what extent such effect could lead to additional longtimescale 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/052782018; PC), the Fonds de Recherche du Québec Nature et Technologies (graduate fellowship number 287277; CT), and a HydroQué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. FerrizMas on the subject matter of this paper. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.
References
 Albert C, FerrizMas A, Gaia F, Ulzega S. 2021. Can stochastic resonance explain recurrence of Grand Minima? Astrophys J Lett 916(2): L9. https://doi.org/10.3847/20418213/ac0fd6. [CrossRef] [Google Scholar]
 Arlt R, Weiss N. 2014. Solar activity in the past and the chaotic behaviour of the dynamo. Space Sci. Rev. 186(1–4): 525–533. https://doi.org/10.1007/s1121401400635. [NASA ADS] [CrossRef] [Google Scholar]
 Babcock H. 1961. The topology of the sun’s magnetic field and the 22year cycle. Astrophys J 133: 572–589. https://doi.org/10.1086/147060. [CrossRef] [Google Scholar]
 Beer J, Tobias SM, Weiss NO. 2018. On longterm modulation of the Sun’s magnetic cycle. Mon Notices Royal Astron Soc 473(2): 1596–1602. https://doi.org/10.1093/mnras/stx2337. [CrossRef] [Google Scholar]
 Brooke J, Moss D, Phillips A. 2002. Deep minima in stellar dynamos. A&A 395: 1013–1022. https://doi.org/10.1051/00046361:20021320. [CrossRef] [EDP Sciences] [Google Scholar]
 Brun AS, Strugarek A, Noraz Q, Perri B, Varela J, Augustson K, Charbonneau P, Toomre J. 2022. Powering stellar magnetism: energy transfers in cyclic dynamos of sunlike stars. Astrophys J 926(1): 21. https://doi.org/10.3847/15384357/ac469b. [CrossRef] [Google Scholar]
 Bushby P. 2006. Zonal flows and grand minima in a solar dynamo model. Mon. Not. R. Astron. Soc. 371: 772–780. https://doi.org/10.1111/j.13652966.2006.10706.x. [CrossRef] [Google Scholar]
 Cameron RH, Schüssler M. 2017. Understanding solar cycle variability. Astrophys J 843: 111. https://doi.org/10.3847/15384357/aa767a. [CrossRef] [Google Scholar]
 Charbonneau P. 2001. Multiperiodicity, chaos, and intermittency in a reduced model of the solar cycle. Solar Phys 199: 385–404. https://doi.org/10.1023/A:1010387509792. [CrossRef] [Google Scholar]
 Charbonneau P. 2020. Dynamo models of the solar cycle. Living Rev in Sol Phys 17(1): 4. https://doi.org/10.1007/s41116020000256. [CrossRef] [Google Scholar]
 Charbonneau P, Barlet G. 2011. The dynamo basis of solar cycle precursor schemes. J Atmos Sol Terr Phys 73(2–3): 198–206. https://doi.org/10.1016/j.jastp.2009.12.020. [CrossRef] [Google Scholar]
 Charbonneau P, BlaisLaurier G, StJean C. 2004. Intermittency and phase persistence in a BabcockLeighton model of the solar cycle. Astrophys J 616: L183–L186. https://doi.org/10.1086/426897. [CrossRef] [Google Scholar]
 Charbonneau P, StJean C, Zacharias P. 2005. Fluctuations in BabcockLeighton dynamos. I. Period doubling and transition to chaos. Astrophys J 619: 613–622. https://doi.org/10.1086/426385. [CrossRef] [Google Scholar]
 DasiEspuig M, Solanki SK, Krivova NA, Cameron R, Peñuela T. 2010. Sunspot group tilt angles and the strength of the solar cycle. A&A 518: A7. https://doi.org/10.1051/00046361/201014301. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dikpati M, Charbonneau P. 1999. A BabcockLeighton flux transport dynamo with solarlike differential rotation. Astrophys J 518(1): 508–520. https://doi.org/10.1086/307269. [CrossRef] [Google Scholar]
 Durney BR. 2000. On the differences between odd and even solar cycles. Solar Phys. 196: 421–426. https://doi.org/10.1023/A:1005285315323. [CrossRef] [Google Scholar]
 Fan Y. 2021. Magnetic fields in the solar convection zone. Living Rev in Sol Phys 18(1): 5. https://doi.org/10.1007/s41116021000312. [CrossRef] [Google Scholar]
 Hale GE, Ellerman F, Nicholson SB, Joy AH. 1919. The magnetic polarity of sunspots. Astrophys J 49: 153. https://doi.org/10.1086/142452. [CrossRef] [Google Scholar]
 Hoyng P, Schmitt D, Teuben LJW. 1994. The effect of random alphafluctuations and the global properties of the solar magnetic field. A&A 289: 265–278. [Google Scholar]
 Jha BK, Karak BB, Mandal S, Banerjee D. 2020. Magnetic field dependence of bipolar magnetic region tilts on the sun: indication of tilt quenching. Astrophys J Lett 889(1): L19. https://doi.org/10.3847/20418213/ab665c. [CrossRef] [Google Scholar]
 Jouve L, Proctor MRE, Lesur G. 2010. Buoyancyinduced time delays in BabcockLeighton fluxtransport dynamo models. A&A 519: A68. https://doi.org/10.1051/00046361/201014455. [CrossRef] [EDP Sciences] [Google Scholar]
 Karak BB. 2023. Models for the longterm variations of solar activity. Living Rev Sol Phys 20: 3. https://doi.org/10.1007/s4111602300037y. [CrossRef] [Google Scholar]
 Küker M, Arlt R, Rüdiger R. 1999. The Maunder minimum as due to magnetic Λquenching. A&A 343: 977–982. [Google Scholar]
 Leighton R. 1964. Transport of magnetic fields on the sun. Astrophys J 140: 1547–1562. https://doi.org/10.1086/148058. [CrossRef] [Google Scholar]
 McClintock BH, Norton AA. 2013. Recovering Joy’s law as a function of solar cycle, hemisphere, and longitude. Sol Phys 287(1–2): 215–227. https://doi.org/10.1007/s1120701303380. [CrossRef] [Google Scholar]
 Moss D, Brandenburg A, Tavakol R, Tuominen I. 1992. Stochastic effects in meanfield dynamos. A&A 265: 843–849. [Google Scholar]
 Ossendrijver A, Hoyng P, Schmitt D. 1996. Stochastic excitation and memory of the solar dynamo. A&A 313: 938–948. [Google Scholar]
 Ossendrijver MAJH. 2000. The dynamo effect of magnetic flux tubes. A&A 359: 1205–1210. [Google Scholar]
 Schatten KH, Scherrer PH, Svalgaard L, Wilcox JM. 1978. Using dynamo theory to predict the sunspot number during solar cycle 21. Geophys Res Lett 5(5): 411–414. https://doi.org/10.1029/GL005i005p00411. [CrossRef] [Google Scholar]
 Simard C, Charbonneau P. 2020. Grand Minima in a spherical nonkinematic α_{2} Ω meanfield dynamo model. J Space Weather Space Clim 10: 9. https://doi.org/10.1051/swsc/2020006. [CrossRef] [EDP Sciences] [Google Scholar]
 Svalgaard L, Cliver EW, Kamide Y. 2005. Sunspot cycle 24: Smallest cycle in 100 years? Geophys Res Lett 32(1): L01104. https://doi.org/10.1029/2004GL021664. [Google Scholar]
 Tobias SM. 1997. The solar cycle: parity interactions and amplitude modulation. A&A 322: 1007–1017. [Google Scholar]
 Tobias SM, Weiss NO, Kirk V. 1995. Chaotically modulated stellar dynamos. Mon Notices Royal Astron Soc 273(4): 1150–1166. https://doi.org/10.1093/mnras/273.4.1150. [CrossRef] [Google Scholar]
 Tripathi B, Nandy D, Banerjee S. 2021. Stellar midlife crisis: subcritical magnetic dynamos of solarlike stars and the breakdown of gyrochronology. Mon Notices Royal Astron Soc Lett 506(1): L50–L54. https://doi.org/10.1093/mnrasl/slab035. [CrossRef] [Google Scholar]
 Usoskin IG. 2023. A history of solar activity over millennia. Liv Rev Sol Phys 20(1): 2. https://doi.org/10.1007/s4111602300036z. [CrossRef] [Google Scholar]
 Usoskin IG, Hulot G, Gallet Y, Roth R, Licht A, Joos F, Kovaltsov GA, Thébault E, Khokhlov A. 2014. Evidence for distinct modes of solar activity. A&A 562: L10. https://doi.org/10.1051/00046361/201423391. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin IG, Sokoloff D, Moss D. 2009. Grand Minima of solar activity and the meanfield dynamo. Sol Phys 254(2): 345–355. https://doi.org/10.1007/s1120700892936. [CrossRef] [Google Scholar]
 Weiss NO, Tobias SM. 2016. Supermodulation of the Sun’s magnetic activity: the effects of symmetry changes. Mon Notices Royal Astron Soc 456: 2654–2661. https://doi.org/10.1093/mnras/stv2769. [CrossRef] [Google Scholar]
 WilmotSmith A, Nandy D, Hornig G, Martens P. 2006. A time delay model for solar and stellar dynamos. Astrophys. J. 652: 696–708. https://doi.org/10.1086/508013. [CrossRef] [Google Scholar]
 Yoshimura H. 1978. Nonlinear astrophysical dynamos: multipleperiod dynamo wave oscillations and longterm modulations of the 22 year solar cycle. Astrophys J 226: 706–719. https://doi.org/10.1086/156653. [CrossRef] [Google Scholar]
Cite this article as: Thibeault C, Miara L & Charbonneau P 2023. Nonlinearity, time delay, and Grand Maxima in supercritical BabcockLeighton dynamos. J. Space Weather Space Clim. 13, 32. https://doi.org/10.1051/swsc/2023029.
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 w_{1} = w_{2} ≃ 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/N_{crit} = 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/N_{crit} = 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 (D_{crit}) 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 midlatitude 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 BabcockLeighton 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 timelatitude “butterfly” diagrams for three sample solutions taken from Figure 5, all using the same BabcockLeighton 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 BabcockLeighton nonlinearity f(B) for the two timelatitude 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 nonlinearity 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 midlatitude and at R = 0.7 R_{⊙} (red). (B) Time series of smoothed fullcycle amplitude A_{k} (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 A_{k}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 A_{k} histogram, plotted in (C). Panels D and E show histograms of strong mode durations and interevent waiting times, equivalent to the duration distribution for normal dynamo operation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.