Issue 
J. Space Weather Space Clim.
Volume 10, 2020
Topical Issue  Space climate: The past and future of solar activity



Article Number  9  
Number of page(s)  15  
DOI  https://doi.org/10.1051/swsc/2020006  
Published online  16 March 2020 
Research Article
Grand Minima in a spherical nonkinematic α^{2}Ω meanfield dynamo model
Département de Physique, Université de Montréal, C.P. 6128, Montréal, QC, Canada
^{*} Corresponding author: corinne.simard2@canada.ca
Received:
27
November
2019
Accepted:
22
January
2020
We present a nonkinematic axisymetric α^{2}Ω meanfield dynamo model in which the complete αtensor and mean differential rotation profile are both extracted from a global magnetohydrodynamical simulation of solar convection producing cycling largescale magnetic fields. The nonlinear backreaction of the Lorentz force on differential rotation is the only amplitudelimiting mechanism introduced in the meanfield model. We compare and contrast the amplitude modulation patterns characterizing this meanfield dynamo, to those already wellstudied in the context of nonkinematic αΩ models using a scalar αeffect. As in the latter, we find that large quasiperiodic modulation of the primary cycle are produced at low magnetic Prandtl number (Pm), with the ratio of modulation period to the primary cycle period scaling inversely with Pm. The variations of differential rotation remain well within the bounds set by observed solar torsional oscillations. In this lowPm regime, moderately supercritical solutions can also exhibit aperiodic Maunder Minimumlike periods of strongly reduced cycle amplitude. The interevent waiting time distribution is approximately exponential, in agreement with solar activity reconstructions based on cosmogenic radioisotopes. Secular variations in lowlatitude surface differential rotation during Grand Minima, as compared to epochs of normal cyclic behavior, are commensurate in amplitude with historical inferences based on sunspot drawings. Our modeling results suggest that the low levels of observed variations in the solar differential rotation in the course of the activity cycle may nonetheless contribute to, or perhaps even dominate, the regulation of the magnetic cycle amplitude.
Key words: Solar cycle / solar dynamo / Grand Minima
© C. Simard & P. Charbonneau, Published by EDP Sciences 2020
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
The past decade has witnessed remarkable progress in the design of global magnetohydrodynamical (MHD) simulations of solar convection and dynamo action. Many different such simulations, relying on distinct physical and algorithmic designs, have now managed to produce magnetic fields wellorganized on large spatial scales, and undergoing moreorless regular polarity reversals resembling to some extent solar behavior (Brown et al., 2010, 2011; Ghizaru et al., 2010; Racine et al., 2011; Käpylä et al., 2012, 2013; Masada et al., 2013; Nelson et al., 2013; Fan & Fang, 2014; Passos & Charbonneau, 2014; Augustson et al., 2015; Simitev et al., 2015; Duarte et al., 2016; Guerrero et al., 2016; Hotta et al., 2016; Käpylä et al., 2016; Strugarek et al., 2018; Guerrero et al., 2019). While some of these simulations are showing occasional episodes of suppressed magnetic activity perhaps analogous to the seventeenth century Maunder Minimum (e.g., Augustson et al., 2015), even the most temporally extended simulation runs carried out yet (e.g., Passos & Charbonneau, 2014; Käpylä et al., 2016) remain akin to “case studies” and cannot be used to systematically assess the occurrence frequency of cycle fluctuation patterns developing on timescales much longer than the primary magnetic cycle period. Toward this end, geometrically and/or dynamicallysimplified dynamo models remain the most practical approach.
Solar cycle models based on meanfield electrodynamics have indeed been used for many years to study longterm cycle variability (see Sect. 5 in Charbonneau, 2010, and references therein). Most often such models solve only for the largescale, axisymmetric component of the solar magnetic field, with the effects of smallscale, turbulent convective motions subsumed in coefficients capturing the effects of turbulent induction and dissipation. Specifying the form of these coefficients, and their dependencies on rotation rate, magnetic field strength, etc, becomes a crucial aspect of this approach. Here global MHD simulations can provide useful guidance, as it is possible to extract the turbulent electromotive force from the output of such simulations, and thus calculate the αtensor and turbulent diffusivity tensor. The former is a key component of many solar cycle models, as it provides a source term for the largescale poloidal magnetic field component, which allows to circumvent Cowling’s theorem. The isotropic part of the latter is an ubiquitous ingredients even in solar cycle models that rely on other inductive effects, most notably perhaps the BabcockLeighton mechanism of dipole regeneration through the surface decay of active regions.
The extraction of the meanfield tensors has been carried out in a subset of the aforecited MHD simulations (Racine et al., 2011; Augustson et al., 2015; Simard et al., 2016; Warnecke et al., 2018). The corresponding simulations are distinct in overall design, treatment of smallscales, parameter regimes, and the techniques used to extract the meanfield coefficients also differ. While results differ at various levels, some common features nonetheless emerge. Most notably perhaps, the αtensor is full, with offdiagonal components of magnitude comparable to diagonal elements. The antisymmetric component of the tensor is also of comparable magnitude to the offdiagonal symmetric components, indicating that turbulent pumping represents a significant transport mechanism for the largescale magnetic field.
In contrast, meanfield dynamo models of the solar cycle are often designed in the socalled αΩ regime: the poloidal largescale magnetic component is generated by the αeffect, while the toroidal magnetic component grows only through the shearing of the poloidal component by differential rotation. In such models the αtensor is thus effectively reduced to its α_{ϕϕ} component, and turbulent diffusivity is deemed isotropic, reducing to a scalar quantity. This αΩ approximation is usually justified through orderormagnitude estimates of the relative importance of turbulent induction and differential rotation shear.
To the extent that the aforecited global MHD simulations are applicable to the sun and stars, the full αtensor extracted from these simulations suggests that dynamos therein are operating in what is usually dubbed the α^{2}Ω regime: both shearing by differential rotation and turbulent inductive effects contribute to the production of the largescale toroidal magnetic component. Even in the kinematic regime, in which largescale flows are assumed steady, α^{2}Ω meanfield models have received far less attention than their αΩ counterparts. Some wellknown solarlike properties of αΩ models do carry over to the α^{2}Ω regime, notably the generation of dynamo waves (e.g., Choudhuri, 1990; Moss et al., 1991; Ruediger & Brandenburg, 1995; Charbonneau & MacGregor, 2001).
Hardly studied thus far in the α^{2}Ω context is the dynamical backreaction of the largescale magnetic field produced by dynamo action on the largescale inductive flows, specifically the differential rotation. This socalled MalkusProctor mechanism (Malkus & Proctor, 1975) has received a lot of attention in the context of αΩ models, as it has been shown capable to produce long timescale modulation of the primary cycle, including recurrent epochs of very low cycle amplitudes, reminiscent of Maunder Minimumlike episodes of suppressed surface magnetic activity (Tobias, 1996, 1997; Küker et al., 1999; Brooke et al., 2002; Bushby, 2006). This long timescale modulation is associated with the suppression of differential rotation by the Lorentz force, which eventually leads to a quenching of the dynamo. The largescale magnetic field then dissipates, eventually allowing the reestablishment of differential rotation, and renewed dynamo action. Investigations of these effects in geometrically simple slab configurations has evidenced rich dynamical behavior, characterized by two basic modes of cycle amplitude variation: parity modulation and energy exchange between the field and flow (Tobias, 1996; Knobloch et al., 1998). These finding appear robust, in that their equivalents have been observed in a variety of geometrically more realistic αΩ model setups.
The primary aim of this paper is to extend these studies to the α^{2}Ω regime. Toward this end we make use of the full αtensor extracted from the millenium simulation of Passos & Charbonneau (2014), as described in Simard et al. (2013, 2016). The mathematical formulation and physical input to our nonkinematic meanfield model are described in Section 2. In Section 3, we present a linear, kinematic reference dynamo solution. Representative high and lowPrandtl number (Pm) nonkinematic dynamo solutions are then presented and analyzed in Sections 4 and 5, respectively. Section 6 focuses on the long timescale modulations materializing in the lowPm regime. We close in Section 7 by examining the similarity and differences with nonkinematic meanfield models of the αΩ variety, and seek the origin of the most striking differences in the interactions between the various inductive effects operating in our α^{2}Ω model.
2 A nonkinematic meanfield model
In this section, we briefly summarize the mathematical development of our nonkinematic meanfield model. Our formulation of α^{2}Ω axisymmetric meanfield dynamo equations is in itself entirely conventional, the model’s novelty lying with the inclusion of dynamical backreaction on the differential rotation and inclusion of an αtensor extracted from a global MHD simulation of solar convection.
2.1 Meanfield dynamo equation
Following the usual meanfied electrodynamics Ansatz, we introduce a scale separation of the flow and magnetic field, their substitution into the MHD induction equation leading straightforwardly, upon averaging, to the meanfield field dynamo equation:(1)where(2)is the turbulent electromotive force (emf) and η is the magnetic diffusivity (see e.g., Moffatt, 1978; Krause & Rädler, 1980; Ossendrijver, 2003; Brandenburg & Subramanian, 2005; Charbonneau, 2010). The prime quantities in the above expression represent the smallscale, fluctuating parts of the total flow and magnetic field. Closure is achieved by developing ℰ in terms of the mean magnetic field :(3)where the tensor a and b appearing in this expression are assumed to depend only on the statistical properties of the smallscale flow and field (see, e.g., Moffatt, 1978). Truncation of the higher order derivatives is justified provided a good separation of scales exists between the fluctuating velocity and magnetic field on one hand, and the largescale magnetic and flow field on the other. The first term in the expansion involves a ranktwo tensor capturing (among other effects) the socalled αeffect, which arises from the symmetric part of the a tensor and can act as a source term in the meanfield equations (5) and (6) (see Rädler, 1980; Schrinner et al., 2007; Simard et al., 2016).
The axisymmetric form of the meanfield dynamo equation (1) is obtained by separating the mean magnetic field 〈B〉into a toroidal and a vector potential components:(4)
Inserting this expression into the meanfield dynamo equation yields the two evolution equations:(5)(6)where ϖ = rsinθ, and Ω_{t} is the total rotational angular velocity. This is the only contribution to the largescale flow in equation (1), i.e., we do not include any largescale meridional circulation in the model. Note that here the mean electromotive force ℰ is retained in equation (6), which yields the socalled α^{2}Ω meanfield dynamo model (see Simard et al., 2013 for more details).
For the purpose of numerical solution, it is convenient to express lengths in units of the solar radius R and time in units of the magnetic dissipation time τ = R^{2}/η_{0}. This leads to the appearance of two dimensionless numbers measuring the inductive strength of the meanelectromotive force and largescale shear, respectively:(7)where α_{0}, η_{0}, and Ω_{0} represent typical values of the corresponding unsubscripted quantities.
2.2 Magnetic backreaction
With the magnetic field generated by the dynamo is associated a Lorentz force opposing the inductive largescale flows. Since the Sun’s convection zone is a strongly turbulent environment, this dynamical backreaction can operate at all spatial scales. Suppression of the smallscale Reynolds stresses powering differential is known as Λquenching and is a potentially efficient mechanism to quench and/or modulate the magnetic cycle (see, e.g., Kitchatinov et al., 1994; Rempel, 2006). In the present nonkinematic dynamo model, we only include the backreaction of the mean magnetic field on the largescale flow, also know as the MalkusProctor mechanism (Malkus & Proctor, 1975). This requires solving the ϕcomponent of the equations of motion, which in principle must include the Reynolds stresses powering differential rotation (see Kitchatinov & Rüdiger, 2005, and references therein). Following Tobias (1996), we circumvent the specification of Reynolds stresses by separating the angular velocity into a temporally steady component Ω(r, θ), and a timedependent contribution Ω′(r, θ, 〈B〉(r, θ, t)) driven by the Lorentz force:(8)
Underlying this separation is the implicit assumption that the Reynolds stresses powering the steady part of the differential rotation remain unaffected by the largescale magnetic field produced by the dynamo, so that Ω represents a stationary solution of the unmagnetized momentum equation. The governing equation for the timevarying contribution Ω′ then becomes:(9)where ν is the viscosity and ρ(r) the density. Again, by developing the mean magnetic field into axisymmetric poloidal (A) and toroidal (B) components, the dimensionless form of equation (9) can be written as:(10)
Two new dimensionless parameters have appeared: the magnetic Prandtl number, Pm = ν/η sets the relative magnitudes of viscous (ν) versus Ohmic (η) dissipation, and the parameter measures the influence of Lorentz force. The coupled axisymmetric dynamo equations (5) and (6) being linear in A, B, the parameter Λ thus sets the scale of the magnetic field amplitude in the model.
2.3 A nonkinematic α^{2}Ω meanfield dynamo model
Our nonkinematic meanfield dynamo model is formulated in spherical polar coordinates (r, θ, ϕ) and solves for the axisymmetric (∂/∂ϕ ≡ 0) vector potential A(r, θ, t), toroidal magnetic field B(r, θ, t) and angular velocity deviation Ω′(r, θ, t) in a full meridional plane restricted in radius to the interval 0.5 ≤ r/R ≤ 1.0. The solution is matched to a potential field in r/R > 1, and the bottom boundary is rigid (Ω′ = 0) and perfectly conducting (A = 0 and B = 0). The nonlinearly coupled system of dimensionless PDEs (5), (6) and (10) is solved using the finite elementbased dynamo model described in Simard et al. (2013), extended to include the Lorentz force backreaction on the zonal flow, as embodied in equation (10).
The steady component Ω(r, θ) is extracted from the output of a purely hydrodynamical EUlerianLAGrangian (EULAG) simulation of solar convection, specifically an unmagnetized version of the millenium simulation discussed in Passos & Charbonneau (2014). The corresponding angular velocity profile is displayed in Figure 1A. It includes a number of solarlike features, including a rapidly rotating equatorial region and slowly rotating poles, and a thin tachoclinelike rotational shear layer immediately beneath the coreenvelope interface, indicated in Figure 1A by the white dashed line. Nonsolar features at low latitudes include alignment of the isocontours with the rotation axis, lowlatitude radial shear peaking within the convective envelope rather than at its base, and lack of a nearsurface shear layer. The largescale meridional flow generated in this simulation is weak, shows strong temporal fluctuations, and lacks a welldefined spatial organization (Racine et al., 2011). We therefore opt to exclude it from our α^{2}Ω meanfield dynamo model.
Fig. 1 (A) Steady angular velocity component (see Eq. (8)) as extracted from a hydrodynamical EULAG simulation, plotted as isocontours of rotational frequency. (B) and (C) Total angular velocity at epoch of Grand Maximum and Minimum, respectively, as produced by our nonkinematic α^{2}Ω meanfield dynamo model (specifically, simulation 4 in Table 1). The interface between the convection zone and underlying radiative zone is indicated by the white dashed line. 
As in Simard et al. (2013), our meanfield dynamo model use an αtensor extracted from a reference EULAGMHD simulation of solar convection (Passos & Charbonneau, 2014) through a leastsquare fit to equation (3), as described in Racine et al. (2011) and Simard et al. (2016). This yields the full 3 × 3 αtensor, in which offdiagonal components are here roughly of the same order of magnitude as the diagonal components (see Fig. 1 in Simard et al., 2013).
The (turbulent) diffusivity η appearing in equations (5) and (6) corresponds to the isotropic part of the turbulent diffusivity tensor embodied in the tensor b on the RHS of equation (3). We use values in the range 2–8 × 10^{11} cm^{2} s^{−1}, commensurate with the values extracted from our reference EULAGMHD simulation (see Simard et al., 2016, Fig. 6). This value is held fixed throughout our convection zone, and then smoothly ramps down to a much smaller value below across a thin layer coinciding with the coreenvelope interface. This variations is parametrized via an error function profile (see Eq. (14) in Simard et al., 2013). All dynamo solutions reported upon in what follow use a coretoenvelope diffusivity ratio of 10^{−2}.
The density ρ(r) in equation (10) follows the same polytropic profile as used in the EULAGMHD from which the αtensor and mean differential rotation are extracted. Likewise, the coreenvelope interface is set at r/R = 0.718. We set a floor on the density at ρ = 0.021 kg m^{−3} in the subsurface layers, in order to help maintain numerical stability in the low Pm regime.
The backreaction of the Lorentz force on the differential rotation, as described by equation (9), is the only nonlinearity present in the model. In order to capture as unambiguously as possible its impact on the dynamo, no additional mechanism (such as αquenching) is introduced to achieve amplitude saturation. No attempt is made to alter the model setup so as to produce solarlooking butterfly diagrams. In fact, we would have very little leeway to do so, since our αtensor component and steady part of differential rotation are all taken as extracted from the parent numerical simulations.
As shown in what follows, this nonkinematic α^{2}Ω meanfield model can generate a rich spectrum of dynamo behavior, including parity modulation, rotational torsional oscillations, and amplitude modulation on long time scales. In the remainder of this paper we discuss in some detail a small subset of representative dynamo solutions that exemplifies these behaviors. Although all simulations make use of the αtensor and differential rotation profile extracted from the EULAGMHD millenium simulation, we emphasize that our primary aim here is not to produce a close meanfield equivalent, but rather to explore the range of possible dynamical behavior associated with magnetic backreaction on differential rotation in a meanfield α^{2}Ω dynamo formulated in spherical geometry. Table 1 lists the defining parameters of the set of representative simulations discussed in the remainder of the paper.
Defining dimensionless parameters for the reference dynamo solutions discussed in the text.
3 Simulation 0: a kinematic, linear solution
As background to the analyses and discussion to follow, we first consider a representative kinematic solution in which the backreaction of the magnetic field on differential rotation has been turned off, by setting Λ = 0 in equation (10). The problem then reduces to solving the linear but coupled system defined by equations (5) and (6), with Ω′ set to zero. Holding C_{Ω} = 35,000 and varying C_{α}, the critical dynamo number is found to be D_{crit} ≡ C_{α} × C_{Ω} = 87,500, for C_{α} ≃ 2.5. Figure 2 displays a moderately supercritical solution (C_{α} = 15, C_{Ω} = 35,000), in which exponential growth has been scaled out. The top panel shows a timelatitude diagram of the toroidal component extracted at middepth (r/R = 0.85). Four snapshots in the meridional plane are plotted at bottom, at times corresponding to the dashed lines on the top panel.
Fig. 2 A kinematic α^{2}Ω dynamo solution. Defining parameters are listed in Table 1 (solution run 0). The toroidal field component is plotted as a timelatitude diagram at depth r/R = 0.85 in the top panel, and meridional plane snapshots extracted at times corresponding to the dashed lines on the top panel, and spanning a magnetic halfcycle, are displayed at bottom. 
This kinematic solution shows concentration of the toroidal magnetic field at low latitudes, consistent with the presence of the strongest differential rotation shear at these latitudes. A clear dynamo wave pattern is also present, characteristic of αΩ meanfield models and also materializing in α^{2}Ω models in which differential rotation contributes significantly to toroidal field induction (Choudhuri, 1990; Charbonneau & MacGregor, 2001), which is the case here. The dynamo waves propagate along the isocontours of angular velocity (cf. Fig. 1), in a direction set by the sign of the product of the α_{ϕϕ} tensor component and angular velocity shear. Here at low latitudes we have ∂Ω/∂r > 0 and α_{ϕϕ} > 0, leading to propagation away from the equatorial plane, as per the ParkerYoshimura sign rule.
We calculate the equatorial parity of our dynamo solutions parity P through(11)with(12)(13)where the toroidal magnetic components B_{N}, B_{S} are sampled in the Northern (N) and Southern (S) hemispheres at r/R = 0.85 and ±15° latitude, corresponding to the region of peak toroidal field in the course of the magnetic cycle. Under these definitions, antisymmetric and symmetric parity correspond to P = −1 and P = +1, respectively. For the solution of Figure 2 the preferred equatorial parity is antisymmetric (i.e., dipolelike).
4 Simulation 1: amplitude saturation
We now turn to a nonkinematic solution, otherwise identical to the linear solution displayed in Figure 2. We now solve the nonlinearly coupled system defined by equations (5), (6), and (10), setting Λ = 1 and magnetic Prandtl number, Pm = 2 in the latter. With C_{α} = 35 and C_{Ω} = 35,000, this is a moderately supercritical solution (D/D_{crit} = 6). In this Pm > 1 regime, the second term on the RHS of equation (10) operates on a time scale shorter than the magnetic cycle period, and thus restricts the perturbation of differential rotation to very low amplitudes. As displayed on Figure 3, the magnetic field grows exponentially but rapidly saturates to a cyclic state of constant cycle amplitude and cycle period similar to that characterizing the linear regime for these parameter values. Dynamo wave propagation also closely resembles the linear regime, except for a faint equatoriallypropagating branch at high latitudes, associated with the perturbation of differential rotation. This nonlinearly saturated dynamo remains locked onto antisymmetric parity (blue time series on Fig. 3A), as in its linear, kinematic counterpart (Fig. 2).
Fig. 3 Solution run 1, a nonkinematic dynamo solution in the Pm > 1 regime. Its defining parameters are listed in the second line of Table 1. (A) Time series of normalized magnetic energy (black), integrated shear proxy (red), and equatorial parity (blue), starting from a seed magnetic field and Ω′ = 0 at t = 0. (B) Timelatitude diagram of the toroidal magnetic component extracted at r/R = 0.85, with a meridional plane snapshot at right. (C) Similar timelatitude diagram at r/R = 0.85 and meridional snapshot, but for the angular velocity perturbation normalized to the background angular velocity profile, i.e., Ω′(r, θ, t)/Ω(r, θ). 
Figure 3C shows a timelatitude diagram of the differential rotation perturbation Ω′/Ω, extracted at the same depth r/R = 0.85 as the timelatitude diagram for the toroidal component (Fig. 3B). In response to the initial exponential growth of the magnetic amplitude, rotational torsional oscillations develop at twice the frequency of the magnetic cycle, as expected, eventually leading to saturation of the magnetic amplitude. In this saturation state, torsional oscillations remain at a very low amplitude, peaking at a mere ~0.3% of the steady differential rotation component. From the inductive point of view, however, what matters is the shear, not the rotation rate per se. Accordingly, we define a local measure of the absolute shear as:(14)
A global version of this proxy is constructed by integrating the above quantity over the solution domain at each time step, and dividing by the integration of equation (14) computed with the steady profile Ω(r, θ) replacing Ω′(r, θ, t). The corresponding time series, plotted in red on Figure 3A, thus gives a measure of the relative variation of the global shear. This measure oscillates with the magnetic cycle frequency, remaining nearly in antiphase with the latter, as measured by total magnetic energy (in black); in other words, the rotational shear is globally reduced at magnetic cycle maximum, reflecting the backreaction of the magnetic force on the inductive largescale flow.
5 Simulations 2, 3, and 4: amplitude modulation
While our nonkinematic α^{2}Ω dynamo solutions stabilize on a constantamplitude cyclic behavior at high (>1) values of the magnetic Prandtl number (Pm), large cycle amplitude fluctuations materialize in the Pm < 1 regime. An analogous behavior has been studied extensively in a wide range of nonkinematic αΩ dynamo models, in which the turbulent electromotive force is zeroed out in equation (6). Geometrically, these models have ranged from spatiallyintegrated dynamical systems, through onedimensional cartesian slab models, all the way to axisymmetric models in spherical geometry (see e.g., Tobias, 1996; Knobloch et al., 1998; Küker et al., 1999; Brooke et al., 2002; Bushby, 2006). In particular, Knobloch et al. (1998) have argued that cycle modulation in such nonkinematic models can be divided into two generic classes: type 1 modulation takes place through parity oscillation, i.e., energy is cyclically exchanged between the fundamental symmetric and antisymmetric modes through nonlinear interactions with the large scale flow, without significantly affecting the kinetic energy of the latter; type 2 modulation, in contrast, involves energy exchange between the magnetic and kinetic energy reservoirs, at fixed parity. The coexistence of both types of modulation in geometrically more complex models can then lead to a wide range of modulation patterns.
The aforecited results were obtained in the context of meanfield models of the αΩ variety. We now turn to nonkinematic α^{2}Ω dynamo solutions (solutions 2, 3 and 4 in Table 1) that exhibit behavior akin to the type 1 and type 2 modulations identified by Knobloch et al. (1998).
Figure 4 shows a representative solution exhibiting parity modulation (type 1). The overall format is the same as on Figure 3, except that the portion of the solution plotted is restricted to the saturated phase of the simulation, i.e., it excludes the initial phase of exponential growth to saturation. This solutions shows a short primary cycle of period ≃ 2.5 years, undergoing an amplitude modulation by a factor ≃2 on a timescale ≃30 years. This amplitude modulation is associated with a clear parity modulation, with transition from antisymmetric (P = −1) to symmetric (P = +1) parity leading to a growth of the modulation envelope, the opposite transition driving instead a decrease of the modulation amplitude. Despite variations by a factor of ~2 in magnetic energy, the zonal flows are hardly affected, with the global shear proxy varying by about 10% between peaks and troughs of the modulating envelope.
Fig. 4 Simulation 2, a nonkinematic α^{2}Ω dynamo exhibiting parity modulation. (A) Normalized magnetic energy in black, the shear proxy (Ψ) integrated in the meridional plane in red, and the parity in light blue. (B) Toroidal magnetic field in a timelatitude diagram extracted at 0.85R and (C) the angular velocity fluctuation (Ω′(r, θ, t)/Ω(r, θ)) in the same configuration. The meridional cuts of the toroidal field B and angular velocity perturbation Ω′ plotted at right of the timelatitude diagrams are extracted at t/τ = 2.7, as indicated by the dashdotted line segments on the corresponding timelatitude diagrams. 
As can be seen on Figure 4B, parity modulation leads here to an alternating pattern of hemispheric dominance in the timelatitude diagram of the internal toroidal magnetic component, evidence that the fundamental symmetric and antisymmetric modes reach comparable energy levels in this dynamo solution. With the symmetric mode now building up to high amplitudes, toroidal magnetic fields now extend across the equatorial plane, even though peak amplitudes are still reached around latitudes ±15°.
Although still modest in absolute terms, the differential rotation perturbation, plotted on Figure 4A, now reaches a much larger amplitude than it did in the highPm regime (cf. Fig. 3), but with only a very faint cyclic signal on the timescales of the primary magnetic cycle. Nonetheless, an overall reduction of the differential shear in the dynamo region remains the mechanism saturating the overall amplitude of this supercritical dynamo solution, even though the magnitude of the flow perturbation is at the ~1% level. As in the case of the Pm = 2 solution of Figure 3, comparing the meridional cuts at bottom right to the mean differential rotation profile on Figure 1 reveals immediately that at low latitudes, the zonal flow perturbation is negative (positive) where the angular velocity is highest (lowest), i.e., the magneticallydriven flow perturbation reduces the shear characterizing the background differential rotation profile.
Figure 5 offers a different view of the amplitude modulation pattern, in the form of a trajectory in a threedimensional space defined by parity, magnetic field and the Ω proxy from the solution run 2 presented in Figure 4. In order to have a clean path in the phase space, we extract only the modulation envelope for the toroidal magnetic field to avoid sign change. The result is a welldefined path where the solution does not change much in magnetic amplitude (~25%) or in shearing (~35%) but varies from −1 to +1 in parity. One lap along this attractor lasts ~0.175 diffusion time (see Fig. 4A). In this global phase space, the attractor is a thin, twisted tube bounding the trajectory. The bottom panel on Figure 5 is a Poincaré section computed by plotting the points where the trajectory on the left panel crosses the P = 0 plane, with the colors now coding time (blue → green → yellow → red). While the four crossing regions in the P = 0 plane each cover only a small surface, as shown on the inset within these regions the trajectory is spacefilling, consistent with the expected chaotic nature of the cycle modulation.
Fig. 5 Top panel: 3D phase portrait between parity, toroidal magnetic field and the shear proxy (see Eq. (14)) from solution 2. Each of these three quantities is extracted at r/R = 0.85 and 15° latitude. Only the modulation envelope of those quantities is used for plotting, i.e., the primary cycle is removed. The color code represents the magnitude of the parity from −1 in blue to +1 in red. The bottom panel is a Poincaré section in the P = 0 plane, with the inset focusing on one of the four crossing locations. The trajectory therein is spacefilling, indicating that the modulation is chaotic (see text). 
Figure 6 shows yet another example of a nonkinematic α^{2}Ω dynamo solution (run 3 from Table 1), this time obtained at a much lower value of the magnetic Prandtl number, Pm = 5 × 10^{−3}. Magnetic activity is now peaking at midlatitudes in both hemisphere, with much less transequatorial coupling than on the Pm = 0.08 solution of Figure 4. Variations of differential rotation are again quite small, well below the 1% level, and are developing on the long timescale of the cycle amplitude modulation, without any significant signal at the primary cycle frequency. Variations in parity are again quite large, but mostly reflect a phase drift between hemisphere, setting in during the rising phase of the magnetic cycle modulation.
Fig. 6 Identical in format to Figure 4, but showing now a very low Pm solution (Pm = 5 × 10^{−3}; solution run 3 in Table 1). 
The two solutions of Figures 4 and 6 exhibit a property already wellknown of nonkinematic αΩ meanfield models operating in the lowPm regime, namely their ability to generate periodic amplitude modulations on timescales much longer than the period of the primary magnetic cycle. More specifically, the ratio between the longterm modulation period (P_{2}) and the primary cycle period (P_{1}) is observed to increase when Pm is decreased (Tobias, 1997; Bushby, 2006). A similarly robust trend also emerges from our nonkinematic α^{2}Ω model. This is shown on Figure 7, for a sequence of solutions of varying Pm but otherwise identical to solution 3 (Fig. 6). The period ratio P_{2}/P_{1} is found to scale as ∝Pm^{−1} over an orderofmagnitude in Pm.
Fig. 7 Variation of the ratio between the longterm modulation period (P_{2}) and the primary magnetic cycle period (P_{1}) versus the magnetic Prandtl number (Pm). Except for the latter, all solutions use the same parameter values as solution 3 in Table 1 (see also Fig. 6). The dashed line indicates a Pm^{−1} scaling. 
The solutions of Figures 4 and 6 show quasiperiodic longterm modulations, but more strongly apediodic modulation also develops in many regions of the model’s parameter space, in particular at higher dynamo numbers and small values of Pm. Figure 8 shows an example (solution 4 in Table 1), in the now usual format. Once again magnetic activity is concentrated to low latitudes, and is characterized by parity modulation, hemispheric decouping, and variations of differential rotation reminiscent of type 2 modulation.
Fig. 8 Identical in format to Figure 4, but for a solution exhibiting strongly aperiodic behavior. This is solution run 4 in Table 1. 
In this solution, extended epochs of high amplitude in magnetic energy (such as beginning here at ≃ 9150 years) tend to be of symmetric equatorial parity, while periods during which antisymmetric parity prevails almost always coincide with low magnetic energy.
Interestingly, in all of our nonkinematic α^{2}Ω solutions the variations in angular velocity remain quite small, within the observed level of solar rotational torsional oscillations (Howe, 2009). This remains the case even when large variations of the shear proxy are taking place over the long modulation period, as on Figure 6, in contrast to nonkinematic αΩ model undergoing type 2 modulation in the low Pm regime (Bushby, 2006). The physical origin of this behavioral difference is discussed in Section 7.
6 Characteristics of Grand Minima
The occurrence of extended periods of strongly suppressed activity, of which the 1645–1715 Maunder Minimum remains the archetype (Eddy, 1976), is one of the most intriguing feature of the solar activity record. With only one such event present in the 400year sunspot record, the statistical study of these Grand Minima is better carried out via indirect proxies of solar activity, namely the abundances of cosmogenic radioisotopes in natural archives (see Usoskin et al., 2007; Usoskin, 2017, and references therein). Reconstructions of solar activity based on ^{14}C or ^{10}Be extends back over 10,000 years, and reveal important variations in the overall level of solar activity. While these reconstructions are model dependent to some degree, some robust features do emerge (Damon & Sonett, 1991; Vasiliev & Dergachev, 2002; Peristykh & Damon, 2003; Usoskin & Kovaltsov, 2004; Usoskin et al., 2007). These include a number of Maunder Minimumlike episodes of suppressed activity, as well as periods of much higher than average activity (dubbed Grand Maxima). Because of the wide disparity of timescales involved, meanfield dynamo models remain the primary modeling tool used to explore the underlying physical mechanisms.
Many dynamobased models have been developed as explanatory frameworks for Grand Minima (see Sect. 5 in Charbonneau, 2010). One class of such models is based on intermittency, namely transition between distinct dynamical states mediated by either deterministic or stochastic processes. In this context, stochastically forced dynamos operating close to criticality can produce very solarlike patterns of Grand Minima (see, e.g., Olemskoy & Kitchatinov, 2013; Passos et al., 2014; Cameron & Schüssler, 2017; Inceoglu et al., 2017). Strong amplitude modulation, as exemplified by solutions 3 and 4 above (Figs. 6 and 8) acting in conjunction with a field strength threshold for the destabilization and emergence of sunspotforming magnetic flux ropes, represent another attractive physical explanation for Grand Minima. In such a scenario, the underlying magnetic cycle never stops, but simply falls below the threshold at which sunspot formation is possible. This notion is supported by the ^{10}Be radioisotope record, which shows a sustained cyclic signal across the Maunder Minimum (see Beer et al., 1998; Cliver et al., 1998; Usoskin et al., 2015), presumably driven by a weak, global solar magnetic field undergoing cyclic polarity reversals even in the near absence of sunspots.
In the remainder of this section we present a detailed analysis of a 50,000 long extension of simulation 4 (Table 1 and Fig. 8), in which Grand Minima recur in an aperiodic fashion. In this simulation the epochs of very high modulation amplitude, referred to as “Grand Maxima” in what follows, are of equatorially symmetric parity, and are almost invariably followed by a Grand Minimum. Parity then usually switches to antisymmetric in the recovery to normal cyclic behavior upon exiting Grand Minima. Observational evidence exists for a qualitatively similar behavior occuring at the end of the Maunder Minimum (Ribes & NesmeRibes, 1993; Sokoloff & NesmeRibes, 1994).
6.1 Parity switch across Grand Minima
Figure 9 focuses on a representative Grand Minimum occuring in solution 4. The top panel is again a timelatitude diagram of the toroidal magnetic field, while the four meridional snapshots at bottom show the evolving structure of the complete axisymmetric magnetic field. The transition from symmetric to antisymmetric equatorial parity going in and out of the Grand Minimum is strikingly apparent on the top panel. For this specific Grand Minimum the return to symmetric parity occurs rather rapidly upon recovery to normal cyclic behavior, but in other cases antisymmetric parity can persist for many tens of cycles (cf. the Grand Minimum on Fig. 8).
Fig. 9 Snapshot of an epoch of a grand maximum followed by a grand minimum for solution run 4. Top panel shows the toroidal magnetic field in a timelatitude diagram extracted at r/R = 0.85. The bottom panel shows four meridional planes snapshots extracted from the simulation at times indicated by the dashed lines on the upper panel. The color scale encodes the toroidal magnetic field, and poloidal fieldlines are plotted as solid (dashed) lines for counterclockwise (clockwise) orientation. 
Figure 10 shows a phase portrait of the dynamo solution during the time interval spanned by Figure 9. The solution is plotted as a trajectory in a 3D space defined by the parity, and local values of the toroidal magnetic field and shear proxy (Eq. (14)) evaluated at (r/R, θ) = (0.85, 15°). The direction of time follows the color shading, from blue to red. The parity transition is quite smooth and follows a welldefined oscillatory path in this threedimensional phase space, in which the solution first undergoes a collapse of the magnetic cycle amplitude (dark to light blue) as the shear is reduced, followed by a slow parity transition with marginal growth of the amplitude (green), which picks up once dominantly antisymmetic parity is restored. The growth of the angular velocity perturbation then saturates the dynamo (yellow to orange), which then transits rapidly back to symmetric parity (orange to red).
Fig. 10 3D phase portrait between parity, toroidal magnetic field and shear proxy from solution 4 taken at the same time span as Figure 9 and extracted at 0.85R radius and 15° latitude. The color encodes the advance of time from blue to red. 
The general dependence of cycle amplitude on parity is not restricted to epochs of Grand Minima and Maxima. The distribution of cycle amplitudes for the full simulation is plotted on Figure 11A, and appears roughly symmetrical about its mean, in gaussianlike fashion. However, if one extracts cycles for which the parity is either predominantly dipolelike or quadrupolelike, the two data subsets show markedly distinct distributions. These are plotted on Figure 11B, using only cycles with P ≤ −0.7 (dominantly antisymmetric, in purple) or P ≥ +0.7 (dominantly symmetric, in yellow). These differences reflect the spatially distinct patterns of Lorentz force altering the internal differential rotation under the two parity states of the dynamo.
Fig. 11 (A) Amplitude distribution of the toroidal magnetic field maxima extracted from the primary cycle in simulation run 4. Panel (B) separates these amplitude by parity, considering only cycles for which P ≥ 0.7. The corresponding distribution for symmetric and antisymmetric parity are plotted in yellow and purple, respectively, with their sum in red. 
6.2 Statistics of Grand Minima and Maxima
Figure 12A shows a 11,000 years segment of a ~50,000 yearlong simulation where episodes of grand maxima/minima are observed. Figure 12A illustrates the modulation envelope of normalized total magnetic energy, where the colored regions in red (blue) represent period of grand maxima (minima). This time series is constructed by plotting only the peak amplitude of each individual energy peak associated with the primary (short) cycle (cf. black line on Fig. 8). The threshold values used to delineate epochs of Grand Minima and Maxima from “normal” cyclic behavior, are indicated by the two horizontal dashed lines. These thresholds are set here so that the fraction of simulation time spent in either state is roughly the same (~10%) as inferred from radioisotopebased reconstructions (see Fig. 3 in Usoskin et al., 2007). As in those reconstructions, the oscillation associated with the primary cycle is lost and only the longterm activity modulation trend remains. The fact that this simulation exhibits markedly aperiodic longterm modulation results entirely from the primary cycle operating in the chaotic regime, and of the relatively weak coupling between hemispheres characterizing the parameter regime in which it operates.
Fig. 12 (A) Time series of total magnetic energy in a 11,000 years segment of a 50,000 years realization of simulation 4. Time is expressed both in units diffusion time (bottom axis) and in years (top axis). The oscillations from the primary cycle have been removed, leaving only the longterm modulation envelope. Epochs of grand maxima and minima are indicated in red and blue, respectively. Compare to Figure 3 in Usoskin et al. (2007). The vertical dashed lines delineate the time interval spanned by Figure 8. The four other panels show histograms of (B) Grand Minima durations, (C) waiting times between successive Grand Minima, (D) Grand Maxima durations, and (E) waiting times between successive Grand Maxima (see text). 
Figure 12B plots the distribution of Grand Minima durations, as extracted from the time series of panel (A). This distribution is built from 132 measured Grand Minima, and its low duration end is influenced to some extent by our adopted lower threshold value. Nonetheless, the distribution shows a distinct, broad peak for durations spanning the range 11–23 periods of the primary cycle. That Grand Minima should be characterized by a moderately welldefined mean value is expected here. As discussed already in the context of Figure 10, Grand Minima are triggered when the cycle amplitude becomes high enough to reduce differential rotation close to criticality. The transition back to normal cyclic behavior is controlled by two processes: dissipation of the magnetic field, accompanied by dissipation of the differential rotation perturbation Ω′, the timescale of latter process being controlled by the magnetic Prandtl number, Pm = ν/η on the RHS of equation (10).
Figure 12C plots the waiting time distribution (WTD) for Grand Minima, waiting time being defined as the time elapsed between the end of a Grand Minimum and the onset of the next. The WTD has an exponential shape, indicative of a stationary memoryless Poisson process, i.e., the onset of a Grand Minimum is not influenced by the occurrence of earlier Grand Minima. This is in partial agreement with inferences based on analyses of radioisotope data, which to first order reveal an exponential WTD, but also small deviations from pure exponential behavior ascribed to a tendency for Grand Minima to cluster (see Usoskin et al., 2007, 2009, 2016).
Similar statistical trends are observed for grand maxima, as shown on Figures 12D and 12E. Grand Maxima durations show again a welldefined peak centered on 25–30 primary cycle periods. This characteristic timescale is governed by the time taken for the growing magnetic field to quench differential rotation, and is controlled by the parameter Λ in equation (10). The interevent waiting time for grand maxima is again approximately exponential, indicating that the triggering of a grand maximum is independent of the time elapsed since the preceding grand maximum.
Our simulation 4 exhibits one robust pattern which has no counterpart in the reconstructions based on radioisotopes: a Grand Maximum is almost always followed by a Grand Minimum. This is again a direct reflection of the manner in which the dynamo enters a Grand Minimum, through a magneticallydriven reduction of differential rotation associate with a high amplitude excursion of the magnetic field.
6.3 Surface differential rotation during Grand Minima and Maxima
Another interesting feature accessible from this simulation, and with a potential observational counterpart, is the variation of surface differential rotation during Grand Minima. Figures 1B and 1C display respectively the total angular velocity at epochs of a Grand Maximum and Minimum in simulation 4. Differences with the steady profile in (A) are hard to distinguish on such a plot, but are present at a significant level at latitudes below ±30° in the outer convection zone. Accordingly, we plot on Figure 13A the surface angular velocity difference between the equator and latitude ±20° at times of Grand Minima and Maxima, as labeled. The corresponding quantity averaged over the whole simulation is indicated by the green asterisk (*). This plot indicates that lowlatitude surface differential rotation is reduced (increased) during Grand Minima (Maxima). Figure 13B displays the latitudinal profiles of this variation for two representative Grand Minima and Maxima. These are obtained by subtracting the latitudinal profile averaged over the whole simulation from the profiles at times of the selected Grand Minima and Maxima. It is noteworthy that these variations, ≈8 nHz peaktopeak, are of the same order of magnitude as inferences based on sunspot drawings during the Maunder Minimum (Casas et al., 2006). This agreement should not be overinterpreted, however, as our internal differential rotation profile shows significant departures from helioseismic inferences, notably the alignment of lowlatitude isocontours with the rotation axis, strong cylindrical shear peaking in the middle of the lowlatitude convection zone rather than at its base, and the lack of a surface shear layer (cf. Fig. 1 herein and Fig. 1 in Howe, 2009).
Fig. 13 (A) Difference in surface angular velocity between the equator and latitude 20° latitude for every Grand Maxima (Minima) in red and yellow (blue and black) measured in simulation, i.e., at extrema of the time series plotted on Figure 12A. The green asterisk indicates the simulation mean for this quantity. (B) Surface differential rotation against latitude for two grand maxima in red and two grand minima in blue, with the simulation mean again subtracted. 
7 Discussion and conclusion
We have developed a nonkinematic axisymmetric meanfield dynamo model of the α^{2}Ω variety, i.e., a mean field model in which both the turbulent electromotive force and shearing by differential rotation contribute to the induction of the toroidal magnetic component. The motivation underlying this work originates with analyses of global MHD simulations of solar convection producing largescale magnetic cycles, which indicate that the turbulent electromotive force characterizing these simulations can be well described by a full αtensor (e.g., Racine et al., 2011; Augustson et al., 2015; Simard et al., 2016; Warnecke et al., 2018). Taken at face value, these inferences from numerical simulations raise doubts on the applicability of the αΩ regime commonly used to model solar and stellar dynamos.
A number of wellknown properties characterizing nonkinematic αΩ meanfield models are found to carry over to the α^{2}Ω regime:
at high magnetic Prandtl number (Pm > 1), nonlinear magneticallymediated backreaction on the differential rotation can produce stable magnetic cycles,
for Pm < 1, modulation of the amplitude of the primarily cycle materializes, with the modulation period increasing with decreasing Pm,
amplitude modulation can take place either through parity changes, or energy exchange with the largescale flow, or a combination of both.
One feature distinguishing our nonkinematic α^{2}Ω model from its αΩ counterparts, however, is that magneticallydriven variations of differential rotation remains quite small, at the level of a few percent or less, independently of the value of Pm and level of supercriticality, at least in the parameter ranges explored. This novel feature demands an explanation. The analyses of Racine et al. (2011) indicates that in the EULAGMHD simulations used here to extract the αtensor, the turbulent contribution to the induction of the mean toroidal field is of opposite sign and almost exactly balances the shearing of the largescale poloidal component by differential rotation; in other words, the two last terms in equation (6) almost cancel each other, with their residual driving the production of the mean toroidal component (see Fig. 14 in Racine et al., 2011). Under such circumstances, even a small change in one of these two contributions – here shearing by differential rotation – can thus lead to a large change in the residual net electromotive force. Consequently, even a very small change in differential rotation can saturate the dynamo, or induce large changes in overall cycle amplitude.
As a mechanisms for generating Grand Minima, our simulations in the Pm < 1 regime do reproduce some features inferred from longterm reconstructions based on the cosmogenic radioisotope records, including in particular an exponential distribution of interevent waiting times (Usoskin et al., 2007), as well as parity changes upon exiting Grand Minima (Ribes & NesmeRibes 1993). The simulation analyzed in detail in Section 6 presents one clearly nonsolar feature, namely the fact that Grand Minima are almost invariably preceeded by Grand Maxima. We are currently pursuing our exploration of the model’s parameter space, to assess whether or not this is an ubiquitous feature of our nonkinematic α^{2}Ω model. Experience with other types of nonlinear dynamo models also suggests that this undesirable characteristic could be eliminated by the introduction of low amplitude stochastic noise (see, e.g., Olemskoy & Kitchatinov 2013).
Acknowledgments
This work was supported by the Natural Sciences and Engineering Research Council of Canada, the Fond Québécois pour la Recherche – Nature et Technologie, the Canadian Foundation for Innovation, and time allocation on the computing infrastructures of Calcul Québec, a member of the Compute Canada consortium. The editor thanks Dmitry Sokoloff and an anonymous reviewer for their assistance in evaluating this paper.
References
 Augustson K, Brun AS, Miesch M, Toomre J. 2015. Grand Minima and Equatorward Propagation in a Cycling Stellar Convective Dynamo. Astrophys J 809: 149. https://doi.org/10.1088/0004637X/809/2/149. [CrossRef] [Google Scholar]
 Beer J, Tobias S, Weiss N. 1998. An active sun throughout the maunder minimum. Solar Phys 181: 237–249. https://doi.org/10.1023/A:1005026001784. [Google Scholar]
 Brandenburg A, Subramanian K. 2005. Astrophysical magnetic fields and nonlinear dynamo theory. Physics Rep 417: 1–209. https://doi.org/10.1016/j.physrep.2005.06.005. [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. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brown BP, Browning MK, Brun AS, Miesch MS, Toomre J. 2010. Persistent magnetic wreaths in a rapidly rotating sun. Astrophys J 711: 424–438. https://doi.org/10.1088/0004637X/711/1/424. [NASA ADS] [CrossRef] [Google Scholar]
 Brown BP, Miesch MS, Browning MK, Brun AS, Toomre J. 2011. Magnetic cycles in a convective dynamo simulation of a young solartype star. Astrophys J 731: 69. https://doi.org/10.1088/0004637X/731/1/69. [CrossRef] [Google Scholar]
 Bushby PJ. 2006. Zonal flows and grand minima in a solar dynamo model. Mon Not Roy Astron Soc 371: 772–780. https://doi.org/10.1111/j.13652966.2006.10706.x. [NASA ADS] [CrossRef] [Google Scholar]
 Cameron RH, Schüssler M. 2017. Understanding solar cycle variability. Astrophys J 843(2): 111. https://doi.org/10.3847/15384357/aa767a. [Google Scholar]
 Casas R, Vaquero JM, Vazquez M. 2006. Solar rotation in the 17th century. Solar Phys 234: 379–392. https://doi.org/10.1007/s1120700600362. [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau P. 2010. Dynamo models of the solar cycle. Living Rev Sol Phys 7: 3. [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau P, MacGregor KB. 2001. Magnetic fields in massive stars. I. Dynamo models. Astrophys J 559: 1094–1107. https://doi.org/10.1086/322417. [NASA ADS] [CrossRef] [Google Scholar]
 Choudhuri AR. 1990. On the possibility of an alphasq omegatype dynamo in a thin layer inside the sun. Astrophys J 355: 733–744. https://doi.org/10.1086/168806. [CrossRef] [Google Scholar]
 Cliver EW, Bounar KH, Boriakoff V. 1998. Geomagnetic activity and the solar wind during the maunder minimum. In: Synoptic solar physics vol. 140 of Astronomical Society of the Pacific Conference Series, Balasubramaniam KS, Harvey J, Rabin D (Eds.), Astronomical Society of the Pacific, San Francisco, pp. 437–445. [Google Scholar]
 Damon PE, Sonett CP. 1991. Solar and terrestrial components of the atmospheric C14 variation spectrum. In: The Sun in time, Sonett CP, Giampapa MS, Matthews MS, (Eds.), University of Arizona Press, Tucson, pp. 360–388. [Google Scholar]
 Duarte LDV, Wicht J, Browning MK, Gastine T. 2016. Helicity inversion in spherical convection as a means for equatorward dynamo wave propagation. Mon Not Roy Astron Soc 456: 1708–1722. https://doi.org/10.1093/mnras/stv2726. [CrossRef] [Google Scholar]
 Eddy JA. 1976. The Maunder minimum. Science 192: 1189–1202. https://doi.org/10.1126/science.192.4245.1189. [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Fan Y, Fang F. 2014. A simulation of convective dynamo in the solar convective envelope: Maintenance of the solarlike differential rotation and emerging flux. Astrophys J 789(1): 35. http://stacks.iop.org/0004637X/789/i=1/a=35. [Google Scholar]
 Ghizaru M, Charbonneau P, Smolarkiewicz PK. 2010. Magnetic cycles in global largeeddy simulations of solar convection. Astrophys J Lett 715: L133–L137. https://doi.org/10.1088/20418205/715/2/L133. [CrossRef] [Google Scholar]
 Guerrero G, Smolarkiewicz PK, de Gouveia Dal Pino EM, Kosovichev AG, Mansour NN. 2016. On the role of tachoclines in solar and stellar dynamos. Astrophys J 819: 104. https://doi.org/10.3847/0004637X/819/2/104. [Google Scholar]
 Guerrero G, Zaire B, Smolarkiewicz PK, de Gouveia Dal Pino EM, Kosovichev AG, Mansour NN. 2019. What sets the magnetic field strength and cycle period in solartype stars? Astrophys J 880(1): 6. https://doi.org/10.3847/15384357/ab224a. [Google Scholar]
 Hotta H, Rempel M, Yokoyama T. 2016. Largescale magnetic fields at high Reynolds numbers in magnetohydrodynamic simulations. Science 351: 1427–1430. https://doi.org/10.1126/science.aad1893. [NASA ADS] [CrossRef] [Google Scholar]
 Howe R. 2009. Solar interior rotation and its variation. Living Rev Sol Phys 6: 1. https://doi.org/10.12942/lrsp20091. [CrossRef] [Google Scholar]
 Inceoglu Fadil, Arlt Rainer, Rempel Matthias. 2017. The nature of grand minima and maxima from fully nonlinear flux transport dynamos. Astrophys J 848(2): 93. https://doi.org/10.3847/15384357/aa8d68. [NASA ADS] [CrossRef] [Google Scholar]
 Käpylä MJ, Käpylä PJ, Olspert N, Brandenburg A, Warnecke J, Karak BB, Pelt J. 2016. Multiple dynamo modes as a mechanism for longterm solar activity variations. A&A 589: A56. https://doi.org/10.1051/00046361/201527002. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Käpylä PJ, Mantere MJ, Brandenburg A. 2012. Cyclic magnetic activity due to turbulent convection in spherical wedge geometry. Astrophys J Lett 755: L22. https://doi.org/10.1088/20418205/755/1/L22. [Google Scholar]
 Käpylä PJ, Mantere MJ, Brandenburg A. 2013. Oscillatory largescale dynamos from Cartesian convection simulations. Geophys Astrophys Fluid Dyn 107: 244–257. https://doi.org/10.1080/03091929.2012.715158. [NASA ADS] [CrossRef] [Google Scholar]
 Kitchatinov LL, Rüdiger G. 2005. Differential rotation and meridional flow in the solar convection zone and beneath. Astron Nachr 326(6): 379–385. https://doi.org/10.1002/asna.200510368. [Google Scholar]
 Kitchatinov LL, Ruediger G, Kueker M. 1994. {LAMBDA}quenching as the nonlinearity in stellarturbulence dynamos. A&A 292: 125–132. [Google Scholar]
 Knobloch E, Tobias SM, Weiss NO. 1998. Modulation and symmetry changes in stellar dynamos. Mon Not Roy Astron Soc 297: 1123–1138. https://doi.org/10.1046/j.13658711.1998.01572.x. [NASA ADS] [CrossRef] [Google Scholar]
 Krause F, Rädler KH. 1980. Meanfield magnetohydrodynamics and dynamo theory, Pergamon Press, Oxford. [Google Scholar]
 Küker M, Arlt R, Rüdiger G. 1999. The Maunder minimum as due to magnetic Lambdaquenching. A&A 343: 977–982. [Google Scholar]
 Malkus WVR, Proctor MRE. 1975. The macrodynamics of αeffect dynamos in rotating fluids. J. Fluid Mech. 67: 417–443. [NASA ADS] [CrossRef] [Google Scholar]
 Masada Y, Yamada K, Kageyama A. 2013. Effects of penetrative convection on solar dynamo. Astrophys J 778: 11. https://doi.org/10.1088/0004637X/778/1/11. [NASA ADS] [CrossRef] [Google Scholar]
 Moffatt HK. 1978. Magnetic field generation in electrically conducting fluids, Cambridge University Press, Cambridge. [Google Scholar]
 Moss D, Brandenburg A, Tuominen I. 1991. Properties of mean field dynamos with nonaxisymmetric alphaeffect. A&A 247: 576–579. [Google Scholar]
 Nelson NJ, Brown BP, Brun AS, Miesch MS, Toomre J. 2013. Magnetic wreaths and cycles in convective dynamos. Astrophys J 762: 73. https://doi.org/10.1088/0004637X/762/2/73. [CrossRef] [Google Scholar]
 Olemskoy SV, Kitchatinov LL. 2013. Grand minima and northsouth asymmetry of solar activity. Astrophys J 777: 71. https://doi.org/10.1088/0004637X/777/1/71. [NASA ADS] [CrossRef] [Google Scholar]
 Ossendrijver M. 2003. The solar dynamo. A&A Rev. 11: 287–367. https://doi.org/10.1007/s0015900300193. [NASA ADS] [CrossRef] [Google Scholar]
 Passos D, Charbonneau P. 2014. Characteristics of magnetic solarlike cycles in a 3D MHD simulation of solar convection. A&A 568: A113. https://doi.org/10.1051/00046361/201423700. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Passos D, Nandy D, Hazra S, Lopes I. 2014. A solar dynamo model driven by meanfield alpha and BabcockLeighton sources: fluctuations, grandminimamaxima, and hemispheric asymmetry in sunspot cycles. A&A 563: A18. https://doi.org/10.1051/00046361/201322635. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Peristykh AN, Damon PE. 2003. Persistence of the Gleissberg 88year solar cycle over the last 12,000 years: Evidence from cosmogenic isotopes. J Geophys Res (Space Phys) 108: 1003. https://doi.org/10.1029/2002JA009390. [NASA ADS] [CrossRef] [Google Scholar]
 Racine É, Charbonneau P, Ghizaru M, Bouchat A, Smolarkiewicz PK. 2011. On the mode of dynamo action in a global largeeddy simulation of solar convection. Astrophys J 735: 46. https://doi.org/10.1088/0004637X/735/1/46. [Google Scholar]
 Rädler KH. 1980. Meanfield approach to spherical dynamo models. Astron Nachr 301: 101–129. [NASA ADS] [CrossRef] [Google Scholar]
 Rempel M. 2006. Fluxtransport dynamos with lorentz force feedback on differential rotation and meridional flow: saturation mechanism and torsional oscillations. Astrophys J 647(1): 662–675. https://doi.org/10.1086/505170. [Google Scholar]
 Ribes JC, NesmeRibes E. 1993. The solar sunspot cycle in the Maunder minimum AD1645 to AD1715. A&A 276: 549. [Google Scholar]
 Ruediger G, Brandenburg A. 1995. A solar dynamo in the overshoot layer: cycle period and butterfly diagram. A&A 296: 557. [Google Scholar]
 Schrinner M, Rädler KH, Schmitt D, Rheinhardt M, Christensen UR. 2007. Meanfield concept and direct numerical simulations of rotating magnetoconvection and the geodynamo. Geophys Astrophys Fluid Dyn 101(2): 81–116. https://doi.org/10.1080/03091920701345707. [CrossRef] [Google Scholar]
 Simard C, Charbonneau P, Bouchat A. 2013. Magnetohydrodynamic simulationdriven kinematic mean field model of the solar cycle. Astrophys J 768: 16. https://doi.org/10.1088/0004637X/768/1/16. [NASA ADS] [CrossRef] [Google Scholar]
 Simard C, Charbonneau P, Dubé C. 2016. Characterisation of the turbulent electromotive force and its magneticallymediated quenching in a global EULAGMHD simulation of solar convection. Adv Space Res 58: 1522–1537. https://doi.org/10.1016/j.asr.2016.03.041. [NASA ADS] [CrossRef] [Google Scholar]
 Simitev RD, Kosovichev AG, Busse FH. 2015. Dynamo effects near the transition from solar to antisolar differential rotation. Astrophys J 810: 80. https://doi.org/10.1088/0004637X/810/1/80. [NASA ADS] [CrossRef] [Google Scholar]
 Sokoloff D, NesmeRibes E. 1994. The Maunder minimum: A mixedparity dynamo mode? A&A 288: 293–298. [Google Scholar]
 Strugarek A, Beaudoin P, Charbonneau P, Brun AS. 2018. On the sensitivity of magnetic cycles in global simulations of solarlike stars. Astrophys J 863(1): 35. https://doi.org/10.3847/15384357/aacf9e. [Google Scholar]
 Tobias SM. 1996. Grand minimia in nonlinear dynamos. A&A 307: L21. [Google Scholar]
 Tobias SM. 1997. The solar cycle: parity interactions and amplitude modulation. A&A 322: 1007–1017. [Google Scholar]
 Usoskin IG. 2017. A history of solar activity over millennia. Living Rev Sol Phys 14(1): 3. https://doi.org/10.1007/s4111601700069. [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin IG, Kovaltsov GA. 2004. Longterm solar activity: direct and indirect study. Solar Phys 224: 37–47. https://doi.org/10.1007/s1120700539977. [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin IG, Solanki SK, Kovaltsov GA. 2007. Grand minima and maxima of solar activity: new observational constraints. A&A 471: 301–309. https://doi.org/10.1051/00046361:20077704. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin IG, Sokoloff D, Moss D. 2009. Grand minima of solar activity and the meanfield dynamo. Solar Phys 254: 345–355. https://doi.org/10.1007/s1120700892936. [NASA ADS] [CrossRef] [Google Scholar]
 Usoskin IG, Arlt R, Asvestari E, Hawkins E, Käpylä M, et al. 2015. The Maunder minimum (1645–1715) was indeed a grand minimum: A reassessment of multiple datasets. A&A 581: A95. https://doi.org/10.1051/00046361/201526652. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Usoskin IG, Gallet Y, Lopes F, Kovaltsov GA, Hulot G. 2016. Solar activity during the Holocene: the Hallstatt cycle and its consequence for grand minima and maxima. A&A 587: A150. https://doi.org/10.1051/00046361/201527295. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vasiliev SS, Dergachev VA. 2002. The 2400year cycle in atmospheric radiocarbon concentration: bispectrum of 14C data over the last 8000 years. Ann Geophys 20: 115–120. https://doi.org/10.5194/angeo201152002. [NASA ADS] [CrossRef] [Google Scholar]
 Warnecke J, Rheinhardt M, Tuomisto S, Käpylä PJ, Käpylä MJ, Brandenburg A. 2018. Turbulent transport coefficients in spherical wedge dynamo simulations of solarlike stars. A&A 609: A51. https://doi.org/10.1051/00046361/201628136. [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cite this article as: Simard C & Charbonneau P. 2020. Grand Minima in a spherical nonkinematic α^{2}Ω meanfield dynamo model. J. Space Weather Space Clim. 10, 9.
All Tables
Defining dimensionless parameters for the reference dynamo solutions discussed in the text.
All Figures
Fig. 1 (A) Steady angular velocity component (see Eq. (8)) as extracted from a hydrodynamical EULAG simulation, plotted as isocontours of rotational frequency. (B) and (C) Total angular velocity at epoch of Grand Maximum and Minimum, respectively, as produced by our nonkinematic α^{2}Ω meanfield dynamo model (specifically, simulation 4 in Table 1). The interface between the convection zone and underlying radiative zone is indicated by the white dashed line. 

In the text 
Fig. 2 A kinematic α^{2}Ω dynamo solution. Defining parameters are listed in Table 1 (solution run 0). The toroidal field component is plotted as a timelatitude diagram at depth r/R = 0.85 in the top panel, and meridional plane snapshots extracted at times corresponding to the dashed lines on the top panel, and spanning a magnetic halfcycle, are displayed at bottom. 

In the text 
Fig. 3 Solution run 1, a nonkinematic dynamo solution in the Pm > 1 regime. Its defining parameters are listed in the second line of Table 1. (A) Time series of normalized magnetic energy (black), integrated shear proxy (red), and equatorial parity (blue), starting from a seed magnetic field and Ω′ = 0 at t = 0. (B) Timelatitude diagram of the toroidal magnetic component extracted at r/R = 0.85, with a meridional plane snapshot at right. (C) Similar timelatitude diagram at r/R = 0.85 and meridional snapshot, but for the angular velocity perturbation normalized to the background angular velocity profile, i.e., Ω′(r, θ, t)/Ω(r, θ). 

In the text 
Fig. 4 Simulation 2, a nonkinematic α^{2}Ω dynamo exhibiting parity modulation. (A) Normalized magnetic energy in black, the shear proxy (Ψ) integrated in the meridional plane in red, and the parity in light blue. (B) Toroidal magnetic field in a timelatitude diagram extracted at 0.85R and (C) the angular velocity fluctuation (Ω′(r, θ, t)/Ω(r, θ)) in the same configuration. The meridional cuts of the toroidal field B and angular velocity perturbation Ω′ plotted at right of the timelatitude diagrams are extracted at t/τ = 2.7, as indicated by the dashdotted line segments on the corresponding timelatitude diagrams. 

In the text 
Fig. 5 Top panel: 3D phase portrait between parity, toroidal magnetic field and the shear proxy (see Eq. (14)) from solution 2. Each of these three quantities is extracted at r/R = 0.85 and 15° latitude. Only the modulation envelope of those quantities is used for plotting, i.e., the primary cycle is removed. The color code represents the magnitude of the parity from −1 in blue to +1 in red. The bottom panel is a Poincaré section in the P = 0 plane, with the inset focusing on one of the four crossing locations. The trajectory therein is spacefilling, indicating that the modulation is chaotic (see text). 

In the text 
Fig. 6 Identical in format to Figure 4, but showing now a very low Pm solution (Pm = 5 × 10^{−3}; solution run 3 in Table 1). 

In the text 
Fig. 7 Variation of the ratio between the longterm modulation period (P_{2}) and the primary magnetic cycle period (P_{1}) versus the magnetic Prandtl number (Pm). Except for the latter, all solutions use the same parameter values as solution 3 in Table 1 (see also Fig. 6). The dashed line indicates a Pm^{−1} scaling. 

In the text 
Fig. 8 Identical in format to Figure 4, but for a solution exhibiting strongly aperiodic behavior. This is solution run 4 in Table 1. 

In the text 
Fig. 9 Snapshot of an epoch of a grand maximum followed by a grand minimum for solution run 4. Top panel shows the toroidal magnetic field in a timelatitude diagram extracted at r/R = 0.85. The bottom panel shows four meridional planes snapshots extracted from the simulation at times indicated by the dashed lines on the upper panel. The color scale encodes the toroidal magnetic field, and poloidal fieldlines are plotted as solid (dashed) lines for counterclockwise (clockwise) orientation. 

In the text 
Fig. 10 3D phase portrait between parity, toroidal magnetic field and shear proxy from solution 4 taken at the same time span as Figure 9 and extracted at 0.85R radius and 15° latitude. The color encodes the advance of time from blue to red. 

In the text 
Fig. 11 (A) Amplitude distribution of the toroidal magnetic field maxima extracted from the primary cycle in simulation run 4. Panel (B) separates these amplitude by parity, considering only cycles for which P ≥ 0.7. The corresponding distribution for symmetric and antisymmetric parity are plotted in yellow and purple, respectively, with their sum in red. 

In the text 
Fig. 12 (A) Time series of total magnetic energy in a 11,000 years segment of a 50,000 years realization of simulation 4. Time is expressed both in units diffusion time (bottom axis) and in years (top axis). The oscillations from the primary cycle have been removed, leaving only the longterm modulation envelope. Epochs of grand maxima and minima are indicated in red and blue, respectively. Compare to Figure 3 in Usoskin et al. (2007). The vertical dashed lines delineate the time interval spanned by Figure 8. The four other panels show histograms of (B) Grand Minima durations, (C) waiting times between successive Grand Minima, (D) Grand Maxima durations, and (E) waiting times between successive Grand Maxima (see text). 

In the text 
Fig. 13 (A) Difference in surface angular velocity between the equator and latitude 20° latitude for every Grand Maxima (Minima) in red and yellow (blue and black) measured in simulation, i.e., at extrema of the time series plotted on Figure 12A. The green asterisk indicates the simulation mean for this quantity. (B) Surface differential rotation against latitude for two grand maxima in red and two grand minima in blue, with the simulation mean again subtracted. 

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.