© C. Simard & P. Charbonneau, Published by EDP Sciences 2020

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 well-organized on large spatial scales, and undergoing more-or-less 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 dynamically-simplified dynamo models remain the most practical approach.

Solar cycle models based on mean-field electrodynamics have indeed been used for many years to study long-term cycle variability (see Sect. 5 in Charbonneau, 2010, and references therein). Most often such models solve only for the large-scale, axisymmetric component of the solar magnetic field, with the effects of small-scale, 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 large-scale 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 Babcock-Leighton mechanism of dipole regeneration through the surface decay of active regions.

The extraction of the mean-field 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 small-scales, parameter regimes, and the techniques used to extract the mean-field coefficients also differ. While results differ at various levels, some common features nonetheless emerge. Most notably perhaps, the α-tensor is full, with off-diagonal components of magnitude comparable to diagonal elements. The antisymmetric component of the tensor is also of comparable magnitude to the off-diagonal symmetric components, indicating that turbulent pumping represents a significant transport mechanism for the large-scale magnetic field.

In contrast, mean-field dynamo models of the solar cycle are often designed in the so-called αΩ regime: the poloidal large-scale 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 order-or-magnitude 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 α²Ω regime: both shearing by differential rotation and turbulent inductive effects contribute to the production of the large-scale toroidal magnetic component. Even in the kinematic regime, in which large-scale flows are assumed steady, α²Ω mean-field models have received far less attention than their αΩ counterparts. Some well-known solar-like properties of αΩ models do carry over to the α²Ω 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 α²Ω context is the dynamical backreaction of the large-scale magnetic field produced by dynamo action on the large-scale inductive flows, specifically the differential rotation. This so-called Malkus-Proctor 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 Minimum-like 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 large-scale magnetic field then dissipates, eventually allowing the re-establishment 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 α²Ω 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 non-kinematic mean-field model are described in Section 2. In Section 3, we present a linear, kinematic reference dynamo solution. Representative high- and low-Prandtl number (Pm) non-kinematic dynamo solutions are then presented and analyzed in Sections 4 and 5, respectively. Section 6 focuses on the long timescale modulations materializing in the low-Pm regime. We close in Section 7 by examining the similarity and differences with non-kinematic mean-field models of the αΩ variety, and seek the origin of the most striking differences in the interactions between the various inductive effects operating in our α²Ω model.

2 A non-kinematic mean-field model

In this section, we briefly summarize the mathematical development of our non-kinematic mean-field model. Our formulation of α²Ω axisymmetric mean-field 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 Mean-field dynamo equation

Following the usual mean-fied 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 mean-field 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 small-scale, 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 small-scale 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 large-scale magnetic and flow field on the other. The first term in the expansion involves a rank-two tensor capturing (among other effects) the so-called α-effect, which arises from the symmetric part of the a tensor and can act as a source term in the mean-field equations (5) and (6) (see Rädler, 1980; Schrinner et al., 2007; Simard et al., 2016).

The axisymmetric form of the mean-field 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 mean-field 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 large-scale flow in equation (1), i.e., we do not include any large-scale meridional circulation in the model. Note that here the mean electromotive force ℰ is retained in equation (6), which yields the so-called α²Ω mean-field 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²/η₀. This leads to the appearance of two dimensionless numbers measuring the inductive strength of the mean-electromotive force and large-scale shear, respectively:(7)where α₀, η₀, and Ω₀ 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 large-scale flows. Since the Sun’s convection zone is a strongly turbulent environment, this dynamical backreaction can operate at all spatial scales. Suppression of the small-scale 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 non-kinematic dynamo model, we only include the backreaction of the mean magnetic field on the large-scale flow, also know as the Malkus-Proctor 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 time-dependent 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 large-scale magnetic field produced by the dynamo, so that Ω represents a stationary solution of the unmagnetized momentum equation. The governing equation for the time-varying 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 non-kinematic α²Ω mean-field dynamo model

Our non-kinematic mean-field 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 element-based 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 EUlerian-LAGrangian (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 solar-like features, including a rapidly rotating equatorial region and slowly rotating poles, and a thin tachocline-like rotational shear layer immediately beneath the core-envelope interface, indicated in Figure 1A by the white dashed line. Non-solar features at low latitudes include alignment of the isocontours with the rotation axis, low-latitude radial shear peaking within the convective envelope rather than at its base, and lack of a near-surface shear layer. The large-scale meridional flow generated in this simulation is weak, shows strong temporal fluctuations, and lacks a well-defined spatial organization (Racine et al., 2011). We therefore opt to exclude it from our α²Ω mean-field 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 non-kinematic α2Ω mean-field 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 mean-field dynamo model use an α-tensor extracted from a reference EULAG-MHD simulation of solar convection (Passos & Charbonneau, 2014) through a least-square 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 off-diagonal 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¹¹ cm² s⁻¹, commensurate with the values extracted from our reference EULAG-MHD 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 core-envelope 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 core-to-envelope diffusivity ratio of 10⁻².

The density ρ(r) in equation (10) follows the same polytropic profile as used in the EULAG-MHD from which the α-tensor and mean differential rotation are extracted. Likewise, the core-envelope interface is set at r/R = 0.718. We set a floor on the density at ρ = 0.021 kg m⁻³ 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 solar-looking 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 non-kinematic α²Ω mean-field 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 EULAG-MHD millenium simulation, we emphasize that our primary aim here is not to produce a close mean-field equivalent, but rather to explore the range of possible dynamical behavior associated with magnetic backreaction on differential rotation in a mean-field α²Ω dynamo formulated in spherical geometry. Table 1 lists the defining parameters of the set of representative simulations discussed in the remainder of the paper.

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 time-latitude diagram of the toroidal component extracted at mid-depth (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 time-latitude 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 half-cycle, 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 αΩ mean-field models and also materializing in α²Ω 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 Parker-Yoshimura 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., dipole-like).

4 Simulation 1: amplitude saturation

We now turn to a non-kinematic 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 equatorially-propagating 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 non-kinematic 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) Time-latitude diagram of the toroidal magnetic component extracted at r/R = 0.85, with a meridional plane snapshot at right. (C) Similar time-latitude 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 time-latitude diagram of the differential rotation perturbation Ω′/Ω, extracted at the same depth r/R = 0.85 as the time-latitude 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 large-scale flow.

5 Simulations 2, 3, and 4: amplitude modulation

While our non-kinematic α²Ω dynamo solutions stabilize on a constant-amplitude 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 non-kinematic αΩ dynamo models, in which the turbulent electromotive force is zeroed out in equation (6). Geometrically, these models have ranged from spatially-integrated dynamical systems, through one-dimensional 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 non-kinematic 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 co-existence 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 mean-field models of the αΩ variety. We now turn to non-kinematic α²Ω 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 non-kinematic α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 time-latitude 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 time-latitude diagrams are extracted at t/τ = 2.7, as indicated by the dash-dotted line segments on the corresponding time-latitude diagrams.

As can be seen on Figure 4B, parity modulation leads here to an alternating pattern of hemispheric dominance in the time-latitude 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 high-Pm 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 magnetically-driven 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 three-dimensional 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 well-defined 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 space-filling, 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 space-filling, indicating that the modulation is chaotic (see text).

Figure 6 shows yet another example of a non-kinematic α²Ω dynamo solution (run 3 from Table 1), this time obtained at a much lower value of the magnetic Prandtl number, Pm = 5 × 10⁻³. Magnetic activity is now peaking at mid-latitudes in both hemisphere, with much less trans-equatorial 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 well-known of non-kinematic αΩ mean-field models operating in the low-Pm 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 long-term modulation period (P₂) and the primary cycle period (P₁) is observed to increase when Pm is decreased (Tobias, 1997; Bushby, 2006). A similarly robust trend also emerges from our non-kinematic α²Ω 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₂/P₁ is found to scale as ∝Pm⁻¹ over an order-of-magnitude in Pm.

Fig. 7 Variation of the ratio between the long-term modulation period (P2) and the primary magnetic cycle period (P1) 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 long-term 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 non-kinematic α²Ω 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 non-kinematic αΩ 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 400-year 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 ¹⁴C or ¹⁰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 Minimum-like 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, mean-field dynamo models remain the primary modeling tool used to explore the underlying physical mechanisms.

Many dynamo-based 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 solar-like 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 sunspot-forming 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 ¹⁰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 & Nesme-Ribes, 1993; Sokoloff & Nesme-Ribes, 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 time-latitude 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 time-latitude 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 well-defined oscillatory path in this three-dimensional 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 gaussian-like fashion. However, if one extracts cycles for which the parity is either predominantly dipole-like or quadrupole-like, 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 year-long 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 radioisotope-based 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 long-term activity modulation trend remains. The fact that this simulation exhibits markedly aperiodic long-term 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 long-term 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 well-defined 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 well-defined 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 inter-event 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 magnetically-driven 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 low-latitude 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 peak-to-peak, 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 low-latitude isocontours with the rotation axis, strong cylindrical shear peaking in the middle of the low-latitude 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 non-kinematic axisymmetric mean-field dynamo model of the α²Ω 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 large-scale 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 well-known properties characterizing non-kinematic αΩ mean-field models are found to carry over to the α²Ω regime:

1. at high magnetic Prandtl number (Pm > 1), nonlinear magnetically-mediated backreaction on the differential rotation can produce stable magnetic cycles,

2. for Pm < 1, modulation of the amplitude of the primarily cycle materializes, with the modulation period increasing with decreasing Pm,

3. amplitude modulation can take place either through parity changes, or energy exchange with the large-scale flow, or a combination of both.

One feature distinguishing our non-kinematic α²Ω model from its αΩ counterparts, however, is that magnetically-driven 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 EULAG-MHD 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 large-scale 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 long-term reconstructions based on the cosmogenic radioisotope records, including in particular an exponential distribution of inter-event waiting times (Usoskin et al., 2007), as well as parity changes upon exiting Grand Minima (Ribes & Nesme-Ribes 1993). The simulation analyzed in detail in Section 6 presents one clearly non-solar 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 non-kinematic α²Ω 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

All Tables

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 non-kinematic α2Ω mean-field 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 time-latitude 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 half-cycle, are displayed at bottom.
In the text

Fig. 3 Solution run 1, a non-kinematic 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) Time-latitude diagram of the toroidal magnetic component extracted at r/R = 0.85, with a meridional plane snapshot at right. (C) Similar time-latitude 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 non-kinematic α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 time-latitude 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 time-latitude diagrams are extracted at t/τ = 2.7, as indicated by the dash-dotted line segments on the corresponding time-latitude 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 space-filling, 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 long-term modulation period (P2) and the primary magnetic cycle period (P1) 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 time-latitude 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 long-term 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